I have LiDAR data that I’d like to convert to radar format (because I want to use pyart package) and then grid it. However, after gridding, I noticed that the ds.rsw values are NaN. Could you please advise on how to retrieve valid values for ds.rsw?
I have also attached the google-drive link to the lidar data.
lidar data
Code:
import pyart
import numpy as np
from datetime import datetime
from netCDF4 import Dataset
import warnings
warnings.filterwarnings('ignore')
# Load data lidar data
file_path = "D:/lidar_ppi.nc"
data = Dataset(file_path)
sweep_group = data.groups["Sweep_860123"]
time = sweep_group.variables["time"][:]
latitude = data.variables["latitude"][:]
longitude = data.variables["longitude"][:]
altitude = data.variables["altitude"][:]
azimuth = sweep_group.variables["azimuth"][:]
elevation = sweep_group.variables["elevation"][:]
range_ = sweep_group.variables["range"][:]
radial_wind_speed = sweep_group.variables["radial_wind_speed"][:]
rsw = np.array(radial_wind_speed)
# Create radar object using pyart
radar = pyart.testing.make_empty_ppi_radar(rsw.shape[1], len(azimuth), 1)
radar.latitude['data'] = np.array([latitude])
radar.longitude['data'] = np.array([longitude])
radar.altitude['data'] = np.array([altitude])
#radar.time['data'] = np.array([(t - time_converted[0]).total_seconds() for t in time_converted])
radar.time = {
'standard_name': 'time',
'long_name': 'time in seconds since volume start',
'calendar': 'gregorian',
'units': 'seconds since 2023-09-23T04:35:05Z',
'comment': 'times are relative to the volume start_time',
'data': np.array([(t - time_converted[0]).total_seconds() for t in time_converted]),
'FillValue':1e+20
}
radar.azimuth['data'] = azimuth
radar.elevation['data'] = elevation
radar.range['data'] = range
radial_wind_speed_dict = {
'long_name': 'radial_wind_speed',
'standard_name': 'radial_wind_speed_of_scatterers_away_from_instrument',
'units': 'm/s',
'sampling_ratio': 1.0,
'_FillValue': -9999 ,
'grid_mapping': 'grid_mapping',
'coordinates': 'time range',
'data': np.ma.masked_invalid(rsw) # Mask invalid data
}
radar.fields = { 'rws': radial_wind_speed_dict}
#success plot this:
import matplotlib.pyplot as plt
from pyart.graph import RadarDisplay
display = RadarDisplay(radar)
fig, ax = plt.subplots(figsize=(10, 8))
display.plot_ppi("rws", sweep=0, ax=ax, cmap="coolwarm")
plt.show()
## Now grid
grid_limits = ((10., 4000.), (-4500., 4500.), (-4500., 4500.))
grid_shape = (20, 50, 50)
grid = pyart.map.grid_from_radars([radar], grid_limits=grid_limits, grid_shape=grid_shape)
ds = grid_dv.to_xarray()
print(ds)
<xarray.Dataset> Size: 641kB
Dimensions: (time: 1, z: 20, y: 50, x: 50, nradar: 1)
Coordinates: (12/16)
* time (time) object 8B 2023-09-23 04:35:05
* z (z) float64 160B 10.0 220.0 ... 3.79e+03 4e+03
lat (y, x) float64 20kB 45.02 45.02 ... 45.1 45.1
lon (y, x) float64 20kB 7.603 7.605 ... 7.715 7.717
* y (y) float64 400B -4.5e+03 -4.316e+03 ... 4.5e+03
* x (x) float64 400B -4.5e+03 -4.316e+03 ... 4.5e+03
...
origin_altitude (time) float64 8B nan
radar_altitude (nradar) float64 8B nan
radar_latitude (nradar) float64 8B 45.06
radar_longitude (nradar) float64 8B 7.66
radar_time (nradar) int64 8B 0
radar_name (nradar) <U10 40B 'fake_radar'
Dimensions without coordinates: nradar
Data variables:
rws (time, z, y, x) float64 400kB nan nan ... nan
ROI (time, z, y, x) float32 200kB 500.0 ... 500.0
Attributes:
radar_name: fake_radar
nradar: 1
instrument_name: fake_radar
np.nanmax(ds.rws)
Out[18]: nan
np.nanmin(ds.rws)
Out[19]: nan