Hello everyone,
I currently working on RADOLAN binary data. I am trying to classify my data according to different Precipitation types available in the documentation Precipitation Types
Please guide me here if it is possible with my data available in the variable “clip_shape” be used to classify different precipitation types. My data variable is “RW”. You can see the file and data format below in the code. Thanks in advance.
import numpy as np
import pandas as pd
import scipy
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
import os
import sys
import tarfile
import fiona
import datetime
dict(xarray=xr.__version__, rioxarray=rioxarray.__version__, geopandas=gpd.__version__, cartopy=cartopy.__version__,
shapely=shapely.__version__, wradlib=wrl.__version__)
Out[60]:
{'xarray': '2023.2.0',
'rioxarray': '0.13.4',
'geopandas': '0.12.2',
'cartopy': '0.21.1',
'shapely': '2.0.1',
'wradlib': '1.19.1'}
# binary data files downloaded from https://opendata.dwd.de/climate_environment/CDC/grids_germany/hourly/radolan/historical/bin/
#Step 1:--------- load radolan raster data
fname = "Extracted Data/raa01-rw_10000-*-dwd---bin.gz"
rad = xr.open_mfdataset(fname, engine="radolan")
rad.RW.encoding["_FillValue"] = 65536 # fix encoding _FillValue
rad
#Step 2:------------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)
#Step 3:-----------Load German federal states shapefile
germany = gpd.read_file(filepath + "VG2500_LAN.shp")
germany.crs
germany.crs = "EPSG:31467"
#Step 4:-----------Extract Brandenburg
brandenburg = germany.loc[[11], "geometry"]
brandenburg.plot()
#Step 5:-----------Clip using rioxarray clip_box
bounds = brandenburg.to_crs(proj_radolan).bounds.iloc[0]
clip_box = rad.rio.clip_box(*bounds)
clip_box
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(projection=ccrs.PlateCarree())
clip_box.RW[0].plot(ax=ax, transform=proj_radolan)
ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)
#Step 6:-----------clip using rioxarray with shape
clip_shape = rad.rio.clip(brandenburg.geometry.apply(mapping), brandenburg.crs, drop=True)
clip_shape
clip_shape
Out[57]:
<xarray.Dataset>
Dimensions: (time: 8756, y: 254, x: 256)
Coordinates:
* time (time) datetime64[ns] 2006-01-01T00:45:00 ... 2006-12-31T23:...
* y (y) float64 -4.171e+06 -4.17e+06 ... -3.919e+06 -3.918e+06
* x (x) float64 8.904e+04 9.004e+04 9.104e+04 ... 3.43e+05 3.44e+05
spatial_ref int32 0
Data variables:
RW (time, y, x) float32 dask.array<chunksize=(1, 254, 256), meta=np.ndarray>
Attributes:
radarid: 10000
formatversion: 2
radolanversion: 01.01.00
radarlocations: ['bln', 'drs', 'eis', 'emd', 'ess', 'fbg', 'fld', 'fra',...