Skip to content

Commit 81910c9

Browse files
kmuehlbauerkmuehlbauer
authored
FIX: dwd-grids (#65)
* FIX: dwd-grids * MNT: add `nc-time-axis` to test environment * MNT: use tempfile
1 parent 5aac6eb commit 81910c9

File tree

7 files changed

+118
-47
lines changed

7 files changed

+118
-47
lines changed

‎ci/requirements/notebooks.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ dependencies:
1616
- h5netcdf
1717
- lat_lon_parser
1818
- nbconvert
19+
- nc-time-axis
1920
- netCDF4
2021
- notebook
2122
- numpy

‎notebooks/fileio/wradlib_odim_backend.ipynb‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -476,7 +476,7 @@
476476
"name": "python",
477477
"nbconvert_exporter": "python",
478478
"pygments_lexer": "ipython3",
479-
"version": "3.9.6"
479+
"version": "3.10.0"
480480
}
481481
},
482482
"nbformat": 4,

‎notebooks/fileio/wradlib_xarray_radial_odim.ipynb‎

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@
4545
" get_ipython().magic(\"matplotlib inline\")\n",
4646
"except:\n",
4747
" pl.ion()\n",
48-
"from wradlib.io.xarray_depr import CfRadial, OdimH5"
48+
"from wradlib.io.xarray_depr import CfRadial, OdimH5\n",
49+
"import tempfile"
4950
]
5051
},
5152
{
@@ -313,8 +314,10 @@
313314
"metadata": {},
314315
"outputs": [],
315316
"source": [
316-
"cf1.to_odim('knmi_odim.h5')\n",
317-
"cf1.to_cfradial2('knmi_odim_as_cfradial.nc')"
317+
"h5_tempfile = tempfile.NamedTemporaryFile(suffix=\".h5\").name\n",
318+
"nc_tempfile = tempfile.NamedTemporaryFile(suffix=\".nc\").name\n",
319+
"cf1.to_odim(h5_tempfile)\n",
320+
"cf1.to_cfradial2(nc_tempfile)"
318321
]
319322
},
320323
{
@@ -330,8 +333,8 @@
330333
"metadata": {},
331334
"outputs": [],
332335
"source": [
333-
"cf1a = OdimH5('knmi_odim.h5', standard='cf', georef=True)\n",
334-
"cf1b = CfRadial('knmi_odim_as_cfradial.nc', georef=True)"
336+
"cf1a = OdimH5(h5_tempfile, standard='cf', georef=True)\n",
337+
"cf1b = CfRadial(nc_tempfile, georef=True)"
335338
]
336339
},
337340
{
@@ -528,8 +531,10 @@
528531
"metadata": {},
529532
"outputs": [],
530533
"source": [
531-
"cf2.to_cfradial2('timrex_cfradial2.nc')\n",
532-
"cf2.to_odim('timrex_cfradial_as_odim.h5')"
534+
"h5_tempfile = tempfile.NamedTemporaryFile(suffix=\".h5\").name\n",
535+
"nc_tempfile = tempfile.NamedTemporaryFile(suffix=\".nc\").name\n",
536+
"cf2.to_cfradial2(nc_tempfile)\n",
537+
"cf2.to_odim(h5_tempfile)"
533538
]
534539
},
535540
{
@@ -554,8 +559,8 @@
554559
"metadata": {},
555560
"outputs": [],
556561
"source": [
557-
"cf2a = CfRadial('timrex_cfradial2.nc')\n",
558-
"cf2b = OdimH5('timrex_cfradial_as_odim.h5', standard='cf')"
562+
"cf2a = CfRadial(nc_tempfile)\n",
563+
"cf2b = OdimH5(h5_tempfile, standard='cf')"
559564
]
560565
},
561566
{
@@ -614,7 +619,7 @@
614619
"name": "python",
615620
"nbconvert_exporter": "python",
616621
"pygments_lexer": "ipython3",
617-
"version": "3.9.6"
622+
"version": "3.10.0"
618623
}
619624
},
620625
"nbformat": 4,

