|
#!/usr/bin/env python
|
|
# coding: utf-8
|
|
|
|
import os
|
|
import pyart
|
|
import xarray as xr
|
|
import numpy as np
|
|
|
|
import matplotlib.pyplot as plt
|
|
import cartopy.crs as ccrs
|
|
import cartopy.feature as cfeature
|
|
|
|
infile = os.environ['HOME'] + '/Downloads/tmp.netcdf'
|
|
|
|
ds = xr.open_dataset(infile)
|
|
|
|
data1 = ds.ReflectivityQC
|
|
|
|
print('Data min: ', data1.min().values, ' max: ', data1.max().values)
|
|
|
|
#-- Create an empty ppi radar object.
|
|
|
|
radar1 = pyart.testing.make_empty_ppi_radar(ds.range.shape[0], ds.azimuth.shape[0],1)
|
|
|
|
#-- Broadcast the data components to the radar object.
|
|
|
|
radar1.latitude['data'] = np.array([float(ds.latitude.values)])
|
|
radar1.longitude['data'] = np.array([float(ds.longitude.values)])
|
|
|
|
radar1.range['data'] = np.array(ds.range)
|
|
radar1.azimuth['data'] = np.array(ds.azimuth.values)
|
|
|
|
radar1.sweep_number['data'] = np.array([float(ds.sweep_number.values)])
|
|
radar1.sweep_start_ray_index['data'] = np.array([float(ds.sweep_start_ray_index)])
|
|
radar1.sweep_end_ray_index['data'] = np.array([float(ds.sweep_end_ray_index)])
|
|
|
|
radar1.fixed_angle['data'] = np.array(ds.fixed_angle)
|
|
radar1.altitude['data'] = np.array([float(ds.altitude)])
|
|
|
|
#-- Compute the elevation by the fixed_angle and azimuth.
|
|
|
|
radar1.elevation['data'] = np.array(ds.fixed_angle) * len(ds.azimuth)
|
|
|
|
#-- Recalculate gate latitude, longitude and altitude.
|
|
|
|
radar1.init_gate_altitude()
|
|
radar1.init_gate_longitude_latitude()
|
|
|
|
print('Longitude: ',radar1.gate_longitude['data'].min(),radar1.gate_longitude['data'].max())
|
|
print('Latitude: ',radar1.gate_latitude['data'].min(), radar1.gate_latitude['data'].max())
|
|
|
|
#-- Create the plot of the ReflectivityQC data using a simple Matplotlib script
|
|
|
|
projection = ccrs.PlateCarree()
|
|
|
|
fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection=projection))
|
|
|
|
ax.add_feature(cfeature.LAND, facecolor='gainsboro')
|
|
ax.gridlines(draw_labels=True)
|
|
ax.coastlines(resolution='10m')
|
|
|
|
plot = plt.pcolormesh(radar1.gate_longitude['data'],
|
|
radar1.gate_latitude['data'],
|
|
data1,
|
|
shading='auto',
|
|
transform=ccrs.PlateCarree(),
|
|
zorder=9)
|
|
cbar = plt.colorbar(plot, ax=ax, shrink=0.4, pad=0.1)
|
|
|
|
fig.savefig('./plot_radar_range_azimuth_convert_to_lonlat_v2.png',
|
|
facecolor='white',
|
|
bbox_inches='tight',
|
|
dpi=100)
|
|
#plt.show()
|
|
|
|
#-- Write data to new netCDF file
|
|
|
|
yc = np.arange(0, ds.azimuth.shape[0])
|
|
xc = np.arange(0, ds.range.shape[0])
|
|
|
|
lat = xr.DataArray(radar1.gate_latitude['data'], coords=[yc,xc], dims=['y', 'x'])
|
|
lat.attrs = {'units': radar1.gate_latitude['units'],
|
|
'standard_name': 'latitude'}
|
|
|
|
lon = xr.DataArray(radar1.gate_longitude['data'], coords=[yc,xc], dims=['y', 'x'])
|
|
lon.attrs = {'units': radar1.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'), data1.data, data1.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'
|
|
|
|
outfile = 'radar_data_lonlat.nc'
|
|
|
|
ds_new.to_netcdf(outfile)
|
|
|
|
#------------------------------------------------------------------------------
|
|
#-- Remap output data with python-cdo
|
|
|
|
from cdo import Cdo
|
|
cdo = Cdo()
|
|
|
|
grid_text = '''gridtype = lonlat
|
|
xsize = 1200
|
|
ysize = 1200
|
|
xfirst = -83.
|
|
xinc = 0.005
|
|
yfirst = 34.
|
|
yinc = 0.005
|
|
'''
|
|
|
|
gridfile = './gridfile_radar_0.005.txt'
|
|
|
|
with open(gridfile, 'w') as gf:
|
|
gf.write(grid_text)
|
|
|
|
cdo.remapnn(gridfile,
|
|
input=' -setmissval,-999. '+outfile,
|
|
output='radar_data_lonlat_0.005.nc')
|
|
|
|
#-- Plot the remapped data
|
|
|
|
ds_test = xr.open_dataset('radar_data_lonlat_0.005.nc')
|
|
print(ds_test)
|
|
|
|
projection = ccrs.PlateCarree()
|
|
|
|
fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection=projection))
|
|
|
|
ax.add_feature(cfeature.LAND, facecolor='gainsboro')
|
|
ax.gridlines(draw_labels=True)
|
|
ax.coastlines(resolution='10m')
|
|
|
|
plot = plt.pcolormesh(ds_test.lon,
|
|
ds_test.lat,
|
|
ds_test.ReflectivityQC,
|
|
shading='auto',
|
|
transform=ccrs.PlateCarree(),
|
|
zorder=9)
|
|
cbar = plt.colorbar(plot, ax=ax, shrink=0.4, pad=0.1)
|
|
|
|
fig.savefig('./plot_radar_range_azimuth_convert_to_lonlat_remapped_0.005.png',
|
|
facecolor='white',
|
|
bbox_inches='tight',
|
|
dpi=100)
|
|
#plt.show()
|
|
|
|
|