|
#!/usr/bin/env python
|
|
|
|
import os
|
|
import pyart
|
|
import xarray as xr
|
|
import numpy as np
|
|
|
|
# Read the radar data from the given file.
|
|
|
|
infile = 'test.netcdf'
|
|
|
|
ds = xr.open_dataset(infile)
|
|
|
|
# Copy the input data and add _FillValue/missing_value attributes.
|
|
|
|
data = ds.DBZH.copy()
|
|
|
|
data.attrs = {'_FillValue':-999.,
|
|
'missing_value':-999.}
|
|
|
|
# Define the x,y dimension arrays for the curvilinear grid (longitude(y,x), latitude(y,x)).
|
|
|
|
y = np.arange(0, ds.azimuth.shape[0])
|
|
x = np.arange(0, ds.range.shape[0])
|
|
|
|
# Create the radar ppi object from the input data. Compute the longitude and latitude values with the Py-Art functions init_gate_altitude() and init_gate_longitude_latitude().
|
|
|
|
radar = pyart.testing.make_empty_ppi_radar(ds.range.shape[0], ds.azimuth.shape[0],1)
|
|
|
|
radar.latitude['data'] = np.array([float(ds.latitude.values)])
|
|
radar.longitude['data'] = np.array([float(ds.longitude.values)])
|
|
radar.range['range'] = np.array(ds.range)
|
|
radar.azimuth['data'] = np.array(ds.azimuth.values)
|
|
radar.sweep_number['data'] = np.array([float(ds.sweep_number.values)])
|
|
radar.sweep_start_ray_index['data'] = np.array([float(ds.sweep_start_ray_index)])
|
|
radar.sweep_end_ray_index['data'] = np.array([float(ds.sweep_end_ray_index)])
|
|
radar.fixed_angle['data'] = np.array(ds.fixed_angle)
|
|
radar.altitude['data'] = np.array([float(ds.altitude)])
|
|
radar.elevation['data'] = np.array(ds.fixed_angle) * len(ds.azimuth)
|
|
|
|
radar.init_gate_altitude()
|
|
radar.init_gate_longitude_latitude()
|
|
|
|
# Assign the coordinate data arrays.
|
|
|
|
lat = xr.DataArray(radar.gate_latitude['data'].copy(), coords=[y,x], dims=['y', 'x'])
|
|
lat.attrs = {'units': radar.gate_latitude['units'],
|
|
'standard_name': 'latitude'}
|
|
|
|
lon = xr.DataArray(radar.gate_longitude['data'].copy(), coords=[y,x], dims=['y', 'x'])
|
|
lon.attrs = {'units': radar.gate_longitude['units'],
|
|
'standard_name': 'longitude'}
|
|
|
|
# Create new dataset for the radar data and its curvilinear coordinates.
|
|
|
|
ds_new = xr.Dataset()
|
|
ds_new['ReflectivityQC'] = (('y','x'), data.data, data.attrs)
|
|
ds_new.coords['latitude'] = (('y','x'), lat.data, lat.attrs)
|
|
ds_new.coords['longitude'] = (('y','x'), lon.data, lon.attrs)
|
|
ds_new.attrs['history'] = 'DKRZ Tutorial example: radar data conversion'
|
|
|
|
# Write the new dataset to a netCDF output file.
|
|
|
|
outfile = 'radar_data_lonlat.nc'
|
|
|
|
ds_new.to_netcdf(outfile)
|
