Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Small update to wradlib notebook #105

Merged
merged 3 commits into from
Sep 5, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/build-book-pullrequest.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ permissions:
contents: read

env:
DOCKER_TAG: pr_100
DOCKER_TAG: pr_105

jobs:
build-book:
Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/build-book.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ on:
type: string # had a lot of trouble with boolean types, see https://github.com/actions/runner/issues/1483

env:
DOCKER_TAG: pr_100
DOCKER_TAG: pr_105

jobs:
build-container:
Expand Down
2 changes: 1 addition & 1 deletion binder/Dockerfile
Original file line number Diff line number Diff line change
@@ -1 +1 @@
FROM ghcr.io/openradar/erad2024:pr_100
FROM ghcr.io/openradar/erad2024:pr_105
175 changes: 44 additions & 131 deletions notebooks/wradlib/wradlib_clutter_beamblockage.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# clutter and beamblockage"
"# wradlib - clutter and beamblockage"
]
},
{
Expand Down Expand Up @@ -223,14 +223,17 @@
"metadata": {},
"outputs": [],
"source": [
"radar_parameters[\"radar_beam_width_h\"]"
"bw = radar_parameters[\"radar_beam_width_h\"]\n",
"bw"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Prepare DEM for Polar Processing"
"## Prepare DEM for Polar Processing\n",
"\n",
"Here the power of [xr.apply_ufunc](https://docs.xarray.dev/en/stable/generated/xarray.apply_ufunc.html) is shown, a wrapper to xarray-ify numpy functions."
]
},
{
Expand All @@ -239,15 +242,31 @@
"metadata": {},
"outputs": [],
"source": [
"bw = radar_parameters[\"radar_beam_width_h\"]\n",
"y, x = xr.broadcast(dem.y, dem.x)\n",
"rastercoords = np.dstack([x, y])\n",
"def interpolate_dem(obj, dem, **kwargs):\n",
" dim0 = obj.wrl.util.dim0()\n",
"\n",
"polcoords = np.dstack([swp.x.values, swp.y.values])\n",
"# Map rastervalues to polar grid points\n",
"polarvalues = wrl.ipol.cart_to_irregular_spline(\n",
" rastercoords, dem.DEM, polcoords, order=3, prefilter=False\n",
")"
" def wrapper(sx, sy, dx, dy, dem, *args, **kwargs):\n",
" y, x = np.meshgrid(dy, dx)\n",
" rastercoords = np.dstack([x, y])\n",
" polcoords = np.dstack([sx, sy])\n",
" return wrl.ipol.cart_to_irregular_spline(rastercoords, dem, polcoords, **kwargs)\n",
"\n",
" out = xr.apply_ufunc(\n",
" wrapper,\n",
" obj.x,\n",
" obj.y,\n",
" dem.x,\n",
" dem.y,\n",
" dem,\n",
" input_core_dims=[[dim0, \"range\"], [dim0, \"range\"], [\"x\"], [\"y\"], [\"y\", \"x\"]],\n",
" output_core_dims=[[dim0, \"range\"]],\n",
" dask=\"parallelized\",\n",
" vectorize=True,\n",
" kwargs=kwargs,\n",
" dask_gufunc_kwargs=dict(allow_rechunk=True),\n",
" )\n",
" out.name = \"DEM\"\n",
" return obj.assign(DEM=out)"
]
},
{
Expand All @@ -256,25 +275,23 @@
"metadata": {},
"outputs": [],
"source": [
"swp = swp.assign(\n",
" DEM=([\"azimuth\", \"range\"], polarvalues),\n",
")"
"swp = interpolate_dem(swp, dem.DEM, order=3, prefilter=False)"
]
},
{
"cell_type": "markdown",
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Plot scan strategy"
"swp.DEM.wrl.vis.plot(cmap=\"terrain\", vmin=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"swp"
"## Plot scan strategy"
]
},
{
Expand Down Expand Up @@ -325,104 +342,6 @@
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Arbitrary cuts through volume"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sweeps = []\n",
"reindex = dict(angle_res=1, direction=1, start_angle=0, stop_angle=360)\n",
"for sn in dtree.groups:\n",
" if \"sweep\" in sn:\n",
" sweeps.append(\n",
" dtree[sn]\n",
" .ds.pipe(xd.util.reindex_angle, **reindex)\n",
" .set_coords(\"sweep_mode\")\n",
" .wrl.georef.georeference()\n",
" .set_coords(\"sweep_fixed_angle\")\n",
" )\n",
"vol = xr.concat(sweeps, dim=\"sweep_fixed_angle\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Pass meta variables to coords to avoid some issues\n",
"vol = vol.set_coords((\"sweep_mode\", \"sweep_number\", \"prt_mode\", \"follow_mode\"))\n",
"# Reduce coordinates so the georeferencing works\n",
"vol[\"elevation\"] = vol[\"elevation\"].median(\"azimuth\")\n",
"vol[\"time\"] = vol[\"time\"].min(\"azimuth\")\n",
"vol[\"sweep_mode\"] = vol[\"sweep_mode\"].min()\n",
"vol"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test extract single azimuth\n",
"azimuth = 180\n",
"rec_rhi = wrl.util.cross_section_ppi(vol, azimuth, method=\"nearest\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display(rec_rhi)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rec_rhi.DBZ.plot.contourf(x=\"range\", y=\"z\", ylim=(0, 12000), levels=21)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rec_cut = wrl.util.cross_section_ppi(\n",
" vol,\n",
" [(8.5, 45.8), (9.75, 45.3)],\n",
" crs=wrl.georef.get_earth_projection(),\n",
" method=\"nearest\",\n",
" npl=501,\n",
")\n",
"display(rec_cut)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rec_cut.DBZ.plot.contourf(\n",
" x=\"xyi\", y=\"z\", levels=21, ylim=(0, 12000), cmap=\"HomeyerRainbow\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -480,15 +399,6 @@
"ax4.set_title(\"Reflectivity clutter removed\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"swp"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -565,7 +475,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## BeamBlockage Calculation"
"## BeamBlockage Calculation\n",
"\n",
"Can you xarray-ify the following, too?"
]
},
{
Expand All @@ -575,7 +487,7 @@
"outputs": [],
"source": [
"beamradius = wrl.util.half_power_radius(swp.range, bw)\n",
"PBB = wrl.qual.beam_block_frac(polarvalues, swp.z.values, beamradius)\n",
"PBB = wrl.qual.beam_block_frac(swp.DEM.values, swp.z.values, beamradius)\n",
"PBB = np.ma.masked_invalid(PBB)\n",
"CBB = wrl.qual.cum_beam_block_frac(PBB)"
]
Expand Down Expand Up @@ -754,7 +666,7 @@
"paax.plot(thetap, ecp + alt[angle, :] - beamradius, \":b\", label=\"-3 dB Beam width\")\n",
"\n",
"# orography\n",
"paax.fill_between(thetap, ecp, ecp + polarvalues[angle, :], color=\"0.75\")\n",
"paax.fill_between(thetap, ecp, ecp + swp.DEM[angle, :], color=\"0.75\")\n",
"\n",
"# shape axes\n",
"cgax.set_xlim(0, np.max(ade))\n",
Expand Down Expand Up @@ -809,6 +721,7 @@
"\n",
"- [xarray](https://docs.xarray.dev)\n",
"- [dask](https://docs.dask.org/)\n",
"- [GDAL](https://gdal.org)\n",
"- [xradar backends](https://docs.openradarscience.org/projects/xradar/en/stable/importers.html)\n",
"- [CfRadial1](https://ncar.github.io/CfRadial/)"
]
Expand All @@ -825,7 +738,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.0"
"version": "3.12.4"
},
"nbdime-conflicts": {
"local_diff": [
Expand Down