Source code for wradlib.verify

#!/usr/bin/env python
# Copyright (c) 2011-2023, wradlib developers.
# Distributed under the MIT License. See LICENSE.txt for more info.


Verification mainly refers to the comparison of radar-based precipitation
estimates to ground truth.

.. autosummary::
   :toctree: generated/

__all__ = ["ErrorMetrics", "PolarNeighbours"]
__doc__ = __doc__.format("\n   ".join(__all__))

from pprint import pprint

import numpy as np
from scipy import spatial, stats

from wradlib import util
from wradlib.georef import polar

[docs]class PolarNeighbours: """For a set of projected point coordinates, extract the neighbouring bin \ values from a data set in polar coordinates. Use as follows: First, create an instance of PolarNeighbours by passing all the information needed to georeference the polar radar data to the points of interest (see parameters). Second, use the method *extract* in order to extract the values from a data array which corresponds to the polar coordinates. Parameters ---------- r : :class:`numpy:numpy.ndarray` (see :mod:`wradlib.georef` for documentation) az : :class:`numpy:numpy.ndarray` (see :mod:`wradlib.georef` for documentation) site : sequence sequence of floats (see :mod:`wradlib.georef` for documentation) crs : :py:class:`gdal:osgeo.osr.SpatialReference` GDAL OSR Spatial Reference Object describing projection (see :mod:`wradlib.georef` for documentation) x : :class:`numpy:numpy.ndarray` array of x coordinates of the points in map projection corresponding to crs y : :class:`numpy:numpy.ndarray` array of y coordinates of the points in map projection corresponding to crs nnear : int number of neighbouring radar bins you would like to find Examples -------- See :ref:`/notebooks/verification/verification.ipynb`. """ def __init__(self, r, az, site, crs, x, y, *, nnear=9): self.nnear = nnear = az self.r = r self.x = x self.y = y # compute the centroid coordinates in crs bin_coords = polar.spherical_to_centroids(r, az, 0, site, crs=crs) self.binx = bin_coords[..., 0].ravel() self.biny = bin_coords[..., 1].ravel() # compute the KDTree tree = spatial.KDTree(list(zip(self.binx, self.biny))) # query the tree for nearest neighbours self.dist, self.ix = tree.query(list(zip(x, y)), k=nnear)
[docs] def extract(self, vals): """Extracts the values from an array of shape (azimuth angles, \ range gages) which correspond to the indices computed during \ initialisation Parameters ---------- vals : :class:`numpy:numpy.ndarray` array of shape (..., number of azimuth, number of range gates) Returns ------- output : :class:`numpy:numpy.ndarray` array of shape (..., number of points, nnear) """ if vals.ndim < 2: raise ValueError( "`vals` array should at least contain an " "azimuth and a range dimension." ) if tuple(vals.shape[-2:]) != (len(, len(self.r)): raise ValueError( "The shape of `vals` array does not correspond with " "the range and azimuths provided." ) vals = vals.reshape(vals.shape[:-2] + (len( * len(self.r),)) return vals[..., self.ix]
[docs] def get_bincoords(self): """Returns all bin coordinates in map projection Returns ------- output : tuple array of x coordinates, array of y coordinates """ return self.binx, self.biny
[docs] def get_bincoords_at_points(self): """Returns bin coordinates only in the neighbourhood of points Returns ------- output : tuple array of x coordinates, array of y coordinates """ return self.binx[self.ix], self.biny[self.ix]
[docs]class ErrorMetrics: """Compute quality metrics from a set of observations (``obs``) and \ estimates (``est``). First create an instance of the class using the set of observations and estimates. Then compute quality metrics using the class methods. A dictionary of all available quality metrics is returned using the ``all`` method, or printed to the screen using the ``pprint`` method. The ``ix`` member variable indicates valid pairs of ``obs`` and ``est``, based on NaNs and ``minval``. Parameters ---------- obs: :class:`numpy:numpy.ndarray` array of observations (e.g. rain gage observations) est: :class:`numpy:numpy.ndarray` array of estimates (e.g. radar, adjusted radar, ...) minval : float threshold value in order to compute metrics only for values larger than minval Examples -------- >>> obs = np.random.uniform(0, 10, 100) >>> est = np.random.uniform(0, 10, 100) >>> metrics = ErrorMetrics(obs, est) >>> metrics.all() #doctest: +SKIP >>> metrics.pprint() #doctest: +SKIP >>> metrics.ix #doctest: +SKIP See :ref:`/notebooks/verification/verification.ipynb` and :ref:`/notebooks/multisensor/gauge_adjustment.ipynb`. """ def __init__(self, obs, est, *, minval=None): # Check input if len(obs) != len(est): raise ValueError( f"`obs` ({len(obs)}) and `est` ({len(est)}) need to have the same length." ) self.est = est self.obs = obs # remember those pairs which both have valid obs and est self.ix = np.intersect1d( util._idvalid(obs, minval=minval), util._idvalid(est, minval=minval) ) self.n = len(self.ix) if self.n == 0: util.warn("No valid pairs of `obs` and `est` available for ErrorMetrics!") self.resids = self.est[self.ix] - self.obs[self.ix]
[docs] def corr(self): """Correlation coefficient""" return np.round(np.corrcoef(self.obs[self.ix], self.est[self.ix])[0, 1], 2)
[docs] def r2(self): """Coefficient of determination""" return np.round( (np.corrcoef(self.obs[self.ix], self.est[self.ix])[0, 1]) ** 2, 2 )
[docs] def spearman(self): """Spearman rank correlation coefficient""" return np.round(stats.spearmanr(self.obs[self.ix], self.est[self.ix])[0], 2)
[docs] def nash(self): """Nash-Sutcliffe Efficiency""" return np.round(1.0 - (self.mse() / np.var(self.obs[self.ix])), 2)
[docs] def sse(self): """Sum of Squared Errors""" return np.round(np.sum(self.resids**2), 2)
[docs] def mse(self): """Mean Squared Error""" return np.round(self.sse() / self.n, 2)
[docs] def rmse(self): """Root Mean Squared Error""" return np.round(self.mse() ** 0.5, 2)
[docs] def mas(self): """Mean Absolute Error""" return np.round(np.mean(np.abs(self.resids)), 2)
[docs] def meanerr(self): """Mean Error""" return np.round(np.mean(self.resids), 2)
[docs] def ratio(self): """Mean ratio between observed and estimated""" return np.round(np.mean(self.est[self.ix] / self.obs[self.ix]), 2)
[docs] def pbias(self): """Percent bias""" return np.round(self.meanerr() * 100.0 / np.mean(self.obs[self.ix]), 1)
[docs] def all(self): """Returns a dictionary of all error metrics""" out = { "corr": self.corr(), "r2": self.r2(), "spearman": self.spearman(), "nash": self.nash(), "sse": self.sse(), "mse": self.mse(), "rmse": self.rmse(), "mas": self.mas(), "meanerr": self.meanerr(), "ratio": self.ratio(), "pbias": self.pbias(), } return out
[docs] def pprint(self): """Pretty prints a summary of error metrics""" pprint(self.all())
if __name__ == "__main__": print("wradlib: Calling module <verify> as main...")