-
-
Save kmuehlbauer/9d4f0e14f8025daafbb86bf7b4712bb9 to your computer and use it in GitHub Desktop.
Dear @kmuehlbauer,
Everything is good until until this line of code:
bld = gpd.read_file("/automount/db01/python/data/ADM/germany/vg250_0101.gk3.shape.ebenen/vg250_ebenen/vg250_bld.shp") bld.crs
When I run:
brandenburg = bld.loc[[9], "geometry"] brandenburg
it throws me the following error:
Complete code:
import numpy as np
import wradlib as wrl
import xarray as xr
import rioxarray
import geopandas as gpd
import cartopy
import cartopy.crs as ccrs
import shapely
from shapely.geometry import mapping
import matplotlib.pyplot as plt
dict(xarray=xr.__version__, rioxarray=rioxarray.__version__, geopandas=gpd.__version__, cartopy=cartopy.__version__, shapely=shapely.__version__, wradlib=wrl.__version__)
Out[9]: {'xarray': '2022.10.0', 'rioxarray': '0.13.2', 'geopandas': '0.12.2', 'cartopy': '0.20.0', 'shapely': '1.8.0', 'wradlib': '1.18.0'}
load radolan raster data
fname = "D:/Sandeep/Thesis/Data/BIN/Extracted/raa01-rw_10000-*-dwd---bin.gz"
rad = xr.open_mfdataset(fname, engine="radolan")
fix encoding _FillValue
rad.RW.encoding["_FillValue"] = 65536
rad
`Out[10]:
<xarray.Dataset>
Dimensions: (time: 8758, y: 900, x: 900)
Coordinates:
- time (time) datetime64[ns] 2006-01-01T00:45:00 ... 2006-12-31T23:...
- y (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
- x (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
spatial_ref int32 0
Data variables:
RW (time, y, x) float32 dask.array<chunksize=(1, 900, 900), meta=np.ndarray>
Attributes:
radarid: 10000
formatversion: 2
radolanversion: 01.01.00
radarlocations: ['bln', 'drs', 'eis', 'emd', 'ess', 'fbg', 'fld', 'fra',...`
setup projection
proj_radolan = ccrs.Stereographic( true_scale_latitude=60.0, central_latitude=90.0, central_longitude=10.0 )
rad.rio.set_spatial_dims(x_dim="x", y_dim="y", inplace=True)
rad.rio.write_crs(proj_radolan, inplace=True)
Load German federal states shapefile
germany = gpd.read_file(r'D:\Sandeep\Thesis\vg2500_12-31.gk3.shape\vg2500\VG2500_LAN.shp')
germany.crs
`Out[11]:
Name: DHDN / 3-degree Gauss-Kruger zone 3
Axis Info [cartesian]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
Area of Use: - name: Germany - former West Germany onshore between 7°30'E and 10°30'E - states of Baden-Wurtemberg, Bayern, Bremen, Hamberg, Hessen, Niedersachsen, Nordrhein-Westfalen, Rhineland-Pfalz, Schleswig-Holstein.
- bounds: (7.5, 47.27, 10.51, 55.09)
Coordinate Operation: - name: 3-degree Gauss-Kruger zone 3
- method: Transverse Mercator
Datum: Deutsches Hauptdreiecksnetz - Ellipsoid: Bessel 1841
- Prime Meridian: Greenwich`
Extract Brandenburg
brandenburg = germany.loc[[11], "geometry"]
brandenburg
`brandenburg
Traceback (most recent call last):
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\IPython\core\formatters.py", line 702, in call
printer.pretty(obj)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\IPython\lib\pretty.py", line 394, in pretty
return _repr_pprint(obj, self, cycle)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\IPython\lib\pretty.py", line 700, in _repr_pprint
output = repr(obj)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\core\series.py", line 1594, in repr
return self.to_string(**repr_params)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\core\series.py", line 1687, in to_string
result = formatter.to_string()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 397, in to_string
fmt_values = self._get_formatted_values()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 381, in _get_formatted_values
return format_array(
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1328, in format_array
return fmt_obj.get_result()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1359, in get_result
fmt_values = self._format_strings()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1662, in _format_strings
fmt_values = format_array(
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1328, in format_array
return fmt_obj.get_result()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1359, in get_result
fmt_values = self._format_strings()
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1422, in _format_strings
fmt_values.append(f" {_format(v)}")
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\pandas\io\formats\format.py", line 1402, in _format
return str(formatter(x))
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\geopandas\array.py", line 1306, in
return lambda geom: shapely.wkt.dumps(geom, rounding_precision=precision)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\shapely\wkt.py", line 62, in dumps
return geos.WKTWriter(geos.lgeos, trim=trim, **kw).write(ob)
File "C:\Users\Admin\anaconda3\envs\wradlib\lib\site-packages\shapely\geos.py", line 401, in write
result = self._lgeos.GEOSWKTWriter_write(self._writer, geom._geom)
OSError: exception: access violation writing 0x0000000000000000`
This is wonderful. For example considering rad_2006 as the radolan dataset that was saved as rad_2006.nc in the local drive. Is it possible to resample the hourly data to quarterly like the below line:
rad_2006_q = rad_2006.resample(time='1Q').mean()
I am little skeptical to use mean() or sum() as the aggregate function. As this is the rainfall data I assume sum() is a better choice over mean(). Please educate me.