-
Notifications
You must be signed in to change notification settings - Fork 0
/
rca_calc_baseline_hdf5_func.py
149 lines (133 loc) · 6.35 KB
/
rca_calc_baseline_hdf5_func.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
import numpy as np
import matplotlib.pyplot as plt
import pyart
def rca_create_baseline(filename, PCT_on_50, vPCT_on_50, uhPCT_on_50, uvPCT_on_50):
'''Calculates the 95th percentile clutter area reflectivity baseline using clutter map'''
#radar = pyart.io.cfradial.read_cfradial(filename)
radar = pyart.aux_io.read_gamic(filename, file_field_names=True)
date_time = radar.time['units'].replace('seconds since ', '')
# Constrain range between 1 - 5 km
r_start_idx = np.where(radar.range['data'] < 1000.)[0][-1]+1
r_start_idx = 0
r_stop_idx = np.where(radar.range['data'] > 10000.)[0][0]
# Using lowest elevation angle of PPI (0.5 deg)
sweep_start_idx = radar.sweep_start_ray_index['data'][0]
sweep_stop_idx = radar.sweep_end_ray_index['data'][0]+1
# Get variables (only the rays/gates needed)
#zh = radar.fields['reflectivity']['data'][sweep_start_idx:sweep_stop_idx,r_start_idx:r_stop_idx]
zh = radar.fields['Zh']['data'][sweep_start_idx:sweep_stop_idx,r_start_idx:r_stop_idx]
zv = radar.fields['Zv']['data'][sweep_start_idx:sweep_stop_idx,r_start_idx:r_stop_idx]
uzh = radar.fields['UZh']['data'][sweep_start_idx:sweep_stop_idx,r_start_idx:r_stop_idx]
uzv = radar.fields['UZv']['data'][sweep_start_idx:sweep_stop_idx,r_start_idx:r_stop_idx]
r = radar.range['data'][r_start_idx:r_stop_idx]
theta = radar.azimuth['data'][sweep_start_idx:sweep_stop_idx]
# Eliminate duplicate azimuths to maintain a total # of azimuths of 360
if len(theta) > 360:
diff = len(theta) - 360
zh = np.delete(zh,-diff,axis=0)
zv = np.delete(zv,-diff,axis=0)
uzh = np.delete(uzh,-diff,axis=0)
uzv = np.delete(uzv,-diff,axis=0)
theta = np.delete(theta,-diff)
# Arrange/sort azimuths to span 0 to 360 deg. from index 0 to 359
sorted_idx = np.argsort(theta)
zh = zh[sorted_idx]
zv = zv[sorted_idx]
uzh = uzh[sorted_idx]
uzv = uzv[sorted_idx]
theta = theta[sorted_idx]
# Create array to store qualifying reflectivities (fall within PCT_on > 0.5)
zh_car = np.empty((zh.shape))
zh_car[:,:] = np.nan
zv_car = np.empty((zv.shape))
zv_car[:,:] = np.nan
uzh_car = np.empty((uzh.shape))
uzh_car[:,:] = np.nan
uzv_car = np.empty((uzv.shape))
uzv_car[:,:] = np.nan
# H POLARIZATION
# Find and store all reflectivity values that fall within the PCT_on > 0.5 grid gate
for i in range(0,len(PCT_on_50[:,0])):
for j in range(0,len(PCT_on_50[0,:])):
if np.isfinite(PCT_on_50[i,j]):
zh_car[i,j*10-10:j*10] = zh[i,j*10-10:j*10]
# Calculate the PDF of the clutter area reflectivity (CAR)
mask = np.where(np.isfinite(zh_car))
n,bins=np.histogram(zh_car[mask],bins=525,range=(-40.,65.))
n,bins,patches=plt.hist(zh_car[mask],bins=525,range=(-40.,65.))
# Calculate CDF of clutter area reflectivity
cdf = np.cumsum(n)
p = cdf/cdf[-1]*100
# Find coefficients of 13th degree polynomial for CDF
x = np.arange(525)*(1/5)-40
coeff = np.polyfit(p,x,13)
poly_func = np.poly1d(coeff)
#x_poly = np.linspace(p[0],p[-1],105)
#y_poly = poly_func(x_poly)
# Find the value of reflectivity at the 95th percentile of CDF
dbz95 = poly_func(95.)
# UZH POLARIZATION
# Find and store all reflectivity values that fall within the PCT_on > 0.5 grid gate
for i in range(0,len(uhPCT_on_50[:,0])):
for j in range(0,len(uhPCT_on_50[0,:])):
if np.isfinite(uhPCT_on_50[i,j]):
uzh_car[i,j*10-10:j*10] = uzh[i,j*10-10:j*10]
# Calculate the PDF of the clutter area reflectivity (CAR)
mask = np.where(np.isfinite(uzh_car))
uhn,uhbins=np.histogram(uzh_car[mask],bins=525,range=(-40.,65.))
#uhn,uhbins,uhpatches=plt.hist(uzh_car[mask],bins=525,range=(-40.,65.))
# Calculate CDF of clutter area reflectivity
cdf = np.cumsum(uhn)
uhp = cdf/cdf[-1]*100
# Find coefficients of 13th degree polynomial for CDF
x = np.arange(525)*(1/5)-40
coeff = np.polyfit(uhp,x,13)
poly_func = np.poly1d(coeff)
#x_poly = np.linspace(p[0],p[-1],105)
#y_poly = poly_func(x_poly)
# Find the value of reflectivity at the 95th percentile of CDF
dbz95_uh = poly_func(95.)
# V POLARIZATION
# Find and store all reflectivity values that fall within the PCT_on > 0.5 grid gate
for i in range(0,len(vPCT_on_50[:,0])):
for j in range(0,len(vPCT_on_50[0,:])):
if np.isfinite(vPCT_on_50[i,j]):
zv_car[i,j*10-10:j*10] = zv[i,j*10-10:j*10]
# Calculate the PDF of the clutter area reflectivity (CAR)
mask = np.where(np.isfinite(zv_car))
vn,vbins=np.histogram(zv_car[mask],bins=525,range=(-40.,65.))
#vn,vbins,vpatches=plt.hist(zv_car[mask],bins=525,range=(-40.,65.))
# Calculate CDF of clutter area reflectivity
cdf = np.cumsum(vn)
vp = cdf/cdf[-1]*100
# Find coefficients of 13th degree polynomial for CDF
x = np.arange(525)*(1/5)-40
coeff = np.polyfit(vp,x,13)
poly_func = np.poly1d(coeff)
#x_poly = np.linspace(p[0],p[-1],105)
#y_poly = poly_func(x_poly)
# Find the value of reflectivity at the 95th percentile of CDF
dbz95_v = poly_func(95.)
# UV POLARIZATION
# Find and store all reflectivity values that fall within the PCT_on > 0.5 grid gate
for i in range(0,len(uvPCT_on_50[:,0])):
for j in range(0,len(uvPCT_on_50[0,:])):
if np.isfinite(uvPCT_on_50[i,j]):
uzv_car[i,j*10-10:j*10] = uzv[i,j*10-10:j*10]
# Calculate the PDF of the clutter area reflectivity (CAR)
mask = np.where(np.isfinite(uzv_car))
uvn,uvbins=np.histogram(uzv_car[mask],bins=525,range=(-40.,65.))
#uvn,uvbins,uvpatches=plt.hist(uzv_car[mask],bins=525,range=(-40.,65.))
# Calculate CDF of clutter area reflectivity
cdf = np.cumsum(uvn)
uvp = cdf/cdf[-1]*100
# Find coefficients of 13th degree polynomial for CDF
x = np.arange(525)*(1/5)-40
coeff = np.polyfit(uvp,x,13)
poly_func = np.poly1d(coeff)
#x_poly = np.linspace(p[0],p[-1],105)
#y_poly = poly_func(x_poly)
# Find the value of reflectivity at the 95th percentile of CDF
dbz95_uv = poly_func(95.)
del radar
return date_time, n, bins, p, dbz95, vn, vbins, vp, dbz95_v, uhn, uhbins, uhp, dbz95_uh, uvn, uvbins, uvp, dbz95_uv