import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import wradlib as wrl
from matplotlib.colors import from_levels_and_colors
from osgeo import osr
warnings.filterwarnings("ignore")
plt.ioff()
BASE_DIR = os.path.dirname(os.path.abspath(__file__))
OUTPUT_DIR = os.path.join(BASE_DIR, "output")
if not os.path.exists(OUTPUT_DIR):
os.makedirs(OUTPUT_DIR)
PRODUCT_TYPE: str = "YW"
TARGET_EPSG: int = 25833
# Read and preprocess the RADOLAN data
file_path = f"{BASE_DIR}/examples/YW2017.002_20170102/raa01-yw2017.002_10000-1701021050-dwd---bin"
ds = wrl.io.open_radolan_mfdataset(file_path)
print(ds)
gridres = ds.x.diff("x")[0].values
print(gridres)
# This is the native RADOLAN projection
# (polar stereographic projection)
# create radolan projection osr object
if ds.attrs["formatversion"] >= 5:
proj_stereo = wrl.georef.create_osr("dwd-radolan-wgs84")
else:
proj_stereo = wrl.georef.create_osr("dwd-radolan-sphere")
# This is our target projection
proj_utm = osr.SpatialReference()
proj_utm.ImportFromEPSG(TARGET_EPSG)
coordinate_system_name = proj_utm.GetAttrValue("PROJCS")
# Get RADOLAN grid coordinates - center coordinates
x_rad, y_rad = np.meshgrid(ds.x, ds.y)
grid_xy_radolan = np.stack([x_rad, y_rad], axis=-1)
# Reproject the RADOLAN coordinates
xy = wrl.georef.reproject(grid_xy_radolan, src_crs=proj_stereo, trg_crs=proj_utm)
# Assign as coordinates
ds = ds.assign_coords(
{
"xc": (
["y", "x"],
xy[..., 0],
dict(long_name=f"{coordinate_system_name} Easting", units="m"),
),
"yc": (
["y", "x"],
xy[..., 1],
dict(long_name=f"{coordinate_system_name} Northing", units="m"),
),
}
)
trg_shp_file = f"{BASE_DIR}/examples/target_vectors/trg.shp"
trg = wrl.io.VectorSource(trg_shp_file, trg_crs=proj_utm, name="trg")
print(f"Found {len(trg)} sub-catchments in shapefile.")
# Clip subgrid from RADOLAN grid
bbox = trg.extent
buffer = 5000.0
bbox = dict(
left=bbox[0] - buffer,
right=bbox[1] + buffer,
bottom=bbox[2] - buffer,
top=bbox[3] + buffer,
)
print(bbox)
ds_clip = ds.where(
(
((ds.yc > bbox["bottom"]) & (ds.yc < bbox["top"]))
& ((ds.xc > bbox["left"]) & (ds.xc < bbox["right"]))
),
drop=True,
)
# Create vertices for each grid cell
# (MUST BE DONE IN NATIVE RADOLAN COORDINATES)
grid_x, grid_y = np.meshgrid(ds_clip.x, ds_clip.y)
grdverts = wrl.zonalstats.grid_centers_to_vertices(grid_x, grid_y, gridres, gridres)
src = wrl.io.VectorSource(grdverts, trg_crs=proj_utm, name="src", src_crs=proj_stereo)
# This object collects our source and target data
# and computes the intersections
zd = wrl.zonalstats.ZonalDataPoly(src, trg=trg, crs=proj_utm)
# zd = wrl.zonalstats.ZonalDataPoly(grdverts, shpfile, srs=proj_utm)
# This object can actually compute the statistics
obj = wrl.zonalstats.ZonalStatsPoly(zd)
# We just call this object with any set of radar data
avg = obj.mean(getattr(ds_clip, PRODUCT_TYPE).values.ravel())
max_level = avg.max(where=~np.isnan(avg), initial=0.0) + 0.5 # Set level limit to max average plus 0.5
# Create discrete colormap
levels = np.arange(0, max_level, 0.1)
colors = plt.cm.inferno(np.linspace(0, 1, len(levels)))
mycmap, mynorm = from_levels_and_colors(levels, colors, extend="max")
fig = plt.figure(figsize=(12, 6))
# Average rainfall sum
ax = fig.add_subplot(121, aspect="equal")
obj.zdata.trg.geo.plot(
column="mean",
ax=ax,
cmap=mycmap,
norm=mynorm,
edgecolor="white",
lw=0.5,
legend=True,
legend_kwds=dict(orientation="horizontal", pad=0.05),
)
plt.xlabel(f"{coordinate_system_name} Easting (m)")
plt.ylabel(f"{coordinate_system_name} Northing (m)")
plt.title("Catchment areal average")
bbox = obj.zdata.trg.extent
plt.xlim(bbox[0] - 5000, bbox[1] + 5000)
plt.ylim(bbox[2] - 5000, bbox[3] + 5000)
plt.grid()
# Original radar data
ax = fig.add_subplot(122, aspect="equal")
pm = getattr(ds_clip, PRODUCT_TYPE).plot(
x="xc",
y="yc",
cmap=mycmap,
norm=mynorm,
ax=ax,
cbar_kwargs=dict(orientation="horizontal", pad=0.05),
)
obj.zdata.trg.geo.plot(ax=ax, facecolor="None", edgecolor="white")
plt.title("RADOLAN rain depth")
plt.grid(color="white")
plt.tight_layout()
# export to vector GeoJSON
obj.zdata.trg.dump_vector(
f"{OUTPUT_DIR}/average_precipitation.geojson", driver="GeoJSON"
)
# export 'mean' to raster netCDF
obj.zdata.trg.dump_raster(
f"{OUTPUT_DIR}/average_precipitation.nc",
driver="netCDF",
attr="mean",
pixel_size=100.0,
)
# Convert geo data to pandas dataframe
df = pd.DataFrame(obj.zdata.trg.geo)
# Headers to show in CSV
header = [
"index",
"mean",
# "geometry"
]
# Export precipitation data to CSV
df.to_csv(f"{OUTPUT_DIR}/average_precipitation.csv", columns=header)
# Plot results in interactive map using geopandas
fmap = obj.zdata.trg.geo.explore(column="mean")
output_fp = f"{OUTPUT_DIR}/average_precipitation.html"
fmap.save(output_fp)
plt.show()