Recipe 6: Zonal Statistics - Polar Grid

Recipe 6: Zonal Statistics - Polar Grid#

import datetime as dt
import warnings

import matplotlib.pyplot as plt
import numpy as np
import wradlib as wrl
import wradlib_data
import xarray as xr
from IPython.display import display
from matplotlib.colors import from_levels_and_colors
from osgeo import osr

warnings.filterwarnings("ignore")

Setup Examples#

def testplot(
    ds,
    obj,
    col="mean",
    levels=[0, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 40, 50, 100],
    title="",
):
    """Quick test plot layout for this example file"""
    colors = plt.cm.viridis(np.linspace(0, 1, len(levels)))
    mycmap, mynorm = from_levels_and_colors(levels, colors, extend="max")

    radolevels = [0, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 40, 50, 100]
    radocolors = plt.cm.viridis(np.linspace(0, 1, len(radolevels)))
    radocmap, radonorm = from_levels_and_colors(radolevels, radocolors, extend="max")

    fig = plt.figure(figsize=(10, 16))

    # Average rainfall sum
    ax = fig.add_subplot(211, aspect="equal")
    obj.zdata.trg.geo.plot(
        column=col,
        ax=ax,
        cmap=mycmap,
        norm=mynorm,
        edgecolor="white",
        lw=0.5,
        legend=True,
        legend_kwds=dict(shrink=0.5),
    )
    ax.autoscale()
    plt.xlabel("UTM Zone 32 Easting")
    plt.ylabel("UTM Zone 32 Northing")
    plt.title(title)
    plt.draw()

    # Original radar data
    ax1 = fig.add_subplot(212, aspect="equal")
    ds.plot(
        x="xc",
        y="yc",
        cmap=radocmap,
        norm=radonorm,
        ax=ax1,
        cbar_kwargs=dict(shrink=0.5),
    )
    obj.zdata.trg.geo.plot(ax=ax1, facecolor="None", edgecolor="white")
    plt.xlabel("UTM Zone 32 Easting")
    plt.ylabel("UTM Zone 32 Northing")
    plt.title("Original radar rain sums")
    plt.draw()
    plt.tight_layout()

# check for GEOS enabled GDAL
if not wrl.util.has_geos():
    print("NO GEOS support within GDAL, aborting...")
    exit(0)

# create radolan projection osr object
proj_stereo = wrl.georef.create_osr("dwd-radolan")

# create UTM Zone 32 projection osr object
proj_utm = osr.SpatialReference()
proj_utm.ImportFromEPSG(32632)

# Source projection of the shape data (in GK2)
proj_gk2 = osr.SpatialReference()
proj_gk2.ImportFromEPSG(31466)

0

def create_center_coords(ds, crs=None):
    # create polar grid centroids in GK2
    center = wrl.georef.spherical_to_centroids(
        ds.range.values,
        ds.azimuth.values,
        0.5,
        (ds.longitude.values, ds.latitude.values, ds.altitude.values),
        crs=crs,
    )
    ds = ds.assign_coords(
        {
            "xc": (["azimuth", "range"], center[..., 0]),
            "yc": (["azimuth", "range"], center[..., 1]),
            "zc": (["azimuth", "range"], center[..., 2]),
        }
    )
    return ds

filename = wradlib_data.DATASETS.fetch("hdf5/rainsum_boxpol_20140609.h5")
ds = xr.open_dataset(filename)
ds = ds.rename_dims({"phony_dim_0": "azimuth", "phony_dim_1": "range"})
ds = ds.assign_coords(
    {
        "latitude": ds.data.Latitude,
        "longitude": ds.data.Longitude,
        "altitude": 99.5,
        "azimuth": ds.data.az,
        "range": ds.data.r,
        "sweep_mode": "azimuth_surveillance",
        "elevation": 0.5,
    }
)

ds = ds.wrl.georef.georeference(crs=proj_utm)
ds = ds.pipe(create_center_coords, crs=proj_utm)
display(ds)

Downloading file 'hdf5/rainsum_boxpol_20140609.h5' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/rainsum_boxpol_20140609.h5' to '/home/docs/.cache/wradlib-data'.

