The idea is to match ground radar (GR) and space-born radar (SR)
measurements in order to create spatially and temporally coicident
samples without interpolation. The procedure had been suggested by
Schwaller and Morris
(2011) and is based on the
adaption by Warren, et. al.
(2017).

The basic principle is illustrated in Fig. 2 of the original paper of
Schwaller and Morris (2011):

image

Quote Warren, et.al: “[…] In this approach, intersections between
indi vidual SR beams and GR elevation sweeps are identified and the
reflectivity values from both instruments are averaged within a spatial
neighborhood around the intersection. Specifically, SR data are averaged
in range over the width of the GR beam at the GR range of the
intersection, while GR data are averaged in the range–azimuth plane
within the footprint of the SR beam. The result is a pair of
reflectivity measurements corresponding to approximately the same volume
of atmosphere. […]”.

The base routines are designed to process one GR sweep at a time. If a
full GR volume with nelev of sweeps is available, you can iterate
over each sweep. In this code, ee is an index that points to one of
the nelev sweeps/elevation angles. Accordingly, a GR data set
will be organised as an array of shape
(nelev_gr,nray_gr,ngate_gr).

A SR data set is typically organised as arrays with dimensions
(nscan_sr,nray_sr,ngate_sr).

Please note: the GR and SR datasets are managed by using convenience
functions from the module file io_func.py located in the same
directory as this notebook.

# Set parameters for this procedurebw_sr=0.71# SR beam widthplatf="gpm"# SR platform/product: one out of ["gpm", "trmm"]zt=sr_pars[platf]["zt"]# SR orbit height (meters)dr_sr=sr_pars[platf]["dr"]# SR gate length (meters)gr_file=sr_pars[platf]["gr_file"]ee=2# Index that points to the GR elevation angle to be used

# number of rays in gr sweepnray_gr=gr_data['nbeam'].astype("i4")[ee]# number of gates in gr beamngate_gr=gr_data['ngate'].astype("i4")[ee]# number of sweepsnelev=gr_data['ntilt']# elevation of sweep (degree)elev_gr=gr_data['elang'][ee]# gate length (meters)dr_gr=gr_data['dr'][ee]# reflectivity array of sweepref_gr=gr_data['refl'][ee]# sweep datetime stampdate_gr=gr_data['sdate'][ee]# range of first gater0_gr=gr_data['r0'][ee]# azimuth angle of first beama0_gr=gr_data['a0'][ee]# Longitude of GRlon0_gr=gr_data['lon']# Latitude of GRlat0_gr=gr_data['lat']# Altitude of GR (meters)alt0_gr=gr_data['alt']# Beam width of GR (degree)bw_gr=1.print(elev_gr,lon0_gr)

# Longitudes of SR scanssr_lon=sr_data['lon']# Latitudes of SR scanssr_lat=sr_data['lat']# Precip flagpflag=sr_data['pflag']# Number of scans on SR datanscan_sr=sr_data['nscan']# Number of rays in one SR scannray_sr=sr_data['nray']# Number of gates in one SR rayngate_sr=sr_data['nbin']

Based on the above criteria (in radar range, precipitating SR profile)
and based on SR elevation angle (with regard to GR).

In [26]:

# First assumption: no valid SR bins (all False)valid=np.asarray(elev_sr,dtype=np.bool)==Falseprint(valid.shape,precip_mask.shape)# SR is inside GR range and is precipitatingiscan=precip_mask.nonzero()[0]iray=precip_mask.nonzero()[1]valid[iscan,iray]=True# SR bins intersect with GR sweepvalid=valid&(elev_sr>=(elev_gr-bw_gr/2.))&(elev_sr<=(elev_gr+bw_gr/2.))# Number of matching SR bins per profilenvalids=np.sum(valid,axis=2)# scan and ray indices for profiles with at least one valid binvscan,vray=np.where(nvalids>0)# number of profiles with at least one valid binnprof=len(vscan)print(vscan.shape)print(valid.shape)

