Zonal Statistics - Polar Grid

Contents

Zonal Statistics - Polar Grid#

[1]:

import wradlib as wrl
import matplotlib.pyplot as pl
import matplotlib as mpl
import warnings

warnings.filterwarnings("ignore")
try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    pl.ion()
import numpy as np
import xarray as xr

/home/runner/micromamba-root/envs/wradlib-tests/lib/python3.11/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Zonal Stats Polar Grid - Points#

[11]:

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

[12]:

###########################################################################
# 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, srs=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())

Approach #1 computation time:
        Create object from scratch: 3.859122 seconds
        Create object from dumped file: 1.302764 seconds
        Compute stats using object: 0.000041 seconds

[13]:

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

# Just a test for plotting results with zero buffer
zd = wrl.zonalstats.ZonalDataPoint(src1, 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

[14]:

# Illustrate results for an example catchment i
i = 0  # try e.g. 5, 2
fig = pl.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])
pl.xlim(bbox["left"] - 2000, bbox["right"] + 2000)
pl.ylim(bbox["bottom"] - 2000, bbox["top"] + 2000)
pl.legend()
pl.title("Catchment #%d: Points considered for stats" % i)

[14]:

Text(0.5, 1.0, 'Catchment #0: Points considered for stats')

../../_images/notebooks_zonalstats_wradlib_zonalstats_polar_grid_17_1.png

[15]:

# Plot average rainfall and original data
testplot(
    ds_clip.data, obj2, col="mean", title="Catchment rainfall mean (ZonalStatsPoint)"
)

../../_images/notebooks_zonalstats_wradlib_zonalstats_polar_grid_18_0.png

[16]:

testplot(
    ds_clip.data,
    obj2,
    col="var",
    levels=np.arange(0, 20, 1.0),
    title="Catchment rainfall variance (ZonalStatsPoint)",
)

../../_images/notebooks_zonalstats_wradlib_zonalstats_polar_grid_19_0.png