<xarray.Dataset> Size: 29MB
Dimensions:     (azimuth: 360, range: 1000)
Coordinates: (12/17)
  * azimuth     (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
  * range       (range) float64 8kB 100.0 200.0 300.0 ... 9.99e+04 1e+05
    x           (azimuth, range) float64 3MB 3.639e+05 3.639e+05 ... 3.657e+05
    y           (azimuth, range) float64 3MB 5.622e+06 5.622e+06 ... 5.722e+06
    z           (azimuth, range) float64 3MB 100.4 101.2 ... 1.559e+03 1.561e+03
    gr          (azimuth, range) float64 3MB 99.98 200.0 ... 9.986e+04 9.996e+04
    ...          ...
    latitude    float64 8B 50.73
    longitude   float64 8B 7.072
    altitude    float64 8B 99.5
    sweep_mode  <U20 80B 'azimuth_surveillance'
    elevation   float64 8B 0.5
    crs_wkt     int64 8B 0
Data variables:
    data        (azimuth, range) float64 3MB ...
def write_prj(filename, proj):
    with open(filename, "w") as f:
        f.write(proj.ExportToWkt())

flist = ["shapefiles/agger/agger_merge.shx", "shapefiles/agger/agger_merge.dbf"]
[wradlib_data.DATASETS.fetch(f) for f in flist]
# reshape
shpfile = wradlib_data.DATASETS.fetch("shapefiles/agger/agger_merge.shp")
write_prj(shpfile[:-3] + "prj", proj_gk2)
trg = wrl.io.VectorSource(shpfile, trg_crs=proj_utm, name="trg")

bbox = trg.extent

# create catchment bounding box
buffer = 5000.0
bbox = dict(
    left=bbox[0] - buffer,
    right=bbox[1] + buffer,
    bottom=bbox[2] - buffer,
    top=bbox[3] + buffer,
)

Warning 1: Layer trg has no spatial index, DROP SPATIAL INDEX failed.

ds_clip = ds.where(
    (
        ((ds.yc > bbox["bottom"]) & (ds.yc < bbox["top"]))
        & ((ds.xc > bbox["left"]) & (ds.xc < bbox["right"]))
    ),
    drop=True,
)
display(ds_clip)

<xarray.Dataset> Size: 5MB
Dimensions:     (azimuth: 86, range: 695)
Coordinates: (12/17)
  * azimuth     (azimuth) float64 688B 0.5 1.5 2.5 3.5 ... 82.5 83.5 84.5 85.5
  * range       (range) float64 6kB 2.8e+03 2.9e+03 3e+03 ... 7.21e+04 7.22e+04
    x           (azimuth, range) float64 478kB 3.64e+05 3.64e+05 ... 4.36e+05
    y           (azimuth, range) float64 478kB 5.624e+06 5.625e+06 ... 5.625e+06
    z           (azimuth, range) float64 478kB 124.4 125.3 ... 1.037e+03
    gr          (azimuth, range) float64 478kB 2.799e+03 2.899e+03 ... 7.217e+04
    ...          ...
    latitude    float64 8B 50.73
    longitude   float64 8B 7.072
    altitude    float64 8B 99.5
    sweep_mode  <U20 80B 'azimuth_surveillance'
    elevation   float64 8B 0.5
    crs_wkt     int64 8B 0
Data variables:
    data        (azimuth, range) float64 478kB nan nan nan nan ... nan nan nan
radar_utmc = np.dstack([ds_clip.xc, ds_clip.yc]).reshape(-1, 2)
radar_utmc.shape

(59770, 2)

Zonal Stats Polar Grid - Points#

trg = wrl.io.VectorSource(shpfile, trg_crs=proj_utm, name="trg", src_crs=proj_gk2)
src = wrl.io.VectorSource(radar_utmc, trg_crs=proj_utm, name="src")

Warning 1: Layer trg has no spatial index, DROP SPATIAL INDEX failed.

Warning 1: Layer src has no spatial index, DROP SPATIAL INDEX failed.

###########################################################################
# Approach #1: Assign grid points to each polygon and compute the average.
#
# - Uses matplotlib.path.Path
# - Each point is weighted equally (assumption: polygon >> grid cell)
# - this is quick, but theoretically dirty
# - for polar grids a range-area dependency has to be taken into account
###########################################################################

t1 = dt.datetime.now()

# Create instance of type ZonalDataPoint from source grid and
# catchment array
zd = wrl.zonalstats.ZonalDataPoint(src, trg=trg, crs=proj_utm, buf=500.0)
# dump to file
zd.dump_vector("test_zonal_points")
# Create instance of type ZonalStatsPoint from zonal data object
obj1 = wrl.zonalstats.ZonalStatsPoint(zd)

isecs1 = obj1.zdata.isecs
t2 = dt.datetime.now()

t3 = dt.datetime.now()

# Create instance of type ZonalStatsPoint from zonal data file
obj1 = wrl.zonalstats.ZonalStatsPoint("test_zonal_points")
# Compute stats for target polygons
avg1 = obj1.mean(ds_clip.data.values.ravel())
var1 = obj1.var(ds_clip.data.values.ravel())

t4 = dt.datetime.now()

print("Approach #1 computation time:")
print("\tCreate object from scratch: %f seconds" % (t2 - t1).total_seconds())
print("\tCreate object from dumped file: %f seconds" % (t4 - t3).total_seconds())
print("\tCompute stats using object: %f seconds" % (t3 - t2).total_seconds())

Warning 1: DeprecationWarning: 'Memory' driver is deprecated since GDAL 3.11. Use 'MEM' onwards. Further messages of this type will be suppressed.

0...10...20...30...40...50...60...70...80...90

Warning 1: Layer dst has no spatial index, DROP SPATIAL INDEX failed.

Warning 1: Layer src has no spatial index, DROP SPATIAL INDEX failed.
Warning 1: Layer trg has no spatial index, DROP SPATIAL INDEX failed.

Warning 1: Layer dst has no spatial index, DROP SPATIAL INDEX failed.

Approach #1 computation time:
	Create object from scratch: 3.372409 seconds
	Create object from dumped file: 1.312206 seconds
	Compute stats using object: 0.000031 seconds

src1 = wrl.io.VectorSource(radar_utmc, trg_crs=proj_utm, name="src")
trg1 = wrl.io.VectorSource(shpfile, trg_crs=proj_utm, name="trg", src_crs=proj_gk2)

# Just a test for plotting results with zero buffer
zd = wrl.zonalstats.ZonalDataPoint(src1, trg=trg1, buf=0)
# Create instance of type ZonalStatsPoint from zonal data object
obj2 = wrl.zonalstats.ZonalStatsPoint(zd)
obj2.zdata.trg.set_attribute("mean", avg1)
obj2.zdata.trg.set_attribute("var", var1)

isecs2 = obj2.zdata.isecs

Warning 1: Layer src has no spatial index, DROP SPATIAL INDEX failed.

Warning 1: Layer trg has no spatial index, DROP SPATIAL INDEX failed.
Warning 1: Layer dst has no spatial index, DROP SPATIAL INDEX failed.

...100 - done.
0...10...20...30...40...50...60...70...80...90

# Illustrate results for an example catchment i
i = 0  # try e.g. 5, 2
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, aspect="equal")

# Target polygon patches
trg_patch = obj2.zdata.trg.get_data_by_idx([i], mode="geo")
trg_patch.plot(ax=ax, facecolor="None", edgecolor="black", linewidth=2)
trg_patch = obj1.zdata.trg.get_data_by_idx([i], mode="geo")
trg_patch.plot(ax=ax, facecolor="None", edgecolor="grey", linewidth=2)

# pips
sources = obj1.zdata.src.geo
sources.plot(ax=ax, label="all points", c="grey", markersize=200)
isecs1 = obj2.zdata.dst.get_data_by_att(attr="trg_index", value=[i], mode="geo")
isecs1.plot(ax=ax, label="buffer=0 m", c="green", markersize=200)
isecs2 = obj1.zdata.dst.get_data_by_att(attr="trg_index", value=[i], mode="geo")
isecs2.plot(ax=ax, label="buffer=500 m", c="red", markersize=50)

cat = trg.get_data_by_idx([i])[0]
bbox = wrl.zonalstats.get_bbox(cat[..., 0], cat[..., 1])
plt.xlim(bbox["left"] - 2000, bbox["right"] + 2000)
plt.ylim(bbox["bottom"] - 2000, bbox["top"] + 2000)
plt.legend()
plt.title("Catchment #%d: Points considered for stats" % i)

Text(0.5, 1.0, 'Catchment #0: Points considered for stats')
../../_images/d2387849b2f54df8b0de5ff3d0f4c8733b09f2b92b9652513ac8a0f536b89571.png 
# Plot average rainfall and original data
testplot(
    ds_clip.data, obj2, col="mean", title="Catchment rainfall mean (ZonalStatsPoint)"
)
../../_images/2869200893e4d4de084b1eda5c368219fdaf28e76790cba1575bf1704834dd42.png 
testplot(
    ds_clip.data,
    obj2,
    col="var",
    levels=np.arange(0, 20, 1.0),
    title="Catchment rainfall variance (ZonalStatsPoint)",
)
../../_images/7d5767df64489a4b8e12c3daa8f103d340e8bd55c2c46b888680456da97d4bda.png