# Lots of containers to store samples (only for one GR sweep angle!)x=np.zeros(nprof)*np.nan# x coordinate of sampley=np.zeros(nprof)*np.nan# y coordinate of samplez=np.zeros(nprof)*np.nan# z coordinate of sampledz=np.zeros(nprof)*np.nan# depth of sampleds=np.zeros(nprof)*np.nan# width of samplers=np.zeros(nprof)*np.nan# range of sample from GRrefsr1=np.zeros(nprof)*np.nan# SR reflectivityrefsr2=np.zeros(nprof)*np.nan# SR reflectivity (S-band, snow)refsr3=np.zeros(nprof)*np.nan# SR reflectivity (S-band, hail)refgr1=np.zeros(nprof)*np.nan# GR reflectivityrefgr2=np.zeros(nprof)*np.nan# GR reflectivity (Ku-band)ntotpr=np.zeros(nprof,dtype="i4")# total number of SR bins in samplenrej1=np.zeros(nprof,dtype="i4")# number of rejected SR bins in samplentotgr=np.zeros(nprof,dtype="i4")# total number of GR bins in samplenrej2=np.zeros(nprof,dtype="i4")# number of rejected GR bins in sampleiref1=np.zeros(nprof)*np.nan# path-integrated SR reflectivityiref2=np.zeros(nprof)*np.nan# path-integrated GR reflectivitystdv1=np.zeros(nprof)*np.nan# std. dev. of SR reflectivity in samplestdv2=np.zeros(nprof)*np.nan# std. dev. of GR reflectivity in samplevolsr=np.zeros(nprof)*np.nan# total volume of SR bins in samplevolgr=np.zeros(nprof)*np.nan# total volume of GR bins in sample

In [41]:

%%time# Loop over relevant SR profilesforii,(ss,rr)inenumerate(zip(vscan,vray)):# Index and count valid bins in each profileip=np.where(valid[ss,rr])[0]numbins=len(ip)ntotpr[ii]=numbinsifnumbins==0:continue# Compute the mean position of these binsx[ii]=np.mean(xyzp_sr[ss,rr,ip,0])y[ii]=np.mean(xyzp_sr[ss,rr,ip,1])z[ii]=np.mean(xyzp_sr[ss,rr,ip,2])# Thickness of the layerdz[ii]=(numbins*dr_sr)/np.cos(np.radians(alpha[ss,rr]))# SR averaging volumevolsr[ii]=np.sum(vol_sr[ss,rr,ip])# Note mean TRMM beam diameterds[ii]=np.radians(bw_sr)*np.mean(((zt-z[ii])/np.cos(np.radians(alpha[ss,rr]))))# Note distance from radars=np.sqrt(x[ii]**2+y[ii]**2)rs[ii]=(re2+z[ii])*np.sin(s/re2)/np.cos(np.radians(elev_gr))# This should not be required because we applied ZonalData### Check that sample is within radar range##if r[ii,jj]+ds[ii,jj]/2. gt rmax then continue## THIS IS THE ORIGINAL IDL CODE - IS THIS A BUG???##ref1[ii,jj]=MEAN(ref_sr1,/nan)##ref3[ii,jj]=MEAN(ref_sr2,/nan)##ref4[ii,jj]=MEAN(ref_sr3,/nan)# Simple linear average of reflectivity# - we can become fancier in the next step# ATTENTION: NEED TO FLIP ARRAYrefsr1[ii]=np.nanmean(ref_sr[ss,rr,ip])refsr2[ii]=np.nanmean(ref_sr_ss[ss,rr,ip])refsr3[ii]=np.nanmean(ref_sr_sh[ss,rr,ip])## Not sure why we need this...### Note the number of rejected bins##nrej1[ii,jj]=ROUND(TOTAL(FINITE(ref_sr1,/nan)))##if FINITE(stdv1[ii,jj]) eq 0 and np-nrej1[ii,jj] gt 1 then STOP# SHOULD WE USE ZONALDATA INSTEAD? COULD BE MORE ACCURATE, BUT ALSO SLOWER# WE COULD BASICALLY START A NEW LOOP HERE AND RUN ZONALDATA BEFORE# Compute the horizontal distance to all the GR binsd=np.sqrt((gr_xyz[...,0]-x[ii])**2+(gr_xyz[...,1]-y[ii])**2)# Find all GR bins within the SR beamaa,bb=np.where(d<=ds[ii]/2.)# Store the number of binsntotgr[ii]=len(aa)iflen(aa)==0:continue# Extract the relevant GR bins# Compute the GR averaging volumevolgr[ii]=np.sum(vol_gr[aa,bb])# Average over those bins that exceed the reflectivity threshold# IDL code does exponential distance and volume weighting# Let's try simple mean first,# THEN ZonalStats!refgr1[ii]=np.nanmean(ref_gr[aa,bb])refgr2[ii]=np.nanmean(ref_gr_ku[aa,bb])