I am trying to plot CAPPI and MAX Z using wradlib. The following is the error unable to decode time units ‘seconds since yyyy-mm-ddThh:mm:ssZ’ with “calendar ‘gregorian’”. Try opening your dataset with decode_times=False or installing cftime if it is not installed. After decode_times=False the the reading is happening through wradlib. The Coordinates for variables are time and range . I am attaching a sample of the file that want to plot ***A notebook for ploting CAPPI and MAX is much apreciated.
The issues are
time unit attribute is broken
time_coverage_start is non-standard
time_coverage_end is empty
sweep_start_ray_index/sweep stop_ray_index with erroneous values
I hope the following code snippets help to understand the issues with the file.
# Just open as simple xarray.Dataset, not decoding times
vol = xr.open_dataset("RCTLS_30NOV2017_000957_L2B_STD.nc", decode_times=False)
# 1./2. fix time unit and time_coverage_start
start_time = "".join(vol.time_coverage_start.astype(str).values)
start_time = dt.datetime.strptime(start_time, "%Y-%m-%dT%H:%M:%SZ")
vol["time_coverage_start"] = start_time
# time-units is wrong: "seconds since yyyy-mm-ddThh:mm:ssZ"
vol["time"].attrs["units"] = f"seconds since {start_time}"
# new correct value: "seconds since 2017-11-30 00:09:57"
# now we can decode time (!)
vol = xr.decode_cf(vol)
# 4. fix ray_index
# we have wrong values here, so that we cannot extract the individual sweeps
# sweep_start_ray_index: array([0, 0, 0], dtype=int32)
# sweep_end_ray_index: array([0, 0, 0], dtype=int32)
# to resolve this we identify these indices with the help of elevation
# get unique elevations (removing "0" in this case, as it is not an elevation)
ele = sorted(list(set(vol.elevation.values) ^ {0}))
vele = vol.elevation
sweep_start_ray_index = [vele.where(vele==el).argmin("time").values.item() for el in ele]
sweep_end_ray_index = [vol.dims["time"] - vele.where(vele==el)[::-1].argmin("time").values.item() for el in ele]
vol["sweep_start_ray_index"].values = sweep_start_ray_index
vol["sweep_end_ray_index"].values = sweep_end_ray_index
# assign data to wradlib RadarVolume()
# now we have everything in place and can extract individual sweeps
# this is done via private functions of wradlib, so take with care
sweeps = [wrl.io.xarray._assign_data_radial(vol, sweep=f"sweep_{i+1}")[0] for i in range(vol.dims["sweep"])]
new_vol = wrl.io.RadarVolume(engine="cfradial1")
new_vol.extend(sweeps)
new_vol.sort(key=lambda x: x.time.min().values)
# write to CfRadial2-like file (which can be read by new package https://github.com/openradar/xradar)
new_vol.to_cfradial2("test_cfradial2.nc")
After that you can interpolate your data into a 3D-grid and plot MAX CAPPI like:
# Iterate over sweeps and process
swp_list = []
for v in new_vol:
# georeference and stack dimensions
swp = v.pipe(wrl.georef.georeference_dataset).stack(points=["azimuth", "range"])
swp_list.append(swp)
# concat sweeps to volume
vol0 = xr.concat(swp_list, dim="points")
# Create XYZ Coordinate DataArray
xyz = xr.concat([vol0.x, vol0.y, vol0.z], dim="xyz").transpose()
# Create Target 3D Grid, choose your wanted resolution
xcnt = 100
ycnt = 100
zcnt = 200
trgx = np.linspace(xyz[:, 0].min(),xyz[:, 0].max(), xcnt)
trgy = np.linspace(xyz[:, 1].min(),xyz[:, 1].max(), ycnt)
trgz = np.linspace(0, 20000, zcnt)
yy, hh, xx = np.meshgrid(trgy,trgz,trgx)
trgxyz = np.stack([xx.flatten(), yy.flatten(), hh.flatten()]).T
# Create Gridder/Interpolator
trgshape=xx.shape
gridder = wrl.vpr.CAPPI(polcoords=xyz,
gridcoords=trgxyz,
gridshape=trgshape,
maxrange=200000,
minelev=0,
maxelev=90,
ipclass=wrl.ipol.Nearest)
# Interpolate Data into 3D Grid
crtd_ref = vol0.DBZ
vol_zh = np.ma.masked_invalid(gridder(crtd_ref.values).reshape(trgshape))
trgx = trgxyz[:, 0].reshape(trgshape)[0, 0, :]
trgy = trgxyz[:, 1].reshape(trgshape)[0, :, 0]
trgz = trgxyz[:, 2].reshape(trgshape)[:, 0, 0]
yy, hh, xx = np.meshgrid(trgy,trgz,trgx)
# Plot MAXCAPPI
wrl.vis.plot_max_plan_and_vert(trgx, trgy, trgz, vol_zh, unit="dBZH",
cmap="jet", levels=range(-10, 60))
@adithiyrnair Unfortunately I can’t reproduce this error with either of your two files.
The error message means, that vol.time_coverage_start.astype(str).values is not iterable (eg. a scalar and not a list/array). Can you print the content?
Please update “xarray” (I have xarray: 2022.10.0). Also consider to install nc_time_axis, which might bring in some convenience code support for xarray/netCDF4.
AttributeError: Module ‘osgeo.osr’ is not installed.
You tried to access function/module/attribute ‘SpatialReference’
from module ‘osgeo.osr’.
This module is optional right now in wradlib.
You need to separately install this dependency.
Please refer to Installation — wradlib
for further instructions.
I already installed GDAL. Do i need to seperately install, I am Using Conda… Where can i download that perticular package. In conda site the package is in the name of GDAL
Are you running the code via jupyter lab/notebook? The environment might be missing the jupyter kernel. For a quick verification just install jupyterlab inside your wradlib environment. Make sure to activate your environment and start lab from there.
hi @kmuehlbauer
The problem got solved after upgrading numpy. I found from stack overflow that installing Gdal will actually downgrade the numpoy version. so after installlation of Gdal we have to upgrade numpy to work OSgeo.osr properly.
I couldn’t replay to you yesterday because of chat limit shown by the site