Zonal Stats Polar Grid - Polygons#

radar_utm = wrl.georef.spherical_to_polyvert(
    ds.range.values,
    ds.azimuth.values,
    0.5,
    (ds.longitude.values, ds.latitude.values, ds.altitude.values),
    crs=proj_utm,
)
radar_utm.shape = (360, 1000, 5, 3)
ds = ds.assign_coords(
    {
        "xp": (["azimuth", "range", "verts"], radar_utm[..., 0]),
        "yp": (["azimuth", "range", "verts"], radar_utm[..., 1]),
        "zp": (["azimuth", "range", "verts"], radar_utm[..., 2]),
    }
)
display(ds)
trg = wrl.io.VectorSource(shpfile, trg_crs=proj_utm, name="trg", src_crs=proj_gk2)
bbox = trg.extent

# create catchment bounding box
buffer = 5000.0
bbox = dict(
    left=bbox[0] - buffer,
    right=bbox[1] + buffer,
    bottom=bbox[2] - buffer,
    top=bbox[3] + buffer,
)
ds_clip = ds.where(
    (
        ((ds.yc > bbox["bottom"]) & (ds.yc < bbox["top"]))
        & ((ds.xc > bbox["left"]) & (ds.xc < bbox["right"]))
    ),
    drop=True,
)
display(ds_clip)

<xarray.Dataset> Size: 72MB
Dimensions:     (azimuth: 360, range: 1000, verts: 5)
Coordinates: (12/20)
  * azimuth     (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
  * range       (range) float64 8kB 100.0 200.0 300.0 ... 9.99e+04 1e+05
    x           (azimuth, range) float64 3MB 3.639e+05 3.639e+05 ... 3.657e+05
    y           (azimuth, range) float64 3MB 5.622e+06 5.622e+06 ... 5.722e+06
    z           (azimuth, range) float64 3MB 100.4 101.2 ... 1.559e+03 1.561e+03
    gr          (azimuth, range) float64 3MB 99.98 200.0 ... 9.986e+04 9.996e+04
    ...          ...
    latitude    float64 8B 50.73
    longitude   float64 8B 7.072
    altitude    float64 8B 99.5
    sweep_mode  <U20 80B 'azimuth_surveillance'
    elevation   float64 8B 0.5
    crs_wkt     int64 8B 0
Dimensions without coordinates: verts
Data variables:
    data        (azimuth, range) float64 3MB ...
<xarray.Dataset> Size: 12MB
Dimensions:     (azimuth: 86, range: 695, verts: 5)
Coordinates: (12/20)
  * azimuth     (azimuth) float64 688B 0.5 1.5 2.5 3.5 ... 82.5 83.5 84.5 85.5
  * range       (range) float64 6kB 2.8e+03 2.9e+03 3e+03 ... 7.21e+04 7.22e+04
    x           (azimuth, range) float64 478kB 3.64e+05 3.64e+05 ... 4.36e+05
    y           (azimuth, range) float64 478kB 5.624e+06 5.625e+06 ... 5.625e+06
    z           (azimuth, range) float64 478kB 124.4 125.3 ... 1.037e+03
    gr          (azimuth, range) float64 478kB 2.799e+03 2.899e+03 ... 7.217e+04
    ...          ...
    latitude    float64 8B 50.73
    longitude   float64 8B 7.072
    altitude    float64 8B 99.5
    sweep_mode  <U20 80B 'azimuth_surveillance'
    elevation   float64 8B 0.5
    crs_wkt     int64 8B 0
Dimensions without coordinates: verts
Data variables:
    data        (azimuth, range) float64 478kB nan nan nan nan ... nan nan nan
radar_utm = np.stack([ds_clip.xp, ds_clip.yp], axis=-1).reshape(-1, 5, 2)
print(radar_utm.shape)
src = wrl.io.VectorSource(radar_utm, trg_crs=proj_utm, name="src")
trg = wrl.io.VectorSource(shpfile, trg_crs=proj_utm, name="trg", src_crs=proj_gk2)

(59770, 5, 2)

###########################################################################
# Approach #2: Compute weighted mean based on fraction of source polygons
# in target polygons
#
# - This is more accurate (no assumptions), but probably slower...
###########################################################################

t1 = dt.datetime.now()

# Create instance of type ZonalDataPoly from source grid and
# catchment array
zd = wrl.zonalstats.ZonalDataPoly(src, trg=trg, crs=proj_utm)
# dump to file
zd.dump_vector("test_zonal_poly")
# Create instance of type ZonalStatsPoint from zonal data object
obj3 = wrl.zonalstats.ZonalStatsPoly(zd)

obj3.zdata.dump_vector("test_zonal_poly")
t2 = dt.datetime.now()


t3 = dt.datetime.now()

# Create instance of type ZonalStatsPoly from zonal data file
obj4 = wrl.zonalstats.ZonalStatsPoly("test_zonal_poly")

avg3 = obj4.mean(ds_clip.data.values.ravel())
var3 = obj4.var(ds_clip.data.values.ravel())


t4 = dt.datetime.now()

print("Approach #2 computation time:")
print("\tCreate object from scratch: %f seconds" % (t2 - t1).total_seconds())
print("\tCreate object from dumped file: %f seconds" % (t4 - t3).total_seconds())
print("\tCompute stats using object: %f seconds" % (t3 - t2).total_seconds())

obj4.zdata.trg.dump_raster(
    "test_zonal_hdr.nc", driver="netCDF", attr="mean", pixel_size=100.0
)

obj4.zdata.trg.dump_vector("test_zonal_shp")
obj4.zdata.trg.dump_vector("test_zonal_json.geojson", driver="GeoJSON")

...100 - done.
0...10...20...30...40..

.50...60...70...80.

..90

Approach #2 computation time:
	Create object from scratch: 2.175803 seconds
	Create object from dumped file: 1.222498 seconds
	Compute stats using object: 0.000032 seconds
...100 - done.
0

# Plot average rainfall and original data
testplot(
    ds_clip.data,
    obj4,
    col="mean",
    title="Catchment rainfall mean (PolarZonalStatsPoly)",
)
../../_images/10a5d568c65cbb6fd991629faa98d2994a4a139f40da3f453537c175c2cddf13.png 
testplot(
    ds_clip.data,
    obj4,
    col="var",
    levels=np.arange(0, 20, 1.0),
    title="Catchment rainfall variance (PolarZonalStatsPoly)",
)
../../_images/4807d896389de5ae90a14e30b76d094061b7fa26c0505a7a7436fecefddc0009.png 
# Illustrate results for an example catchment i
i = 0  # try e.g. 5, 2
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, aspect="equal")

