In [143]:
import pyart
import dask.array as da
import dask.bag as db
from dask import delayed
from distributed import Client
from dask_jobqueue import PBSCluster
from glob import glob
from datetime import datetime
%pylab inline
Use calibration technique of Bringi and Chandrashekar (2001) that estimates bias from theoretically expected ZDR in light rain. Need good light rain case. Use 2017-10-04 1600001Z
In [41]:
file_path = '/lustre/or-hydra/cades-arm/proj-shared/sgpxsaprcmacsurI5.c1/20171004/sgpxsaprcmacsurI5.c1.20171004.160022.nc'
radar = pyart.io.read(file_path)
In [42]:
radar.fields.keys()
Out[42]:
In [4]:
def zdr_adjusted(zdr0, theta):
theta_rad = np.deg2rad(theta)
zdr_adj = zdr0/(zdr0**0.5*np.sin(theta_rad)**2+np.cos(theta_rad)**2)**2
return zdr_adj
In [140]:
def get_zdr_values(file):
try:
radar = pyart.io.read(file, delay_field_loading=True)
except OSError:
return np.nan
# Check for 40 dBZ echoes above 4 km
gate_z = radar.gate_z['data']
ref = radar.fields['reflectivity']['data']
snr = radar.fields['SNR']['data']
rhohv = radar.fields['cross_correlation_ratio']['data']
elevation = radar.elevation['data']
zdr = radar.fields['differential_reflectivity']['data']
gate_id = radar.fields['gate_id']['data']
r = radar.range['data']
which_ranges = np.where(r > 20000.0)[0]
which_elevations = np.where(elevation > 1)[0]
above_melting = gate_z > 4000.0
below_melting = gate_z <= 4000.0
# First test for convective sample - filter out
num_heavy_gates_aloft = len(np.where(ref[above_melting] > 40.0)[0])
if(num_heavy_gates_aloft > 40):
return np.nan
snr_rho_filter = np.logical_and(snr > 20, rhohv > 0.98)
ref_filtered = ref[snr_rho_filter]
len_1030 = len(np.where(np.logical_and(ref_filtered >= 10, ref_filtered <= 30))[0])
len_30 = len(np.where(ref_filtered > 30)[0])
if(len_30/len_1030 > 0.2):
return np.nan
zdr = np.ma.masked_where(gate_id != 1, zdr)
ref_filtered = np.ma.masked_where(gate_id != 1, ref)
zdr = zdr[which_elevations, :]
zdr = zdr[:, which_ranges]
ref_filtered = ref_filtered[which_elevations, :]
ref_filtered = ref_filtered[:, which_ranges]
del radar
z_bins = np.linspace(19, 30.5, 7)
mean_zdr = np.zeros(len(z_bins)-1)
for i in range(len(z_bins)-1):
inds = np.where(np.logical_and(ref_filtered >= z_bins[i], ref_filtered <= z_bins[i+1]))
print(inds)
if(len(inds[0]) > 200):
mean_zdr[i] = np.ma.mean(zdr[inds])
else:
mean_zdr[i] = np.nan
return np.nanmean(mean_zdr)
In [6]:
get_zdr_values(file_path)
Out[6]:
In [7]:
theta = np.linspace(0, 90, 120)
zdr05 = zdr_adjusted(0.5, theta)
plt.plot(theta, zdr05)
Out[7]:
In [8]:
radar.fields['gate_id']
Out[8]:
In [9]:
display = pyart.graph.RadarDisplay(radar)
gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above('gate_id', 1)
gatefilter.exclude_equal('gate_id', 0)
display.plot_ppi('differential_reflectivity', 0, cmap='pyart_HomeyerRainbow', vmin=-4, vmax=3, gatefilter=gatefilter)
In [77]:
#jobqueue:
# pbs:
# name: dask-worker
# cores: 36
# memory: 270GB
# processes: 6
# interface: ib0
# local-directory: $localscratch
# queue: high_mem # Can also select batch or gpu_ssd
# project: arm
# walltime: 00:30:00 #Adjust this to job size
# job-extra: ['-W group_list=cades-arm']
cluster = PBSCluster(cores=36, processes=20, threads=1, name='dask-worker', ppn=36, memory='8GB', interface='ib0',
local_directory='$localscratch', queue='batch',
walltime='2:00:00', job_extra=["-W group_list=cades-arm", "-l qos=std"])
cluster.scale(1)
In [78]:
client = Client(cluster)
In [75]:
client
Out[75]:
In [79]:
print(cluster.job_script())
In [154]:
file_list = sorted(glob('/lustre/or-hydra/cades-arm/proj-shared/sgpxsaprcmacsurI5.c1/20171004/*.nc'))
the_bag = db.from_sequence(file_list)
biases = the_bag.map(get_zdr_values).compute()
In [183]:
biases = np.array(biases)
time = [parse_time(name) for name in file_list]
#biases_mean = np.mean(rolling_window(biases, 12), -1)
plt.figure(figsize=(10,5))
plt.plot(time, biases)
plt.plot(time, np.nanmean(biases)*np.ones_like(biases), color='k')
plt.title('ZDR bias daily mean = ' + "%4.2f" % np.nanmean(biases) + ' dB')
plt.xlabel('Time')
plt.ylabel('estimated ZDR bias')
Out[183]:
In [150]:
def rolling_window(a, window):
shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
strides = a.strides + (a.strides[-1],)
return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
#sgpxsaprcmacsurI5.c1.20171004.160022.nc
def parse_time(fname):
year = int(fname[-18:-14])
month = int(fname[-14:-12])
day = int(fname[-12:-10])
hour = int(fname[-9:-7])
minute = int(fname[-7:-5])
second = int(fname[-5:-3])
return datetime.datetime(year, month, day, hour, minute, second)
In [135]:
get_zdr_values(file_list[8])
Out[135]:
In [170]:
# Add a 5.74 dB offset
from copy import deepcopy
radar.fields['corrected_differential_reflectivity'] = deepcopy(radar.fields['differential_reflectivity'])
radar.fields['corrected_reflectivity'] = deepcopy(radar.fields['reflectivity'])
radar.fields['corrected_reflectivity']['data'] = np.ma.masked_where(gatefilter.gate_excluded,
radar.fields['corrected_reflectivity']['data'])
# Get specific differential attenuation
# Need height over 0C isobar
iso0 = np.ma.mean(radar.fields['height']['data'][np.where(np.abs(radar.fields['sounding_temperature']['data']) < 0.1)])
radar.fields['height_over_iso0'] = deepcopy(radar.fields['height'])
radar.fields['height_over_iso0']['data'] -= iso0
spec_at, pia_dict, cor_z, spec_diff_at, pida_dict, cor_zdr = pyart.correct.calculate_attenuation_zphi(
radar, temp_field='sounding_temperature', iso0_field='height_over_iso0', zdr_field='corrected_differential_reflectivity',
pia_field='path_integrated_attenuation', refl_field='corrected_reflectivity', c=0.05, d=1, a_coef=0.313, beta=1)
cor_zdr['data'] += 3.05
radar.add_field('specific_attenuation', spec_at, replace_existing=True)
radar.add_field('path_integrated_attenuation', pia_dict, replace_existing=True)
radar.add_field('corrected_reflectivity', cor_z, replace_existing=True)
radar.add_field('specific_differential_attenuation', spec_diff_at, replace_existing=True)
radar.add_field('path_integrated_differential_attenuation', pida_dict, replace_existing=True)
radar.add_field('corrected_differential_reflectivity', cor_zdr, replace_existing=True)
In [181]:
display = pyart.graph.RadarDisplay(radar)
gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above('gate_id', 1)
gatefilter.exclude_equal('gate_id', 0)
display.plot_ppi('corrected_differential_reflectivity', 0, cmap='pyart_HomeyerRainbow', vmin=-3, vmax=8, gatefilter=gatefilter)
In [182]:
display = pyart.graph.RadarDisplay(radar)
gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above('gate_id', 1)
gatefilter.exclude_equal('gate_id', 0)
display.plot_ppi('differential_reflectivity', 0, cmap='pyart_HomeyerRainbow', vmin=-3, vmax=8, gatefilter=gatefilter)
In [175]:
display = pyart.graph.RadarDisplay(radar)
gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above('gate_id', 1)
gatefilter.exclude_equal('gate_id', 0)
display.plot_ppi('specific_differential_attenuation', 0, cmap='pyart_HomeyerRainbow', vmin=0, vmax=0.1, gatefilter=gatefilter)
In [ ]: