Sweep to Raster#
In this notebook we show the production of a maximum reflectivity composite from 3 neighboring radars for projected as well as geographic raster targets.
import tempfile
import warnings
import cmweather
import xarray as xr
import xradar as xd
import matplotlib.pyplot as plt
import wradlib as wrl
import wradlib_data
warnings.filterwarnings("ignore")
/home/docs/checkouts/readthedocs.org/user_builds/wradlib/conda/latest/lib/python3.13/site-packages/tqdm/auto.py:21: 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
Get radar data#
First, we import measurements from three belgian radars. This is done using xradar.io.backends.odim.open_odim_datatree.
filenames = ["bejab.pvol.hdf", "bewid.pvol.hdf", "behel.pvol.hdf"]
paths = [wradlib_data.DATASETS.fetch(f"hdf5/{f}") for f in filenames]
volumes = [xd.io.open_odim_datatree(p) for p in paths]
Downloading file 'hdf5/bejab.pvol.hdf' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/bejab.pvol.hdf' to '/home/docs/.cache/wradlib-data'.
Downloading file 'hdf5/bewid.pvol.hdf' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/bewid.pvol.hdf' to '/home/docs/.cache/wradlib-data'.
Downloading file 'hdf5/behel.pvol.hdf' from 'https://github.com/wradlib/wradlib-data/raw/main/data/hdf5/behel.pvol.hdf' to '/home/docs/.cache/wradlib-data'.
Projected coordinate raster#
Please see wradlib.georef.raster.create_raster_xarray for details.
europe_crs = 3035
bounds = [3614000, 2783000, 4200000, 3338000]
resolution = 1000
raster_projected = wrl.georef.create_raster_xarray(
crs=europe_crs,
bounds=bounds,
resolution=resolution,
)
display(raster_projected)
<xarray.Dataset> Size: 9kB
Dimensions: (x: 586, y: 555)
Coordinates:
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
spatial_ref int64 8B 0
Data variables:
*empty*Geographic coordinate raster#
Please see wradlib.georef.raster.create_raster_geographic for details.
bounds = [0, 48, 9, 53]
resolution = 1000
raster_geographic = wrl.georef.create_raster_geographic(
bounds=bounds, resolution=resolution, resolution_in_meters=True
)
display(raster_geographic)
<xarray.Dataset> Size: 10kB
Dimensions: (x: 648, y: 600)
Coordinates:
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
spatial_ref int64 8B 0
Data variables:
*empty*Transform with Sweep to Raster#
This transforms the lowest sweep into a raster for each radar. See wradlib.comp.sweep_to_raster for details.
raster = {}
raster["projected"] = raster_projected
raster["geographic"] = raster_geographic
rasters = {}
swp = 0
fig, axes = plt.subplots(3, 2, figsize=(13, 10))
for p, name in enumerate(["projected", "geographic"]):
metadata = xd.model.required_sweep_metadata_vars
rasters_radar = []
for volume in volumes:
sweep = volume[f"sweep_{swp}"].to_dataset(inherit="all_coords")
sweep = sweep[["DBZH"] + list(metadata)]
sweep = sweep.sel(range=slice(0, 200e3))
sweep = sweep.wrl.georef.georeference()
raster_radar = sweep.wrl.comp.sweep_to_raster(raster[name])
rasters_radar.append(raster_radar.drop_vars("crs_wkt"))
for i, raster_radar in enumerate(rasters_radar):
ax = axes[i, p]
raster_radar["DBZH"].plot(ax=ax, vmin=0, vmax=50, cmap="HomeyerRainbow")
ax.set_aspect("equal", "box")
ax.set_title(f"{name} - Radar {i}")
rasters[name] = rasters_radar
fig.tight_layout()
display(rasters["projected"])
display(rasters["geographic"])
[<xarray.Dataset> Size: 3MB
Dimensions: (y: 555, x: 586)
Coordinates:
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
latitude float64 8B 51.19
longitude float64 8B 3.064
altitude float64 8B 50.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan,
<xarray.Dataset> Size: 3MB
Dimensions: (y: 555, x: 586)
Coordinates:
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
latitude float64 8B 49.91
longitude float64 8B 5.506
altitude float64 8B 590.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan,
<xarray.Dataset> Size: 3MB
Dimensions: (y: 555, x: 586)
Coordinates:
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
latitude float64 8B 51.07
longitude float64 8B 5.406
altitude float64 8B 140.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan]
[<xarray.Dataset> Size: 3MB
Dimensions: (y: 600, x: 648)
Coordinates:
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
latitude float64 8B 51.19
longitude float64 8B 3.064
altitude float64 8B 50.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan,
<xarray.Dataset> Size: 3MB
Dimensions: (y: 600, x: 648)
Coordinates:
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
latitude float64 8B 49.91
longitude float64 8B 5.506
altitude float64 8B 590.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan,
<xarray.Dataset> Size: 3MB
Dimensions: (y: 600, x: 648)
Coordinates:
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
latitude float64 8B 51.07
longitude float64 8B 5.406
altitude float64 8B 140.0
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan]
Combine rasters#
rasters = {k: xr.concat(v, dim="sweep") for k, v in rasters.items()}
display(rasters["projected"])
display(rasters["geographic"])
<xarray.Dataset> Size: 8MB
Dimensions: (sweep: 3, y: 555, x: 586)
Coordinates:
latitude (sweep) float64 24B 51.19 49.91 51.07
longitude (sweep) float64 24B 3.064 5.506 5.406
altitude (sweep) float64 24B 50.0 590.0 140.0
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
spatial_ref int64 8B 0
Dimensions without coordinates: sweep
Data variables:
DBZH (sweep, y, x) float64 8MB nan nan nan nan ... nan nan nan nan<xarray.Dataset> Size: 9MB
Dimensions: (sweep: 3, y: 600, x: 648)
Coordinates:
latitude (sweep) float64 24B 51.19 49.91 51.07
longitude (sweep) float64 24B 3.064 5.506 5.406
altitude (sweep) float64 24B 50.0 590.0 140.0
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
spatial_ref int64 8B 0
Dimensions without coordinates: sweep
Data variables:
DBZH (sweep, y, x) float64 9MB nan nan nan nan ... nan nan nan nanCalculate maximum reflectivity#
rasters = {k: v.max(dim="sweep", keep_attrs=True) for k, v in rasters.items()}
display(rasters["projected"])
display(rasters["geographic"])
<xarray.Dataset> Size: 3MB
Dimensions: (y: 555, x: 586)
Coordinates:
* y (y) int64 4kB 3337500 3336500 3335500 ... 2784500 2783500
* x (x) int64 5kB 3614500 3615500 3616500 ... 4198500 4199500
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nan<xarray.Dataset> Size: 3MB
Dimensions: (y: 600, x: 648)
Coordinates:
* y (y) int64 5kB 190785 190755 190725 ... 172875 172845 172815
* x (x) int64 5kB 25 75 125 175 225 ... 32225 32275 32325 32375
spatial_ref int64 8B 0
Data variables:
DBZH (y, x) float64 3MB nan nan nan nan nan ... nan nan nan nan nanPlot results#
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
for i, (name, raster) in enumerate(rasters.items()):
raster["DBZH"].plot(ax=axes[i], vmin=0, vmax=50, cmap="HomeyerRainbow")
axes[i].set_title(name)
fig.tight_layout()