# Grid cell patches
src_index = obj3.zdata.get_source_index(i)
trg_patch = obj3.zdata.src.get_data_by_idx(src_index, mode="geo")
trg_patch.plot(ax=ax, facecolor="None", edgecolor="black")

# Target polygon patches
trg_patch = obj3.zdata.trg.get_data_by_idx([i], mode="geo")
trg_patch.plot(ax=ax, facecolor="None", edgecolor="red", linewidth=2)

# intersections
isecs1 = obj3.zdata.dst.get_data_by_att(attr="trg_index", value=[i], mode="geo")
isecs1.plot(column="src_index", ax=ax, cmap=plt.cm.plasma, alpha=0.5)

cat = trg.get_data_by_idx([i])[0]
bbox = wrl.zonalstats.get_bbox(cat[..., 0], cat[..., 1])
plt.xlim(bbox["left"] - 2000, bbox["right"] + 2000)
plt.ylim(bbox["bottom"] - 2000, bbox["top"] + 2000)
plt.legend()
plt.title("Catchment #%d: Polygons considered for stats" % i)

Text(0.5, 1.0, 'Catchment #0: Polygons considered for stats')
../../_images/e871f3baed0521ccb8c91d3002c0f29dfa44ead6d92083f2ed37c9734ab59e39.png 
# Compare estimates
maxlim = np.max(np.concatenate((avg1, avg3)))
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, aspect="equal")
plt.scatter(avg1, avg3, edgecolor="None", alpha=0.5)
plt.xlabel("Average of points in or close to polygon (mm)")
plt.ylabel("Area-weighted average (mm)")
plt.xlim(0, maxlim)
plt.ylim(0, maxlim)
plt.plot([-1, maxlim + 1], [-1, maxlim + 1], color="black")
plt.show()
../../_images/607401e3612d0b591ea72b311f01b212a7f0748b77fe3646a02b4bf9259185dc.png