Decoding error reading CfRadial1 file

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

  1. time unit attribute is broken
  2. time_coverage_start is non-standard
  3. time_coverage_end is empty
  4. sweep_start_ray_index/sweep stop_ray_index with erroneous values

Thank you

Hi @adithiyrnair,

Would be great, if you could attach the data file too.

Update: Attaching data might not work as of now. Need to check with the backend.

Hi @adithiyrnair,

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))

You would have to tweak the plot as needed.

You also might have a look at @syedhamidali’s nice package syedhamidali/PyScanCf: 1.0.21 | Zenodo which is based on Py-ART and dedicated to IMD radars.

HTH,
Kai

Thank you Kai

Let me check it

start_time = “”.join(vol.time_coverage_start.astype(str).values)
error in this line showing

error : can only join an iterable

@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?

print(vol.time_coverage_start.astype(str).values)

this is the output

2017-11-30T00:01:00Z

This value is already correct. So this means in that case no correction is needed. Did you by chance run that particular code block twice?

Yeah by accident… Anyway thank you

new_vol.to_cfradial2(“test.nc”)

error: unsupported dtype for netCDF4 variable: datetime64[ns]

any help

Please send the output of wrl.show_versions(). This smells like dependency issue.

INSTALLED VERSIONS

commit: None
python: 3.8.13 (default, Oct 21 2022, 23:50:54)
[GCC 11.2.0]
python-bits: 64
OS: Linux
OS-release: 5.15.0-52-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_IN
LOCALE: (‘en_IN’, ‘ISO8859-1’)
libhdf5: 1.10.6
libnetcdf: 4.8.1

xarray: 0.20.1
pandas: 1.4.4
numpy: 1.23.3
scipy: 1.9.1
netCDF4: 1.5.7
pydap: None
h5netcdf: 1.0.2
h5py: 3.1.0
Nio: None
zarr: None
cftime: 1.5.1.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.3.5
dask: None
distributed: None
matplotlib: 3.1.3
cartopy: 0.18.0
seaborn: None
numbagg: None
fsspec: None
cupy: None
pint: None
sparse: None
setuptools: 65.5.0
pip: 22.2.2
conda: None
pytest: None
IPython: 8.4.0
sphinx: None

wradlib: 1.17.0

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.

yes it worked. Thanks again

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

GDAL is correct and this should work if GDAL is installed.

But are you sure you are working inside a conda environment? Your above output doesn’t mention conda.

But i am actually running in conda, I created a new enviroment only for wradlib in anaconda…and from that i am using

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.

(base) $ conda activate wradlib
(wradlib) $ conda install jupyterlab
(wradlib) $ conda list
(wradlib) $ jupyter lab

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 :laughing: