Created
December 16, 2022 09:01
-
-
Save kmuehlbauer/9d4f0e14f8025daafbb86bf7b4712bb9 to your computer and use it in GitHub Desktop.
Clip RADOLAN xarray dataset using rioxarray with shapefile
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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:
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]:
Area of Use:
Coordinate Operation:
Datum: Deutsches Hauptdreiecksnetz
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`