‎notebooks/zonalstats/wradlib_zonal_stats_polar_grid.ipynb‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -541,7 +541,7 @@
541541
"name": "python",
542542
"nbconvert_exporter": "python",
543543
"pygments_lexer": "ipython3",
544-
"version": "3.9.9"
544+
"version": "3.10.0"
545545
},
546546
"toc": {
547547
"colors": {

‎notebooks/zonalstats/wradlib_zonal_stats_rectangular_grid.ipynb‎

Lines changed: 24 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -124,7 +124,8 @@
124124
"# Read and prepare the actual data (RADOLAN)\n",
125125
"f = wrl.util.get_wradlib_data_file(\n",
126126
" 'radolan/misc/raa01-sf_10000-1406100050-dwd---bin.gz')\n",
127-
"ds = wrl.io.open_radolan_dataset(f)"
127+
"ds = wrl.io.open_radolan_dataset(f)\n",
128+
"display(ds)"
128129
]
129130
},
130131
{
@@ -134,7 +135,10 @@
134135
"outputs": [],
135136
"source": [
136137
"# create radolan projection osr object\n",
137-
"proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
138+
"if ds.attrs[\"formatversion\"] >= 5:\n",
139+
" proj_stereo = wrl.georef.create_osr(\"dwd-radolan-wgs84\")\n",
140+
"else:\n",
141+
" proj_stereo = wrl.georef.create_osr(\"dwd-radolan-sphere\")\n",
138142
"\n",
139143
"# create UTM Zone 32 projection osr object\n",
140144
"proj_utm = osr.SpatialReference()\n",
@@ -151,7 +155,17 @@
151155
"metadata": {},
152156
"outputs": [],
153157
"source": [
154-
"print(proj_gk2)"
158+
"gridres = ds.x.diff(\"x\")[0].values\n",
159+
"gridres"
160+
]
161+
},
162+
{
163+
"cell_type": "code",
164+
"execution_count": null,
165+
"metadata": {},
166+
"outputs": [],
167+
"source": [
168+
"print(proj_stereo)"
155169
]
156170
},
157171
{
@@ -236,7 +250,7 @@
236250
"\n",
237251
"# Get RADOLAN center grid points for each grid cell\n",
238252
"# (MUST BE DONE IN NATIVE RADOLAN COORDINATES)\n",
239-
"grid_x, grid_y = np.meshgrid(ds_clip.x + 0.5, ds_clip.y + 0.5)\n",
253+
"grid_x, grid_y = np.meshgrid(ds_clip.x, ds_clip.y)\n",
240254
"grdpoints = np.dstack([grid_x, grid_y]).reshape(-1, 2)\n",
241255
"\n",
242256
"src = wrl.io.VectorSource(grdpoints, srs=proj_utm, name=\"src\", projection_source=proj_stereo)\n",
@@ -367,10 +381,10 @@
367381
"\n",
368382
"# Create vertices for each grid cell\n",
369383
"# (MUST BE DONE IN NATIVE RADOLAN COORDINATES)\n",
370-
"grid_x, grid_y = np.meshgrid(ds_clip.x + 0.5, ds_clip.y + 0.5)\n",
384+
"grid_x, grid_y = np.meshgrid(ds_clip.x, ds_clip.y)\n",
371385
"grdverts = wrl.zonalstats.grid_centers_to_vertices(grid_x,\n",
372386
" grid_y, \n",
373-
" 1., 1.)\n",
387+
" gridres, gridres)\n",
374388
"# And reproject to Cartesian reference system (here: UTM Zone 32)\n",
375389
"src = wrl.io.VectorSource(grdverts, srs=proj_utm, name=\"src\", projection_source=proj_stereo)\n",
376390
"trg = wrl.io.VectorSource(shpfile, srs=proj_utm, name=\"trg\", projection_source=proj_gk2)\n",
@@ -452,6 +466,9 @@
452466
"isecs1 = obj3.zdata.dst.get_data_by_att(attr=\"trg_index\", value=[i], mode=\"geo\")\n",
453467
"isecs1.plot(column=\"src_index\", ax=ax, cmap=pl.cm.plasma, alpha=0.5)\n",
454468
"\n",
469+
"# scatter center points\n",
470+
"ds_clip.plot.scatter(x=\"xc\", y=\"yc\")\n",
471+
"\n",
455472
"cat = trg.get_data_by_idx([i])[0]\n",
456473
"bbox = wrl.zonalstats.get_bbox(cat[..., 0], cat[..., 1])\n",
457474
"pl.xlim(bbox[\"left\"] - 2000, bbox[\"right\"] + 2000)\n",
@@ -492,7 +509,7 @@
492509
"name": "python",
493510
"nbconvert_exporter": "python",
494511
"pygments_lexer": "ipython3",
495-
"version": "3.9.9"
512+
"version": "3.10.0"
496513
},
497514
"toc": {
498515
"colors": {

‎notebooks/zonalstats/wradlib_zonalstats_example.ipynb‎

Lines changed: 57 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -134,14 +134,27 @@
134134
"ds = wrl.io.open_radolan_dataset(f)"
135135
]
136136
},
137+
{
138+
"cell_type": "code",
139+
"execution_count": null,
140+
"metadata": {},
141+
"outputs": [],
142+
"source": [
143+
"gridres = ds.x.diff(\"x\")[0].values\n",
144+
"gridres"
145+
]
146+
},
137147
{
138148
"cell_type": "code",
139149
"execution_count": null,
140150
"metadata": {},
141151
"outputs": [],
142152
"source": [
143153
"# create radolan projection osr object\n",
144-
"proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
154+
"if ds.attrs[\"formatversion\"] >= 5:\n",
155+
" proj_stereo = wrl.georef.create_osr(\"dwd-radolan-wgs84\")\n",
156+
"else:\n",
157+
" proj_stereo = wrl.georef.create_osr(\"dwd-radolan-sphere\")\n",
145158
"\n",
146159
"# create UTM Zone 32 projection osr object\n",
147160
"proj_utm = osr.SpatialReference()\n",
@@ -236,7 +249,7 @@
236249
"\n",
237250
"# Get RADOLAN center grid points for each grid cell\n",
238251
"# (MUST BE DONE IN NATIVE RADOLAN COORDINATES)\n",
239-
"grid_x, grid_y = np.meshgrid(ds_clip.x + 0.5, ds_clip.y + 0.5)\n",
252+
"grid_x, grid_y = np.meshgrid(ds_clip.x, ds_clip.y)\n",
240253
"grdpoints = np.dstack([grid_x, grid_y]).reshape(-1, 2)\n",
241254
"\n",
242255
"src = wrl.io.VectorSource(grdpoints, srs=proj_utm, name=\"src\", projection_source=proj_stereo)\n",
@@ -358,10 +371,10 @@
358371
"\n",
359372
"# Create vertices for each grid cell\n",
360373
"# (MUST BE DONE IN NATIVE RADOLAN COORDINATES)\n",
361-
"grid_x, grid_y = np.meshgrid(ds_clip.x + 0.5, ds_clip.y + 0.5)\n",
374+
"grid_x, grid_y = np.meshgrid(ds_clip.x, ds_clip.y)\n",
362375
"grdverts = wrl.zonalstats.grid_centers_to_vertices(grid_x,\n",
363376
" grid_y, \n",
364-
" 1., 1.)\n",
377+
" gridres, gridres)\n",
365378
"# And reproject to Cartesian reference system (here: UTM Zone 32)\n",
366379
"src = wrl.io.VectorSource(grdverts, srs=proj_utm, name=\"src\", projection_source=proj_stereo)\n",
367380
"trg = wrl.io.VectorSource(shpfile, srs=proj_utm, name=\"trg\", projection_source=proj_gk2)\n",
@@ -421,6 +434,15 @@
421434
" title=\"Catchment rainfall variance (ZonalStatsPoly)\")"
422435
]
423436
},
437+
{
438+
"cell_type": "code",
439+
"execution_count": null,
440+
"metadata": {},
441+
"outputs": [],
442+
"source": [
443+
"ds_clip"
444+
]
445+
},
424446
{
425447
"cell_type": "code",
426448
"execution_count": null,
@@ -445,12 +467,17 @@
445467
"isecs1 = obj3.zdata.dst.get_data_by_att(attr=\"trg_index\", value=[i], mode=\"geo\")\n",
446468
"isecs1.plot(column=\"src_index\", ax=ax, cmap=pl.cm.plasma, alpha=0.5)\n",
447469
"\n",
470+
"# scatter center points\n",
471+
"ds_clip.plot.scatter(x=\"xc\", y=\"yc\", s=10)\n",
472+
"\n",
448473
"cat = trg.get_data_by_idx([i])[0]\n",
449474
"bbox = wrl.zonalstats.get_bbox(cat[..., 0], cat[..., 1])\n",
450475
"pl.xlim(bbox[\"left\"] - 2000, bbox[\"right\"] + 2000)\n",
451476
"pl.ylim(bbox[\"bottom\"] - 2000, bbox[\"top\"] + 2000)\n",
452477
"pl.legend()\n",
453-
"pl.title(\"Catchment #%d: Polygons considered for stats\" % i)"
478+
"pl.title(\"Catchment #%d: Polygons considered for stats\" % i)\n",
479+
"#pl.gca().set_xlim(402000, 404000)\n",
480+
"#pl.gca().set_ylim(5642000, 5644000)"
454481
]
455482
},
456483
{
@@ -487,9 +514,9 @@
487514
"source": [
488515
"def create_center_coords(ds, proj=None):\n",
489516
" # create polar grid centroids in GK2\n",
490-
" center = wrl.georef.spherical_to_centroids(ds.range.values, \n",
517+
" center = wrl.georef.spherical_to_centroids(ds.data.r, \n",
491518
" ds.azimuth.values, \n",
492-
" 0.5, \n",
519+
" ds.elevation.values, \n",
493520
" (ds.longitude.values, ds.latitude.values, ds.altitude.values),\n",
494521
" proj=proj)\n",
495522
" ds = ds.assign_coords({\"xc\": ([\"azimuth\", \"range\"], center[..., 0]),\n",
@@ -512,11 +539,19 @@
512539
" \"longitude\": ds.data.Longitude, \n",
513540
" \"altitude\": 99.5,\n",
514541
" \"azimuth\": ds.data.az,\n",
515-
" \"range\": ds.data.r,\n",
542+
" # bin centers\n",
543+
" \"range\": ds.data.r - np.median(np.diff(ds.data.r)) / 2.,\n",
516544
" \"sweep_mode\": \"azimuth_surveillance\",\n",
517545
" \"elevation\": 0.5}\n",
518-
" )\n",
519-
"\n",
546+
" )"
547+
]
548+
},
549+
{
550+
"cell_type": "code",
551+
"execution_count": null,
552+
"metadata": {},
553+
"outputs": [],
554+
"source": [
520555
"ds = ds.pipe(wrl.georef.georeference_dataset, proj=proj_utm)\n",
521556
"ds = ds.pipe(create_center_coords, proj=proj_utm)\n",
522557
"display(ds)"
@@ -699,7 +734,7 @@
699734
"metadata": {},
700735
"outputs": [],
701736
"source": [
702-
"radar_utm = wrl.georef.spherical_to_polyvert(ds.range.values, \n",
737+
"radar_utm = wrl.georef.spherical_to_polyvert(ds.range.values + np.median(np.diff(ds.range.values)) / 2., \n",
703738
" ds.azimuth.values, \n",
704739
" 0.5, \n",
705740
" (ds.longitude.values, ds.latitude.values, ds.altitude.values),\n",
@@ -708,7 +743,7 @@
708743
"ds = ds.assign_coords({\"xp\": ([\"azimuth\", \"range\", \"verts\"], radar_utm[..., 0]),\n",
709744
" \"yp\": ([\"azimuth\", \"range\", \"verts\"], radar_utm[..., 1]),\n",
710745
" \"zp\": ([\"azimuth\", \"range\", \"verts\"], radar_utm[..., 2])})\n",
711-
"display(ds)\n",
746+
"\n",
712747
"trg = wrl.io.VectorSource(shpfile, srs=proj_utm, name=\"trg\", projection_source=proj_gk2)\n",
713748
"bbox = trg.extent\n",
714749
"\n",
@@ -779,11 +814,11 @@
779814
"print(\n",
780815
"\"\\tCompute stats using object: %f seconds\" % (t3 - t2).total_seconds())\n",
781816
"\n",
782-
"obj3.zdata.trg.dump_raster('test_zonal_hdr.nc', 'netCDF', 'mean',\n",
817+
"obj4.zdata.trg.dump_raster('test_zonal_hdr.nc', 'netCDF', 'mean',\n",
783818
" pixel_size=100.)\n",
784819
"\n",
785-
"obj3.zdata.trg.dump_vector('test_zonal_shp')\n",
786-
"obj3.zdata.trg.dump_vector('test_zonal_json.geojson', 'GeoJSON')\n",
820+
"obj4.zdata.trg.dump_vector('test_zonal_shp')\n",
821+
"obj4.zdata.trg.dump_vector('test_zonal_json.geojson', 'GeoJSON')\n",
787822
"\n",
788823
"# Target polygon patches\n",
789824
"trg_patches = [patches.Polygon(item, True) for item in obj3.zdata.trg.data]"
@@ -843,12 +878,17 @@
843878
"isecs1 = obj3.zdata.dst.get_data_by_att(attr=\"trg_index\", value=[i], mode=\"geo\")\n",
844879
"isecs1.plot(column=\"src_index\", ax=ax, cmap=pl.cm.plasma, alpha=0.5)\n",
845880
"\n",
881+
"# scatter center points\n",
882+
"ds_clip.plot.scatter(x=\"xc\", y=\"yc\", s=10)\n",
883+
"\n",
846884
"cat = trg.get_data_by_idx([i])[0]\n",
847885
"bbox = wrl.zonalstats.get_bbox(cat[..., 0], cat[..., 1])\n",
848886
"pl.xlim(bbox[\"left\"] - 2000, bbox[\"right\"] + 2000)\n",
849887
"pl.ylim(bbox[\"bottom\"] - 2000, bbox[\"top\"] + 2000)\n",
850888
"pl.legend()\n",
851-
"pl.title(\"Catchment #%d: Polygons considered for stats\" % i)"
889+
"pl.title(\"Catchment #%d: Polygons considered for stats\" % i)\n",
890+
"pl.gca().set_xlim(402000, 404000)\n",
891+
"pl.gca().set_ylim(5654000, 5656000)\n"
852892
]
853893
},
854894
{
@@ -883,7 +923,7 @@
883923
"name": "python",
884924
"nbconvert_exporter": "python",
885925
"pygments_lexer": "ipython3",
886-
"version": "3.9.9"
926+
"version": "3.10.0"
887927
},
888928
"toc": {
889929
"colors": {

0 commit comments

Comments
 (0)