diff --git a/.gitignore b/.gitignore index df140b2..564148d 100755 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,4 @@ *.ipynb -*.nc *.cdf *.png *.tar.gz @@ -7,3 +6,9 @@ *.pyc *.p *.npy +*.log +output/* +other_samples/* +outputfig/* +configtxt/* +nexrad* diff --git a/GeneralFunctions.py b/GeneralFunctions.py index ade2f0e..10af8c0 100644 --- a/GeneralFunctions.py +++ b/GeneralFunctions.py @@ -35,7 +35,8 @@ #import analysis_tools as AT #import lightning_tools as LT -from CSU_RadarTools.csu_radartools import csu_fhc +#from CSU_RadarTools.csu_radartools import csu_fhc +import csu_fhc import general_tools as gentools import RadarConfig from matplotlib.colors import from_levels_and_colors diff --git a/README b/README index 8adff29..207e733 100644 --- a/README +++ b/README @@ -7,5 +7,11 @@ Colorado State University, Dept. Atmos. Sci. bdolan@atmos.colostate.edu May 2017 +Contributions by Anthony Di Stefano +The University of British Columbia, Dept. EOAS +adistefa@eoas.ubc.ca +June 2021 + based on code from Brody Fuchs +Last updated on June 23, 2021 diff --git a/RadarConfig.py b/RadarConfig.py index 401c3b1..022844e 100644 --- a/RadarConfig.py +++ b/RadarConfig.py @@ -12,7 +12,8 @@ class RadarConfig(object): def __init__(self, dz='DZ', zdr='DR', kdp='KD', ldr='LH', rho='RH', hid = 'HID',conv='Con', temp='T', x='x', y='y', z='z', u='U', v='V',rr='RR', w='Wvar',vr='VR',mphys='None',exper = 'Case', - band = 'C',lat_0 = 0,lon_0=90.0,lat_r=None,lon_r=None,lat=None,lon=None,tm = None,radar_name = None): + band = 'C',lat_0 = 0,lon_0=90.0,lat_r=None,lon_r=None,lat=None,lon=None,tm = None,radar_name = None, + color_blind = False): # ******** first the polarimetric stuff ************* self.dz_name = dz self.zdr_name = zdr @@ -56,21 +57,23 @@ def __init__(self, dz='DZ', zdr='DR', kdp='KD', ldr='LH', rho='RH', hid = 'HID', self.species = np.array(['DZ','RN','CR','AG','WS','VI','LDG','HDG','HA','BD']) - self.hid_colors = ['White','LightBlue','MediumBlue','Darkorange','LightPink','Cyan','DarkGray',\ - 'Lime','Yellow','Red','Fuchsia'] + #self.hid_colors = ['White','LightBlue','MediumBlue','Darkorange','LightPink','Cyan','DarkGray',\ + # 'Lime','Yellow','Red','Fuchsia'] + self.hid_colors = ['LightBlue','MediumBlue','Darkorange','LightPink','Cyan','DarkGray',\ + 'Lime','Yellow','Red','Fuchsia'] self.pol_vars = np.array([self.dz_name, self.zdr_name, self.kdp_name, self.ldr_name, self.rho_name, self.hid_name]) self.cs_colors = ['#FFFFFF', 'DodgerBlue', 'Red', 'Khaki'] self.cs_labels = ['', 'Strat', 'Conv', 'Mixed'] - self.set_dbz_colorbar() + self.set_dbz_colorbar(color_blind=color_blind) self.set_hid_colorbar() self.set_cs_colorbar() # Now just set some defaults self.lims = {dz: [0,80], zdr: [-1, 3], kdp: [-0.5, 3], ldr: [-35, -20], rho: [0.95, 1.00], hid: [0, len(self.species)+1],w:[-25,25],vr:[-25,25],self.cs_name:[0,4],self.rr_name:[0.01,150]} self.delta = {dz: 10, zdr: 1, kdp: 1, ldr: 5, rho: 0.005, hid: 1,w:5,vr:5,self.cs_name:1,self.rr_name:10} - self.units = {dz: '(dBZ)', zdr: '(dB)', kdp: '($^{\circ}$/km)', ldr: '(dB)', rho: '', hid: '',w:'m s$^{-1}$',vr:'m s$^{-1}$',self.cs_name:'',self.rr_name:'mm hr$^{-1}$'} + self.units = {dz: '(dBZ)', zdr: '(dB)', kdp: '($^{\circ}$/km)', ldr: '(dB)', rho: '', hid: '',w:'(m s$^{-1}$)',vr:'(m s$^{-1}$)',self.cs_name:'',self.rr_name:'(mm hr$^{-1}$)'} self.names = {dz: 'Z', zdr: 'Z$_{DR}$', kdp: 'K$_{dp}$', ldr: 'LDR', rho: r'$\rho_{hv}$', hid: '',w:'',vr:'V$_r$',self.cs_name:'',self.rr_name:'RR'} self.longnames = {dz: 'Reflectivity', zdr: 'Differntial reflectivity', kdp: 'Specific differential phase',\ ldr: 'Linear depolarization ratio', rho: 'Correlation coefficient', hid: 'Hydrometeor identification',w:'Vertical Velocity',vr:'Radial Velocity',\ @@ -126,13 +129,19 @@ def sav_title(self,tm = None): return extra ############################################################################################################# - def set_dbz_colorbar(self, color_list=None): + def set_dbz_colorbar(self, color_list=None, color_blind=False): if color_list is None: # just use the default here - radarcbar = ['PeachPuff','Aqua','DodgerBlue','MediumBlue','Lime', \ - 'LimeGreen','Green','Yellow','Orange','OrangeRed', \ - 'Red', 'Crimson','Fuchsia','Indigo','DarkCyan','White'] - else: + if color_blind is not True: + radarcbar = ['PeachPuff','Aqua','DodgerBlue','Blue','Lime', \ + 'LimeGreen','Green','Yellow','Orange','DarkOrange','Red', \ + 'Crimson','Fuchsia','Purple','Indigo','MidnightBlue'] + #radarcbar = ['PeachPuff','Aqua','DodgerBlue','MediumBlue','Lime', \ + # 'LimeGreen','Green','Yellow','Orange','OrangeRed','Red', \ + # 'Crimson','Fuchsia','Indigo','DarkCyan','White'] + else: + radarcbar = ['Lavender', 'Thistle', 'Plum', 'MediumPurple', 'CornFlowerBlue', 'SkyBlue', 'PaleTurquoise', 'LightCyan', 'Yellow', 'Gold', 'Orange', 'DarkOrange', 'Chocolate', 'IndianRed', 'FireBrick', 'Maroon'] + else: radarcbar = deepcopy(color_list) temp_cmap = colors.ListedColormap(radarcbar) @@ -151,7 +160,8 @@ def set_hid_colorbar(self, color_list=None): hidcbar = deepcopy(color_list) self.hid_cmap = colors.ListedColormap(hidcbar) - self.boundshid = np.arange(0,12) + #self.boundshid = np.arange(0,12) + self.boundshid = np.arange(self.hid_cmap.N+1) self.normhid = colors.BoundaryNorm(self.boundshid, self.hid_cmap.N) ############################################################################################################# def set_cs_colorbar(self, color_list=None): diff --git a/RadarData.py b/RadarData.py index 1a3974a..61f0d0f 100644 --- a/RadarData.py +++ b/RadarData.py @@ -55,10 +55,10 @@ class RadarData(RadarConfig.RadarConfig): def __init__(self, data,times, ddata = None,dz='DZ', zdr='DR', kdp='KD', ldr='LH', rho='RH', hid='HID',conv='Con', temp='T', x='x', y='y', z='z', u='U', v='V', w='Wvar', rr='RR',vr='VR',lat=None, lon=None, band='C',exper='CASE',lat_r=None,lon_r=None, radar_name= None,mphys=None,dd_data = None,z_thresh=-10.0,cs_z = 2.0,zconv = 41.,zdr_offset=0, remove_diffatt = False,lat_0 = 0.0,lon_0=90.0, - conv_types = ['CONVECTIVE'],strat_types = ['STRATIFORM'],mixed_types = ['UNCERTAIN'],mixr=['qr','qs','qc','qi','qh','qg'],return_scores=False): + conv_types = ['CONVECTIVE'],strat_types = ['STRATIFORM'],mixed_types = ['UNCERTAIN'],mixr=['qr','qs','qc','qi','qh','qg'],return_scores=False,color_blind=False): super(RadarData, self).__init__(dz=dz, zdr=zdr, kdp=kdp, ldr=ldr, rho=rho, hid=hid, conv=conv,temp=temp, x=x, y=y,lat_0=lat_0,lon_0=lon_0,lat_r=lat_r,lon_r=lon_r, - z=z, u=u, v=v, w=w,rr=rr,vr=vr,mphys=mphys,exper=exper,lat=lat,lon=lon,tm = times,radar_name = radar_name) + z=z, u=u, v=v, w=w,rr=rr,vr=vr,mphys=mphys,exper=exper,lat=lat,lon=lon,tm = times,radar_name = radar_name,color_blind=color_blind) # ********** initialize the data ********************* # self.data = {} @@ -557,7 +557,8 @@ def interp_sounding(self): self.gridded_height = np.zeros(self.data[self.dz_name].shape) #print(self.x) for i in range(self.data[self.z_name].shape[0]): - self.gridded_height[:,:,i,...] = self.data[self.z_name][i] + #self.gridded_height[:,:,i,...] = self.data[self.z_name][i] + self.gridded_height[:,i,:,:] = self.data[self.z_name][i] self.T = np.interp(self.gridded_height, self.snd_height, self.snd_temp) @@ -608,8 +609,8 @@ def set_hid(self, band=None, use_temp=False, name='HID',zthresh = -9999.0,return # print('shape holds',np.shape(dzhold)) if use_temp and hasattr(self, 'T'): #print ('Using T!') - tdum = self.T[v,...] -# print('shape tdum',type(tdum)) + #tdum = self.T[v,...] + tdum = self.T[v,:,:,:] #print(type(tdum),'tdum is') #print('T:',np.shape(tdum)) @@ -621,7 +622,9 @@ def set_hid(self, band=None, use_temp=False, name='HID',zthresh = -9999.0,return # scores.append(scoresdum) #hiddum = np.argmax(scoresdum,axis=0)+1 # print(np.shape(tdum),'tdum') - whbad = np.where(np.logical_and(hiddum ==1,tdum <-5.0)) + #whbad = np.where(np.logical_and(hiddum ==1,tdum <-5.0)) + if tdum.any() == None: whbad = np.where(np.logical_and(hiddum == 1,tdum == None)) + else: whbad = np.where(np.logical_and(hiddum == 1,tdum < -5.0)) dzmask = np.where(np.isnan(dzhold)) hiddum[whbad] = -1 hiddum = np.array(hiddum,dtype='float64') @@ -1447,7 +1450,7 @@ def cappi(self, var, z=1.0, xlim=None, ylim=None, ax=None,ts = None, title_flag= # try: # cs = np.squeeze(csdats.sel(x=slice(xmini,xmaxi),y=slice(ymini,ymaxi)).values) # print('cs shape',np.shape(cs)) - ax.contour(xdat, ydat, csvals, levels=[1,2], colors=['k'], linewidths=[3], alpha=0.8,zorder=10) + cb = ax.contour(xdat, ydat, csvals, levels=[1,2], colors=['k'], linewidths=[3], alpha=0.8,zorder=10) # except: # cs = np.squeeze(csdats.sel(x=slice(xmini, xmaxi), y=slice(ymini, ymaxi)).values) # ax.contour(xdat,ydat,csvals,levels = [0,2],colors=['blue','k'],linewidths = [2],alpha = 0.8) @@ -1461,7 +1464,7 @@ def cappi(self, var, z=1.0, xlim=None, ylim=None, ax=None,ts = None, title_flag= if range_lim >= 10: cb_format = '%d' - + ''' if labels: cb = fig.colorbar(dummy, ax=ax, fraction=0.03, pad=0.03, format=cb_format) if var in self.lims.keys(): @@ -1480,26 +1483,33 @@ def cappi(self, var, z=1.0, xlim=None, ylim=None, ax=None,ts = None, title_flag= else: # if this variable is not included in the defaults, have a lot less customization # can get around this with the **kwargs dummy = ax.pcolormesh(self.data[self.x_name], self.data[self.y_name], self.data[var][z_ind,:,:], **kwargs) - cb = fig.colorbar(dummy, ax=ax, fraction=0.03, pad=0.03) if cblabel is not None: cb.set_label(cblabel) - - + ''' + lur,bur,wur,hur = ax.get_position().bounds + cbar_ax_dims = [lur+wur+0.02,bur-0.001,0.03,hur] + cbar_ax = fig.add_axes(cbar_ax_dims) + cbt = fig.colorbar(dummy,cax=cbar_ax) + cbt.ax.tick_params(labelsize=16) + if var.startswith('REF'): labtxt = 'Reflectivity (dBZ)' + elif var.startswith('RRP'): labtxt = 'Rain Rate (mm/hour)' + else: labtxt = var + cbt.set_label(labtxt, fontsize=16, rotation=270, labelpad=20) ####### plotting limits getting set here ###### if self.x_name == 'longitude': #print('setting min and max',xmin,xmax,ymin,ymax) ax.axis([xmin, xmax, ymin, ymax]) -# if labels: -# ax.set_xlabel('Longitude') -# ax.set_ylabel('Latitude') + if labels: + ax.set_xlabel('Longitude') + ax.set_ylabel('Latitude') else: ax.axis([xmini, xmaxi, ymini, ymaxi]) if labels: - ax.set_xlabel('Distance E of radar (km)') - ax.set_ylabel('Distance N of radar (km)') - + ax.set_xlabel('Distance E of radar (km)',fontsize=16) + ax.set_ylabel('Distance N of radar (km)',fontsize=16) + ax.tick_params(axis='both', which='major', labelsize=16) # Now check for the vectors flag, if it's there then plot it over the radar stuff if vectors is not None: @@ -2035,6 +2045,7 @@ def cfad(self, var, value_bins=None, above=None, below=15.0,tspan=None, pick=Non # tei = self.get_ind(te,np.array(self.date)) # # print 'cscfad',cscfad + if cscfad == 'convective': #mask = np.where(self.raintype != 2) mask= np.where(self.raintype != 2) @@ -2050,7 +2061,10 @@ def cfad(self, var, value_bins=None, above=None, below=15.0,tspan=None, pick=Non mask = np.where(self.raintype > 100) # print('entering deep copy') holddat = deepcopy(self.data[var].values) - self.data[var].values[mask] = np.nan + holddat2 = deepcopy(self.data[var].values) + holddat2[mask] = np.nan + #self.data[var].values[mask] = np.nan + self.data[var].values = holddat2 #print('ready to go in loop!') # if left blank, check the whole thing for ivl, vl in (enumerate(tqdm(looped[:-1]))): @@ -2096,7 +2110,7 @@ def cfad(self, var, value_bins=None, above=None, below=15.0,tspan=None, pick=Non ############################################################################################################# def cfad_plot(self, var, nbins=20, ax=None, maxval=10.0, above=None, below=15.0, bins=None, - log=False, pick=None, z_resolution=1.0,levels=None,tspan =None,cont = False,cscfad = False, **kwargs): + log=False, pick=None, z_resolution=1.0,levels=None,tspan =None,cont = False,cscfad = False, cbar=None, ylab=None, **kwargs): from matplotlib.colors import from_levels_and_colors if bins is None: @@ -2104,7 +2118,6 @@ def cfad_plot(self, var, nbins=20, ax=None, maxval=10.0, above=None, below=15.0, else: pass - multiple = np.int(z_resolution/self.dz) # print self.dz # print 'multiple: {}'.format(multiple) @@ -2125,7 +2138,6 @@ def cfad_plot(self, var, nbins=20, ax=None, maxval=10.0, above=None, below=15.0, else: norm = None - # plot the CFAD cfad_ma = np.ma.masked_where(cfad==0, cfad) # print np.max(cfad_ma),var @@ -2138,7 +2150,6 @@ def cfad_plot(self, var, nbins=20, ax=None, maxval=10.0, above=None, below=15.0, #pc = ax.contourf(bins[0:-1],self.data[self.z_name].data[::multiple],(cfad_ma/np.sum(cfad_ma))*100.,levs,color=cols) pc = ax.contourf(bins[0:-1],hts,(cfad_ma),levs,color=cols,cmap=cmap,norm=norm,extend='both') else: - if levels is not None: cmap, norm = from_levels_and_colors([0.02,0.05,0.1,0.2,0.5,1.0,2.0,5.0,10.0,15.0,20.,25.], ['silver','darkgray','slategrey','dimgray','blue','mediumaquamarine','yellow','orange','red','fuchsia','violet']) # mention levels and colors here #print cmap @@ -2149,17 +2160,24 @@ def cfad_plot(self, var, nbins=20, ax=None, maxval=10.0, above=None, below=15.0, # print np.shape(cfad_ma) - cb = fig.colorbar(pc, ax=ax) - cb.set_label('Frequency (%)') - ax.set_ylabel('Height (km MSL)') + if cbar is not None: + cb = fig.colorbar(pc,ax=ax,pad=0.03) + cb.set_label('Frequency (%)',fontsize=22,rotation=270,labelpad=20) + cb.ax.tick_params(labelsize=20) + + if ylab is not None: + ax.set_ylabel('Height (km MSL)',fontsize=22) + ax.tick_params(axis='y',labelsize=20) + else: ax.tick_params(axis='y',labelsize=0) # # try: - ax.set_xlabel('%s %s' %(var, self.units[var])) + ax.set_xlabel('%s %s' %(var, self.units[var]),fontsize=22) + ax.tick_params(axis='x',labelsize=20) # ax.set_title("{d} {r} {v}".format(d=self.date,r=self.radar_name,v=self.longnames[var])) # ax.set_title('%s %s %s CFAD' % (self.print_date(), self.radar_name, self.longnames[var])) # except: # pass - return fig, ax + return fig, ax, pc ############################################################################################################# @@ -2206,9 +2224,11 @@ def plot_2dhist(self, hist,edge,ax=None,cbon = True): cb = ax.contourf(edge[0][:-1],edge[1][:-1],hist.T,norm=colors.Normalize(vmin=0, vmax=np.max(hist)),levels=np.arange(0.01,np.max(hist),0.01)) if cbon == True: #print ' making colorbar' - plt.colorbar(cb,ax=ax) - + col = plt.colorbar(cb,ax=ax) + col.ax.tick_params(labelsize=24) + ax.tick_params(axis='both',labelsize=24) + # This will just look at the whole volume # if above is None: return fig, ax @@ -2505,11 +2525,13 @@ def HID_barplot_colorbar(self, figure, location = [0.9, 0.1, 0.03, 0.8]): scalarMap = plt.cm.ScalarMappable(norm=self.normhid,cmap=self.hid_cmap) axcb = figure.add_axes(location) # x pos, y pos, x width, y width - cb = mpl.colorbar.ColorbarBase(axcb, cmap=self.hid_cmap, norm=self.normhid, boundaries=self.boundshid,\ - orientation = 'vertical') - cb.set_ticks(np.arange(0,11)) + cb = mpl.colorbar.ColorbarBase(axcb, cmap=self.hid_cmap, norm=self.normhid, boundaries=self.boundshid, orientation = 'vertical') + #cb.set_ticks(np.arange(0,10)) + cb.set_ticks(np.arange(len(self.species))+0.5) # need to add a blank at the beginning of species to align labels correctly - labs = np.concatenate((np.array(['']), np.array(self.species))) + #labs = np.concatenate((np.array(['']), np.array(self.species))) + labs = np.array(self.species) + print(labs) cb.set_ticklabels(labs) return cb @@ -2525,9 +2547,9 @@ def plot_hid_cdf(self, data=None, z_resolution=1.0, ax=None, pick=None,cscfad = if ax is None: fig, ax = plt.subplots(1,1) else: - fig = ax.get_figure() + fig = ax.get_figure() - fig.subplots_adjust(left = 0.07, top = 0.93, right = 0.87, bottom = 0.1) + #fig.subplots_adjust(left = 0.07, top = 0.93, right = 0.87, bottom = 0.1) multiple = np.int(z_resolution/self.dz) # if 'd' in self.data[self.z_name].dims: # hgt = self.data[self.z_name].sel(d=0).values @@ -2535,26 +2557,34 @@ def plot_hid_cdf(self, data=None, z_resolution=1.0, ax=None, pick=None,cscfad = # hgt = self.data[self.z_name].valeus print(len(hgt),'heights!') for i, vl in enumerate(np.arange(0, len(hgt), multiple)): - #print vl,i + #print(vl,i) # print self.data[self.z_name].data[vl] #print data[0,:] # print('in plotting cfad',i, vl,np.shape(data),hgt[i])#,np.shape(data[0,:])) - ax.barh(hgt[i], data[0, i], left = 0., edgecolor = 'none', color = self.hid_colors[1]) + ax.barh(hgt[i], data[0, i], left = 0., align = 'center', color = self.hid_colors[0]) for spec in range(1, len(self.species)): # now looping thru the species to make bar plot - #print spec, np.max(data[spec,i]) - - ax.barh(vl, data[spec, i], left = data[spec-1, i], \ - color = self.hid_colors[spec+1], edgecolor = 'none') + #print spec, np.max(data[spec,i]) + ax.barh(hgt[i], data[spec, i], left = data[spec-1, i], color = self.hid_colors[spec], edgecolor = 'none') + #ax.barh(vl, data[spec, i], left = data[spec-1, i], color = self.hid_colors[spec+1], edgecolor = 'none') + ax.set_xlim(0,100) ax.set_xlabel('Cumulative frequency (%)') ax.set_ylabel('Height (km MSL)') # now have to do a custom colorbar? - self.HID_barplot_colorbar(fig) # call separate HID colorbar function for bar plots + + lur,bur,wur,hur = ax.get_position().bounds + cbar_ax_dims = [lur+wur+0.02,bur-0.001,0.03,hur] + #cbar_ax = fig.add_axes(cbar_ax_dims) + #cbt = plt.colorbar(pc,cax=cbar_ax) + #cbt.ax.tick_params(labelsize=16) + #cbt.set_label('Frequency (%)', fontsize=16, rotation=270, labelpad=20) + + self.HID_barplot_colorbar(fig,cbar_ax_dims) # call separate HID colorbar function for bar plots #fig.suptitle('%04d/%02d/%02d - %02d:%02d:%02d %s, cell %d, HID CDF' \ # %(self.year,self.month,self.date,self.hour,self.minute,self.second, \ # self.radar, self.cell_num), fontsize = 14) - ax.set_title('%s %s HID CDF' % (self.print_date(), self.radar_name)) + #ax.set_title('%s %s HID CDF' % (self.print_date(), self.radar_name)) return fig, ax @@ -2730,12 +2760,16 @@ def updraft_width_profile(self, thresh=5.0, temps=np.arange(20,-60,-5),thresh_dz #temps = self.data[self.z_name].data[0,:] # basically just loop thru the Z and get the associated temperature and area - data = self.data[self.w_name].data + #data = self.data[self.w_name].data + #data_dz = self.data[self.dz_name].data + data = self.data[self.w_name].values + data_dz = self.data[self.dz_name].values if thresh_dz == True: - data[self.data[self.dz_name].data < self.z_thresh]=np.nan + data[data_dz < self.z_thresh] = np.nan # print np.shape(data),'ln2075' for iz, z in enumerate(self.data[self.z_name].data): - values_above = np.where(data[iz,...] >= thresh)[0] + #values_above = np.where(data[iz,...] >= thresh)[0] + values_above = np.where(data[:,iz,:,:] >= thresh)[0] num_above = len(values_above) uw[iz] = num_above*self.dx*self.dy/self.ntimes if self.data[self.x_name].units == "[deg]": @@ -2745,12 +2779,18 @@ def updraft_width_profile(self, thresh=5.0, temps=np.arange(20,-60,-5),thresh_dz #print np.shape(uw) #print np.shape(self.T[0,:,0,0]) # now inerpolate this to the temps listed + print(self.T[0,:,0,0]) + print(temps) + print(uw) + self.T = xr.DataArray(data=self.T,dims=['d','z','y','x']) if 'd' in self.T.dims: print('shapes in updraft width',np.shape(uw),np.shape(self.T.sel(x=0,y=0))) f_temp_u = sint.interp1d(self.T.sel(d=0,x=0,y=0), uw, bounds_error=False) else: f_temp_u = sint.interp1d(self.T[:,0,0], uw, bounds_error=False) uwp_interp = f_temp_u(temps) + print(uwp_interp) + input() #return temps, uw return temps,uwp_interp @@ -2772,7 +2812,7 @@ def def_convstrat(self): def calc_timeseries_stats(self,var,ht_lev = 3,thresh=-99.,cs_flag=False,make_zeros=False,areas=False): #First calculate the domain averaged rain rates at a given level. #if self.rr_name is not None: - # rr_timeseries_uncond = rr.mean(dim=['x','y','z'],skipna=True) + # rr_timeseries_uncond = rr.mean(dim=['x','y','z'],skipna=True) data = deepcopy(self.data[var].sel(z=slice(ht_lev,ht_lev+1))) whbad= np.where(data.values0) - if areas==False: + if areas==False: datstrat_ts= datstrat.mean(dim=['z','y','x'],skipna=True) datconv_ts= datconv.mean(dim=['z','y','x'],skipna=True) datall_ts= datall.mean(dim=['z','y','x'],skipna=True) @@ -2798,7 +2838,6 @@ def calc_timeseries_stats(self,var,ht_lev = 3,thresh=-99.,cs_flag=False,make_zer datstrat_ts = cssum.where(cssum==1).count(dim=['y','x']) datconv_ts = cssum.where(cssum==2).count(dim=['y','x']) datall_ts = cssum.where(cssum>0).count(dim=['y','x']) - return datstrat_ts,datconv_ts,datall_ts else: @@ -3021,4 +3060,4 @@ def composite(self, var,ts=0,map_on = True,res='10m'): gl.xlabels_top = False gl.ylabels_right = False - # ax.set_title('MC3E CSAPR {d: \ No newline at end of file + # ax.set_title('MC3E CSAPR {d: diff --git a/accessing_so_files.txt b/accessing_so_files.txt new file mode 100644 index 0000000..765f752 --- /dev/null +++ b/accessing_so_files.txt @@ -0,0 +1,3 @@ +When running on MAC or Linux OS, perform the following commands before running `python run_ipolarris.py my_config.txt`: +1) Put /usr/bin ahead of your "default" executable directories (i.e. before /opt/local/bin) in $PATH. It does not need to be ahead of your custom executable directories (i.e. ~/../anaconda3/bin). +2) Run: f2py -c calc_kdp_ray_fir.f -m calc_kdp_ray_fir. This will allow you to use the Fortran compiler in /usr/bin to convert your .f file into a readable .so file for your MAC or Linux OS. diff --git a/calc_kdp_ray_fir.cpython-39-darwin.so b/calc_kdp_ray_fir.cpython-39-darwin.so new file mode 100755 index 0000000..f3345c9 Binary files /dev/null and b/calc_kdp_ray_fir.cpython-39-darwin.so differ diff --git a/calc_kdp_ray_fir.cpython-39-x86_64-linux-gnu.so b/calc_kdp_ray_fir.cpython-39-x86_64-linux-gnu.so new file mode 100755 index 0000000..3a3c698 Binary files /dev/null and b/calc_kdp_ray_fir.cpython-39-x86_64-linux-gnu.so differ diff --git a/calc_kdp_ray_fir.so b/calc_kdp_ray_fir.so index 2622771..907725e 100755 Binary files a/calc_kdp_ray_fir.so and b/calc_kdp_ray_fir.so differ diff --git a/doppfiles.txt b/doppfiles.txt new file mode 100644 index 0000000..de1510e --- /dev/null +++ b/doppfiles.txt @@ -0,0 +1,2 @@ +samples/KLGX_KATX_20151203_010200_V06_dopvel_gridded.nc +samples/KLGX_KATX_20151203_010607_V06_dopvel_gridded.nc diff --git a/dpolfiles.txt b/dpolfiles.txt new file mode 100644 index 0000000..af9d0ae --- /dev/null +++ b/dpolfiles.txt @@ -0,0 +1,2 @@ +samples/KLGX_20151203_010200_V06.nc +samples/KLGX_20151203_010607_V06.nc diff --git a/env.yml b/env.yml new file mode 100644 index 0000000..8dabfc3 --- /dev/null +++ b/env.yml @@ -0,0 +1,21 @@ +name: pol +channels: +- defaults +- conda-forge +dependencies: +- pandas +- dill +- netcdf4 +- xarray +- python +- matplotlib +- pyproj +- tqdm +- libgfortran +- cartopy +- pip +- toolz +- pip: + - dask + - click + - scipy diff --git a/my_config.txt b/my_config.txt new file mode 100644 index 0000000..d968e9c --- /dev/null +++ b/my_config.txt @@ -0,0 +1,189 @@ +#################### MY_CONFIG.TXT ################# + +#------------------- +#------ Input ------ +#------------------- + +sdatetime == '20151203_0100' == # Start time of analysis of interest +sdatetime_format == %Y%m%d_%H%M == # Start time format +edatetime == '20151203_0200' == # End time of analysis of interest +edatetime_format == %Y%m%d_%H%M == # End time format + +date == '20151203' == # Date of the case +extra == '' == # This will be appended to all output files as an identifier. +extrax == '' == # This will be appended to all output files as an identifier. +rfiles == './dpolfiles.txt' == # Input directory of radar files to read in +rdate_format == %Y%m%d_%H%M%S == # Format for the radar file date +rdend == 20 == # End of date timestamp in radar filename +rdstart == 5 == # Start of date timestamp in radar filename +type == 'obs' == # Type of input data: 'obs' OR 'wrf' (obs + simulated) + +#-------------------- +#------ Output ------ +#-------------------- + +image_dir == './output/' == # Output figure directory +ptype == 'png' == # Output figure file extenstion (i.e. .png, .jpg, ...) +cb_friendly == True == # Use color-blindness palette for output (i.e. reflectivity) + +#-------------------------- +#------ Output Flags ------ +#-------------------------- + +# Base Plots +compo_ref == False == # FA) Turn on plotting for spatial composite reflectivity plotting +cappi_ref == False == # FB) Turn on plotting for reflectivity CAPPIs +cappi_rr == False == # FC) Turn on plotting for rain rate CAPPIs +rrstats_txt == False == # FD) +rrhist_txt == False == # FE) +rrstats_areas_txt == False == # FF) +rr_timeseries == False == # FG) +vv_profiles == False == # FH) +percentiles_txt == False == # FI) +vert_ref == False == # FJ) +refcfad == False == # FK) + +# Advanced Plots +plot_int == False == # FA) Turn on plotting of integrated parameters over the whole time frame +plot_cs == False == # FB) Plot seperate Convective and stratiform CFADs +####The following are for EACH radar / model time specified in the radar_files list. +cfad_mpanel_flag == False == # FC) 4 panel CFADs of Dbz, Zdr, Kdp and W +hid_cfad_flag == False == # FD) HID CFAD +joint_flag == False == # FE) Z vs. Zdr +cfad_individ_flag == True == # FF) Make separate images for dbz, Zdr, Kdp, W and Rho +hid_prof == False == # FG) Profile of HID species with height +up_width == False == # FJ) Make a plot of the vertical velocity profile with temperature. +tms_on == True == # +cappi_multi == False == # FM) Make cappi cross section for all times. change parameters in plot_driver.py +xsec_multi == True == # FN) make RHI cross sections for all times. change parmaters in plot_driver.py +qr_cappi == False == # FK) Make cappi cross section of mixing ratios. change parameters in plot_driver.py (only valid for model) +qr_rhi == False == # FL) Make rhis of the mixing ratios (only valid for model) + +ks == ['cfad_mpanel_flag','hid_cfad_flag','joint_flag','cfad_individ_flag','hid_prof','cappi_multi','xsec_multi','up_width','qr_cappi','qr_rhi'] == # FP) list of flags + +#------------------------ +#------ Radar Info ------ +#------------------------ + +alt == 354 == # Altitude of the radar +band == S == # Radar band: X, C OR S. Note: needs to be capital letter. +exper == 'KLGX' == # Radar location +lat == 47.116806 == # Latitude of the radar station +lon == -124.10625 == # Longitude of the radar station +radarname == S-band == # + +#------------------------ +#------ Radar File ------ +#------------------------ + +dz_name == REF == # Name of the reflectivity field +dr_name == ZDR == # Name of the differential reflectivity field +kd_name == KDP == # Name of the Kdp field +rh_name == RHO == # Name of the RhoHV field +rr_name == RRP == # Name of the Rain rate / precipitation field +vr_name == VEL == # Name of radial velocity field + +#----------------------------- +#------ WRF Output File ------ +#----------------------------- + +mphys == obs == # Type of microphysics used in model: 'obs' OR '' if type = 'wrf' + +#-------------------------- +#------ Doppler Info ------ +#-------------------------- + +dd_on == True == # Set this to loop through separate dual-Doppler files (e.g. for obs) + +dfiles == ./doppfiles.txt == # Location of dual-Doppler files +ddate_format == %Y%m%d_%H%M%S == # Format for the dual-Doppler file date +ddstart == 10 == # Offset for dual-Doppler timestamp in filename +ddend == 25 == # Offset for dual-Doppler timestamp in filename + +#--------------------------- +#------ Doppler File ------- +#--------------------------- + +convname == None == # Name of +uname == eastward_wind == # Name of the zonal wind field +vname == northward_wind == # Name of the meridional wind field +wname == upward_air_velocity == # Name of the vertical wind field +xname == x == # Name of the zonal directional variable +yname == y == # Name of the meridional directional variables +zname == z == # Name of the vertical level field +xlim == [-100,100] ==#lat / lon zoom for CAPPIS +ylim == [-100,100] ==# lat / lon zoom for CAPPIS +y==0.5 ==#Lat of the y cross-section for RHI +#xlim == [-127.0,-121.0] ==#lat / lon zoom for CAPPIS +#ylim == [45.0,49.0] ==# lat / lon zoom for CAPPIS +#y==35.7 ==#Lat of the y cross-section for RHI +z == 2.5 ==#Height of the CAPPIs + +#--------------------------- +#------ Sounding Info ------ +#--------------------------- + +snd_on == True == # Set this to loop through separate dual-Doppler files (e.g. for obs) + +ad == config['date']+'_' == # Extra characters +sdate_format == %Y%m%d_%H%M%S == # Sounding file date format +sfiles == ./soundfiles.txt == # Sounding file directory +sdstart == 4 == # Offset for date timestamp in radar filename +sdend == 19 == # Number of characters in timestamp on radar file +sstat == UIL == # Sounding station identification + +#------------------------------------ +#------ Flags if type == 'wrf' ------ +#------------------------------------ +convert_Tk_Tc == False == # Convert temperature in K to deg C + + +###############Set up variable names and how to read the data############### +t_name == T ==#Name of the temperature field +latname == lat == +lonname == lon == +pol_on == True ==#Calculate the pol data such as HID and RR +z_resolution == 1.0 ==#Vertical resolution for CFADs. If comparing 2, they need to be the same. +zthresh == -10. ==#Threshold for good data +wthresh == 5. ==#Threshold for 'updraft' statistics. +trange == np.arange(20,-60,-5) ==#Range of thresholds for plotting temperatures. +cs_z == 2.0 ==#Level to determine Convective / stratiform designation. +zconv == 40 ==#Zconv threshold in raintyping algorithm. +conv_types == ['ISO_CONV_CORE','CONVECTIVE','ISO_CS_CORE'] ==#Which Powell et al. types to consider in convective CFADS +strat_types == ['WEAK_ECHO','STRATIFORM','ISO_CONV_FRINGE'] ==#Which Powell et al. types to consider in stratiform CFADS +mixed_types == ['UNCERTAIN'] ==#Which types to not include in either convective or stratiform but will be considered in totals). +zdr_offset == 0.6 ==#Add any Zdr offset here. Value will be SUBTRACTED from the zdr values. +mask_model == False == +drop_vars == False == +# +#######Set up some variables related to the observations ################### +## lat/lon of the radar +removediffatt == True ==#Remove differential attenuation by Zdr < -1 and dBZ < 35. +# +############Select the types of plots to see on the output########################## +#############Set up some plotting configurations ######################### +wbins == np.arange(-25,26,0.5) ==#Histogram bins for vertical velocity +dzbins == np.arange(-10,60,1) ==#Histogram bins for reflectivity +drbins == np.arange(-2,6,0.05) ==#Histogram bins for differential reflectivity +kdbins == np.arange(-2,6,0.05) ==#Histogram bins for specific differential phase +rrbins == np.logspace(0.01,100.01,30) ==#Histogram bins for rain rates +# +###Set up the types of hid to group together for water, graupel, hail and snow. +hidwater == [1,2,10] ==#Group drizzle, rain and big drops +hidgraup == [7,8] ==#Group low and high density graupel +hidhail == [9] ==#Hail +hidsnow ==[3,4,5,6] ==#Group ice crystals, snow, wet snow and VI +# +# +####Set up some specifics for the cross-sections.#################### +cvectors == [None,None,None,None,None,True] ==#Turn on vectors in the plots. +rvectors == [None,None,None,None,None,True] ==#Turn on vectors in the plots. +skip == 2 ==#Number of vectors to skip for the vector plots. +mix_vars ==['qc','qr','qg','qi','qh',config['rr_name'],config['vr_name'],'HID'] ==#Mixing ratios from model to plot. +rhi_vars ==['HID',config['dz_name'],config['dr_name'],config['kd_name'],config['rh_name'],config['wname']] == #Names of vars for RHI plots +cappi_vars == ['HID',config['dz_name'],config['dr_name'],config['kd_name'],config['rh_name'],config['vr_name']] == #Names of vars for CAPPI plots +comb_vicr == True ==# Combine VI with CR for plotting. +cappi_contours == ['CS',None,None,None,None,None] ==#What contours to apply to the CAPPI images. +cappi_vectres == 5 ==#Defined the vector skip for cappi plots. +rhi_vectres == [6,2] ==#Defines the [x,z] skip for rhi plots +# diff --git a/plot_driver.py b/plot_driver.py index 0c8b010..5a888d0 100644 --- a/plot_driver.py +++ b/plot_driver.py @@ -17,7 +17,8 @@ import RadarData import GeneralFunctions as GF from matplotlib import colors -plt.style.use('presentation') +plt.style.use('presentation.mplstyle') +#plt.style.use('default') from matplotlib.dates import DateFormatter,HourLocator dayFormatter = DateFormatter('%H%M') # e.g., 12 @@ -760,7 +761,13 @@ def plot_upstat(dat1,dat2,config,typ='hid',n1 = None,n2 = None): def make_single_pplots(rdat,flags,config,y=None): - print ('in make_singl_pplots') + + print('\n#######################################') + print('####### Starting plot_driver.py #######') + print('#######################################\n') + + print('ADVANCED PLOTTING: plot_driver.make_single_pplots') + print('Initiating...') tspan= [rdat.date[0],rdat.date[-1]] tms = np.array(rdat.date) #print('DATES',np.array(rdat.date)) @@ -778,21 +785,32 @@ def make_single_pplots(rdat,flags,config,y=None): if flags['cfad_mpanel_flag'] == True: print ('Working on Cfad mpanel') - fig, ax = plt.subplots(2,2,figsize=(18,12)) + fig, ax = plt.subplots(2,2,figsize=(16,12),gridspec_kw={'wspace': 0.05, 'top': 1., 'bottom': 0., 'left': 0., 'right': 1.}) axf = ax.flatten() if config['wname'] in rdat.data.variables.keys(): - dum =rdat.cfad_plot(rdat.w_name,ax = axf[0],bins=config['wbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan) + dum =rdat.cfad_plot(rdat.w_name,ax = axf[0],bins=config['wbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan,ylab=True) dum =rdat.cfad_plot(rdat.dz_name,ax = axf[1],bins=config['dzbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan) - dum =rdat.cfad_plot(rdat.zdr_name,ax= axf[2],bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan) + dum =rdat.cfad_plot(rdat.zdr_name,ax= axf[2],bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan, ylab=True) dum =rdat.cfad_plot(rdat.kdp_name,ax = axf[3],bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan) - plt.tight_layout() + #plt.tight_layout() # print "{s:%Y%m%d%H%M%S}".format(s=ts[0]) - plt.savefig('{d}{p}_CFAD_4panel_{s:%Y%m%d%H%M%S}_{r}_{m}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,m=rdat.mphys,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + + lur1,bur1,wur1,hur1 = axf[1].get_position().bounds + lur2,bur2,wur2,hur2 = axf[-1].get_position().bounds + cbar_ax_dims = [lur2+wur2+0.02,bur2-0.001,0.03,bur1+hur1] + cbar_ax = fig.add_axes(cbar_ax_dims) + cbt = plt.colorbar(dum[-1],cax=cbar_ax) + cbt.ax.tick_params(labelsize=20) + cbt.set_label('Frequency (%)', fontsize=22, rotation=270, labelpad=20) + + axf[0].text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=axf[0].transAxes) # (a) Top-left + + plt.savefig('{d}{p}_CFAD_4panel_{s:%Y-%m-%d_%H%M%S}_{r}_{m}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,m=rdat.mphys,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() - + if config['plot_cs'] == True: - fig, ax = plt.subplots(2,2,figsize=(18,12)) + fig, ax = plt.subplots(2,2,figsize=(18,12),constrained_layout=True) axf = ax.flatten() if config['wname'] in rdat.data.variables.keys(): @@ -824,82 +842,91 @@ def make_single_pplots(rdat,flags,config,y=None): plt.clf() if flags['cfad_individ_flag'] == True: - fig, ax = plt.subplots(1,1,figsize=(18,12)) # axf = ax.flatten() if config['wname'] in rdat.data.variables.keys(): - rdat.cfad_plot(rdat.w_name,ax = ax,bins=config['wbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan) - plt.tight_layout() - plt.savefig('{d}{p}_CFAD_W_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + fig, ax = plt.subplots(1,1,figsize=(14,12)) + rdat.cfad_plot(rdat.w_name,ax = ax,bins=config['wbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan,cbar=True,ylab=True) + #plt.tight_layout() + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) + plt.savefig('{d}{p}_CFAD_W_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() - - fig, ax = plt.subplots(1,1,figsize=(18,12)) - rdat.cfad_plot(rdat.dz_name,ax = ax,bins=config['dzbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan) - plt.tight_layout() - plt.savefig('{d}{p}_CFAD_dBZ_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + fig, ax = plt.subplots(1,1,figsize=(14,12)) + rdat.cfad_plot(rdat.dz_name,ax = ax,bins=config['dzbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan,cbar=True,ylab=True) + #plt.tight_layout() + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) + plt.savefig('{d}{p}_CFAD_dBZ_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() - fig, ax = plt.subplots(1,1,figsize=(18,12)) - rdat.cfad_plot(rdat.zdr_name,ax= ax,bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan) - plt.tight_layout() - plt.savefig('{d}{p}_CFAD_Zdr_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + fig, ax = plt.subplots(1,1,figsize=(14,12)) + rdat.cfad_plot(rdat.zdr_name,ax= ax,bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan= tspan,cbar=True,ylab=True) + #plt.tight_layout() + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) + plt.savefig('{d}{p}_CFAD_Zdr_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() - fig, ax = plt.subplots(1,1,figsize=(18,12)) - rdat.cfad_plot(rdat.kdp_name,ax = ax,bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan) - plt.savefig('{d}{p}_CFAD_Kdp_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) - plt.tight_layout() + fig, ax = plt.subplots(1,1,figsize=(14,12)) + rdat.cfad_plot(rdat.kdp_name,ax = ax,bins=config['drbins'],z_resolution=config['z_resolution'],levels='levs',tspan = tspan,cbar=True,ylab=True) + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) + plt.savefig('{d}{p}_CFAD_Kdp_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') + #plt.tight_layout() plt.clf() - fig, ax = plt.subplots(1,1,figsize=(18,12)) - rdat.cfad_plot(rdat.rho_name,ax = ax,z_resolution=config['z_resolution'],levels='levs',tspan = tspan) - plt.tight_layout() - plt.savefig('{d}{p}_CFAD_RHO_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + fig, ax = plt.subplots(1,1,figsize=(14,12)) + rdat.cfad_plot(rdat.rho_name,ax = ax,z_resolution=config['z_resolution'],levels='levs',tspan = tspan,cbar=True,ylab=True) + #plt.tight_layout() + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) + plt.savefig('{d}{p}_CFAD_RHO_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() if flags['hid_cfad_flag'] == True: fig, ax = rdat.plot_hid_cdf() - plt.savefig('{d}{p}_CFAD_HID_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=12, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + + plt.savefig('{d}{p}_CFAD_HID_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() if flags['joint_flag'] == True: - fig, ax = plt.subplots(2,2,figsize=(12,12)) + fig, ax = plt.subplots(2,2,figsize=(16,14),gridspec_kw={'wspace': 0.25, 'hspace': 0.2, 'top': 1., 'bottom': 0., 'left': 0., 'right': 1.}) axf = ax.flatten() zzdr_wrf,ed = rdat.hist2d(varx=rdat.dz_name,vary=rdat.zdr_name,binsx=config['dzbins'],binsy=config['drbins']) rdat.plot_2dhist(zzdr_wrf,ed,ax=axf[0]) - axf[0].set_xlabel('Zdr') - axf[0].set_ylabel('dBZ') - axf[0].set_title(title_string) + axf[0].set_xlabel(rdat.zdr_name+' '+rdat.units[rdat.zdr_name],fontsize=26) + axf[0].set_ylabel(rdat.dz_name+' '+rdat.units[rdat.dz_name],fontsize=26,labelpad=0) + #axf[0].set_title(title_string) + + axf[0].text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=12, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left zkdp_wrf,edk = rdat.hist2d(varx=rdat.dz_name,vary=rdat.kdp_name,binsx=config['dzbins'],binsy=config['kdbins']) rdat.plot_2dhist(zkdp_wrf,edk,ax=axf[1]) - axf[1].set_title(title_string) - axf[1].set_xlabel('Kdp') - axf[1].set_ylabel('dBZ') - + #axf[1].set_title(title_string) + axf[1].set_xlabel(rdat.kdp_name+' '+rdat.units[rdat.kdp_name],fontsize=26) + axf[1].set_ylabel(rdat.dz_name+' '+rdat.units[rdat.dz_name],fontsize=26,labelpad=0) zw_wrf,edw = rdat.hist2d(varx=rdat.dz_name,vary=rdat.w_name,binsx=config['dzbins'],binsy=config['wbins']) rdat.plot_2dhist(zw_wrf,edw,ax=axf[2]) - axf[2].set_title(title_string) - axf[2].set_xlabel('W') - axf[2].set_ylabel('dBZ') + #axf[2].set_title(title_string) + axf[2].set_xlabel('W '+rdat.units[rdat.w_name],fontsize=26) + axf[2].set_ylabel(rdat.dz_name+' '+rdat.units[rdat.dz_name],fontsize=26,labelpad=0) zr_wrf,edr = rdat.hist2d(varx=rdat.rr_name,vary=rdat.w_name,binsx=config['rrbins'],binsy=config['wbins'],xthr=0.00000) - cb6 = rdat.plot_2dhist(zr_wrf,edr,ax=axf[3],cbon=True) - axf[3].set_title(title_string) - axf[3].set_xlabel('W') - axf[3].set_ylabel(rdat.rr_name) + cb6 = rdat.plot_2dhist(zr_wrf,edr,ax=axf[3]) + #axf[3].set_title(title_string) + axf[3].set_xlabel('W '+rdat.units[rdat.w_name],fontsize=26) + axf[3].set_ylabel(rdat.rr_name+' '+rdat.units[rdat.rr_name],fontsize=26,labelpad=10) axf[3].set_ylim(0,50) - plt.savefig('{d}{p}_2dPDF_4panel_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + plt.savefig('{d}{p}_2dPDF_4panel_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() if flags['hid_prof'] == True: + fig, ax = plt.subplots(1,1,figsize=(18,12)) + hts, mwrf_water_vert = rdat.hid_vertical_fraction(config['hidwater'],z_resolution =config['z_resolution']) hts, mwrf_graup_vert = rdat.hid_vertical_fraction(config['hidgraup'],z_resolution =config['z_resolution']) hts, mwrf_hail_vert = rdat.hid_vertical_fraction(config['hidhail'],z_resolution =config['z_resolution']) @@ -909,33 +936,46 @@ def make_single_pplots(rdat,flags,config,y=None): plt.plot(mwrf_graup_vert,hts,color='g',label='graupel',lw=lw) plt.plot(mwrf_hail_vert,hts,color='r',label='hail',lw=lw) plt.plot(mwrf_snow_vert,hts,color = 'yellow',label='snow',lw=lw) - plt.xlabel('Frequency (%)') - plt.ylabel('Height (km)') - plt.title(title_string) - plt.legend(loc = 'best') - plt.savefig('{d}{p}_HID_prof_{s:%Y%m%d%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=300) + ax.tick_params(axis='both',labelsize=22) + plt.xlabel('Frequency (%)',fontsize=24) + plt.ylabel('Height (km)',fontsize=24) + #plt.title(title_string) + plt.legend(loc='best',fontsize=22) + + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=24, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + plt.savefig('{d}{p}_HID_prof_{s:%Y-%m-%d_%H%M%S}_{r}_{x}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype']),dpi=400,bbox_inches='tight') plt.clf() + if config['wname'] in rdat.data.variables.keys(): + if flags['up_width'] == True: tmp, m_warea_wrf = rdat.updraft_width_profile(thresh_dz=True) #print np.max(m_warea_wrf) + + fig, ax = plt.subplots(1,1,figsize=(18,12)) + plt.plot(m_warea_wrf,tmp,color='k',label='Obs',lw=5) plt.ylim(20,-60) - plt.xlabel('Updraft Width (km$^2$)') - plt.ylabel('Temperature (deg C)') - plt.title(title_string) - plt.savefig('{d}{p}_upwidth_{s:%Y%m%d%H%M%S}_{r}_{x}_{y}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],y=config['y']),dpi=300) + plt.xlabel('Updraft Width (km$^2$)',fontsize=24) + plt.ylabel('Temperature (deg C)',fontsize=24) + ax.tick_params(axis='both',labelsize=22) + #plt.title(title_string) + + ax.text(0, 1, '{e} {r}'.format(e=rdat.exper,r=rdat.radar_name), horizontalalignment='left', verticalalignment='bottom', size=26, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + + plt.savefig('{d}{p}_upwidth_{s:%Y-%m-%d_%H%M%S}_{r}_{x}_{y}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=tstart,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],y=config['y']),dpi=400,bbox_inches='tight') plt.clf() # for k in flags.keys(): # print('flag keys:',k,flags[k]) - for ts in tms: - if flags['all_cappi']== True: - #z=2.0 - #print xlim -# print (config['cappi_vars']) -# print (config['cappi_multi']) + if config['tms_on'] is True: + + for ts in tms: + #z=2.0 + #print xlim + # print (config['cappi_vars']) + # print (config['cappi_multi']) if config['cappi_multi'] is True: #print('cappi mulit, vars',config['cappi_vars']) #print config['cappi_vectres'],eval(config['cvectors']),eval(config['cappi_contours']),config['ylim'],config['xlim'],config['z'],rdat.date,eval(config['cappi_vars']) @@ -963,10 +1003,8 @@ def make_single_pplots(rdat,flags,config,y=None): #label_subplots(fig,yoff=0.01,xoff=0.01,size=16,nlabels=1) plt.savefig('{d}{p}_polcappi_{v}_{s:%Y%m%d%H%M%S}_{r}_{x}_{z}km.{t}'.format(d=config['image_dir'],v=v,p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],z=config['z']),dpi=300) plt.clf() - - if flags['all_xsec']== True: - #y=-12.5 - print ("xsec_multi") + + #y=-12.5 if config['xsec_multi'] == True: if rdat.w_name is not None: fig = rdat.xsec_multiplot(ts=ts,y=config['y'],vectors=eval(config['rvectors']),res = config['rhi_vectres'],xlim=config['xlim'],varlist=eval(config['rhi_vars'])) @@ -978,7 +1016,7 @@ def make_single_pplots(rdat,flags,config,y=None): yof = 0.01 else: yof=-0.02 - + #plt.tight_layout() label_subplots(fig,yoff=yof,xoff=0.01,size=16,nlabels=nvars) @@ -993,19 +1031,18 @@ def make_single_pplots(rdat,flags,config,y=None): plt.savefig('{d}{p}_polrhi_{v}_{s:%Y%m%d%H%M%S}_{r}_{x}_{y}.{t}'.format(d=config['image_dir'],v=v,p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],y=config['y']),dpi=300) plt.clf() - - if flags['qr_cappi'] == True: - print ("qr_cappi") - fig = rdat.cappi_multiplot(z=config['z'],ts=ts,xlim=config['xlim'],ylim=config['ylim'],varlist=eval(config['mix_vars'])) - plt.savefig('{d}{p}_qcappi_6panel_{s:%Y%m%d%H%M%S}_{r}_{x}_{z}km.{t}'.format(d=config['image_dir'],p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],z=config['z']),dpi=300) - plt.clf() - - if flags['qr_rhi'] == True: - print ("qr_rhi") - fig = rdat.xsec_multiplot(ts=ts,y=config['y'],xlim=config['xlim'],varlist=eval(config['mix_vars'])) - plt.savefig('{d}{p}_qrhi_6panel_{s:%Y%m%d%H%M%S}_{r}_{x}_{y}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],y=config['y']),dpi=300) - plt.clf() - + if flags['qr_cappi'] == True: + print ("qr_cappi") + fig = rdat.cappi_multiplot(z=config['z'],ts=ts,xlim=config['xlim'],ylim=config['ylim'],varlist=eval(config['mix_vars'])) + plt.savefig('{d}{p}_qcappi_6panel_{s:%Y%m%d%H%M%S}_{r}_{x}_{z}km.{t}'.format(d=config['image_dir'],p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],z=config['z']),dpi=300) + plt.clf() + + if flags['qr_rhi'] == True: + print ("qr_rhi") + fig = rdat.xsec_multiplot(ts=ts,y=config['y'],xlim=config['xlim'],varlist=eval(config['mix_vars'])) + plt.savefig('{d}{p}_qrhi_6panel_{s:%Y%m%d%H%M%S}_{r}_{x}_{y}.{t}'.format(d=config['image_dir'],p=rdat.exper,s=ts,r=rdat.radar_name,x=config['extrax'],t=config['ptype'],y=config['y']),dpi=300) + plt.clf() + def subset_convstrat(data,rdata,zlev=1): cssum =(rdata.data[rdata.cs_name].max(dim='z')) stratsub=data.sel(z=slice(zlev,zlev+1)).where(cssum==1) @@ -1039,15 +1076,16 @@ def plot_timeseries(data,tm,ax,ls = '-',cs=False,rdata=None,thresh=-50,typ='',zl ax.plot(np.array(tm),cdat.mean(dim=['z','y','x'],skipna=True),color='r',label='Conv {e}'.format(e=typ),ls=ls) ax.plot(np.array(tm),sdat.mean(dim=['z','y','x'],skipna=True),color='b',label='strat {e}'.format(e=typ),ls=ls) - ax.legend(loc='best') + ax.legend(loc='best',fontsize=14) else: ax.plot(np.array(tm),data.where(data>thresh).mean(dim=['z','y','x'],skipna=True),color='k',label='Total') - ax.xaxis.set_major_formatter(hourFormatter) - ax.xaxis.set_major_locator(HourLocator(interval=1)) - d=plt.setp(plt.gca().get_xticklabels(), rotation=45, horizontalalignment='right') - ax.set_xlabel('Time (UTC)') + #ax.xaxis.set_major_formatter(hourFormatter) + #ax.xaxis.set_major_locator(HourLocator(interval=1)) + #d=plt.setp(plt.gca().get_xticklabels(), rotation=45, horizontalalignment='right') + ax.tick_params(axis='both',labelsize=14) + ax.set_xlabel('Time (UTC)',fontsize=16) return ax#,adat,cdat,sdat @@ -1075,7 +1113,7 @@ def plot_quartiles(data,q1,q2,q3,z,ax,c1='goldenrod',c2='r',c3='k',split_updn=Fa ax.plot(wup50,zdat,color=c2,label='50th {e}'.format(e=typ),ls=ls) ax.plot(wup90,zdat,color=c1,label='90th {e}'.format(e=typ),ls=ls) ax.plot(wup99,zdat,color=c3,label='99th {e}'.format(e=typ),ls=ls) - ax.legend(loc='best') + ax.legend(loc='best',fontsize=14) ax.plot(wdn50,zdat,color=c2) ax.plot(wdn90,zdat,color=c1) ax.plot(wdn99,zdat,color=c3) @@ -1096,16 +1134,18 @@ def plot_quartiles(data,q1,q2,q3,z,ax,c1='goldenrod',c2='r',c3='k',split_updn=Fa ax.plot(wup50,zdat,color=c2,label='50th {e}'.format(e=typ),ls=ls) ax.plot(wup90,zdat,color=c1,label='90th {e}'.format(e=typ),ls=ls) ax.plot(wup99,zdat,color=c3,label='99th {e}'.format(e=typ),ls=ls) - ax.legend(loc='best') + ax.legend(loc='best',fontsize=14) + + ax.tick_params(axis='both',labelsize=14) + ax.set_ylabel('Height (km)',fontsize=16) - ax.set_ylabel('Height (km)') return ax def plot_verprof(data,z,ax,c='r',lab='',split_updn=False,ls = '-',typ='',thresh=-50): if split_updn == True: - pdat=data.load() + pdat=data.load().copy() pdat.values[pdat.values<-100] = np.nan wup = pdat.where(data>0) @@ -1121,10 +1161,10 @@ def plot_verprof(data,z,ax,c='r',lab='',split_updn=False,ls = '-',typ='',thresh= zdat = z.values ax.plot(wup50,zdat,color=c,label='{l} {e}'.format(l=lab,e=typ),ls=ls) ax.plot(wdn50,zdat,color=c,label='{l} {e}'.format(l=lab,e=typ),ls=ls) - ax.legend(loc='best') + ax.legend(loc='best',fontsize=14) else: - pdat =data.load() + pdat =data.load().copy() pdat.values[pdat.values= config['etime']) and (dates <= config['stime']): + radcdate=np.str(base[config['rdstart']:config['rdend']]) + #dates=datetime.datetime.strptime(radcdate,config['rdate_format']) + dates = datetime.datetime.strptime('{r}'.format(r=radcdate),config['rdate_format']) + sdt = datetime.datetime.strptime(config['sdatetime'],config['sdatetime_format']) + edt = datetime.datetime.strptime(config['edatetime'],config['edatetime_format']) + #if (dates >= config['etime']) and (dates <= config['stime']): + if (dates <= edt) and (dates >= sdt): #print cname #now find a sounding match mv = match_snd(dates,sdates) @@ -124,7 +136,7 @@ def reduce_dim(ds): try: t1= ds['time'][0].values except KeyError as ke: - print(f"{ke} skipping preprocessing") + #print(f"{ke} skipping preprocessing") return(ds) for v in ds.data_vars.keys(): try: @@ -146,10 +158,21 @@ def hasNumbers(inputString): def polarris_driver(configfile): - config = {} - print('ready to roll') + # ===== + # (1) Read in config file line by line. + # ===== + + config = {} # Load variable for config file data + #print('ready to roll') + print('\n#################################################') + print('############ Entering polarris_driver_new.py ####') + print('#################################################') + time.sleep(3) + + print('\nReading '+str(configfile[0])+'...') with open(configfile[0]) as f: - for line in f: + lines = [l for l in (line.strip() for line in f) if l] # NEW! Allow new lines in config file - can be skipped over! + for line in lines: #f: #print line if not line.startswith("#"): #print('line',line) @@ -164,7 +187,7 @@ def polarris_driver(configfile): #print vval,key if key.replace(" ", "") == 'image_dir': numck = True - if key.replace(" ", "") == 'radar_files': + if key.replace(" ", "") == 'rfiles': numck = True if numck is True or vval == 'None' or vval == 'True' or vval == 'False': @@ -175,17 +198,28 @@ def polarris_driver(configfile): config[(key.replace(" ", ""))] = vval else: config[(key.replace(" ", ""))] = vval - - print(config['radar_files']) + + time.sleep(3) + print('Read-in complete.\n') + + # ===== + # (2) Find input radar files and concatenate the data. Rename x, y, z variables. + # ===== + + print('Finding and concatenating radar files in '+config['rfiles']+'...') drop_vars=config['drop_vars'] - with open(config['radar_files'], 'r') as f: + with open(config['rfiles'], 'r') as f: + #print(config['rfiles']) rfiles = f.read().splitlines() - #rfiles= glob.glob('*.nc') - print((config['exper']),(config['mphys'])) + + #print((config['exper']),(config['mphys'])) + print('Station/experiment: '+config['exper']) + print('Input: '+config['mphys']) + time.sleep(3) if config['exper'] == 'MC3E' and config['mphys'] == 'obs': print("special handling for ",config['exper']) - file = open(config['radar_files'], "r") + file = open(config['rfiles'], "r") rf1=[] rf2=[] for line in file: @@ -205,59 +239,59 @@ def polarris_driver(configfile): # print('trying to read normally') # rvar = xr.open_mfdataset(rfiles,autoclose=True,concat_dim='d',preprocess=reduce_dim,combine='by_coords') # except ValueError as ve: - print("trying nesting") - rvar = xr.open_mfdataset(rfiles,autoclose=True,combine='nested',concat_dim='d',preprocess=reduce_dim) + #rfiles = glob.glob(config['rfiles']+"*") + #rvar = xr.open_mfdataset(rfiles,autoclose=True,combine='nested',concat_dim='d',preprocess=reduce_dim) + rvar = xr.open_mfdataset(rfiles,autoclose=True,concat_dim='d',preprocess=reduce_dim) + try: rvar = rvar.rename({'x0':'x'}) rvar = rvar.rename({'y0':'y'}) rvar = rvar.rename({'z0':'z'}) except: print('Dims do not need renaming') - print('Current dimensions:',rvar.dims) + #print('Current dimensions:',rvar.dims) if drop_vars == True: print("dropping extra variables for memory!") rvar= rvar.drop(['vrad03','vdop02','elev03','elev02','vdop03','vang02','vang03','vrad02','zhh02','zhh03','zdr02','zdr03','kdp02','kdp03','rhohv02','rhohv03']) - lon_0 = config['lon'] - lat_0 = config['lat'] - - lat_r = config['lat'] - lon_r = config['lon'] - - if config['snd_on'] == True: - smatch = find_snd_match(config) - #print("rfiles",rfiles[0]) - sfile = smatch[rfiles[0]] - print('matching sounding') - else: - smatch = None - + + print('Radar files ready.\n') + time.sleep(3) + # ===== + # (3) Get datetime objects from radar file names. + # ===== tm = [] for d in rfiles: - print(d) - dformat = config['wdate_format'] + #print(d) + dformat = config['rdate_format'] base = os.path.basename(d) - radcdate=np.str(base[config['time_parse'][0]:config['time_parse'][1]]) + radcdate=np.str(base[config['rdstart']:config['rdend']]) date=datetime.datetime.strptime(radcdate,dformat) tm.append(date) + # ===== + # (4) + # ===== + if config['dd_on']==True: - with open(config['dd_files'], 'r') as f: - dfiles1 = f.read().splitlines() + print('In your config file, dd_on is set to True.') + time.sleep(3) + with open(config['dfiles'], 'r') as g: + dfiles1 = g.read().splitlines() + #dfiles1 = glob.glob(config['dd_files']+"*") tmd = [] for d in dfiles1: dformat = config['ddate_format'] base = os.path.basename(d) # print('dd base',base,config['ddoff'],config['ddadd']) - radcdate = base[config['ddoff']:config['ddadd']] -# print (radcdate) + radcdate = base[config['ddstart']:config['ddend']] # print('dformat is',dformat,radcdate) if dformat == '%H%M': - hr=int(base[config['ddoff']:config['ddoff']+2]) - mn=int(base[config['ddoff']+2:config['ddoff']+4]) + hr=int(base[config['ddstart']:config['ddstart']+2]) + mn=int(base[config['ddstart']+2:config['ddstart']+4]) #print('hr','mn',hr,mn) # print radcdate #date=datetime.datetime.strptime(radcdate,dformat) @@ -268,7 +302,7 @@ def polarris_driver(configfile): dat2=datetime.datetime.strptime(radcdate,dformat) #dstart=datetime.datetime.strptime(config['date'],'%Y%m%d') tmd.append(dat2) - + print('Matching Dual-Doppler') dmatch = find_dd_match(rfiles,dfiles1,tm,tmd) #print('dmatch is ',dmatch) @@ -285,6 +319,13 @@ def polarris_driver(configfile): conv = np.zeros([rvar.dims['d'],rvar.dims['z'],rvar.dims['y'],rvar.dims['x']]) conv.fill(np.nan) + + # NEW! MultiDop only works if distance values are in metres, not km. Need a condition to convert back to km so that doppler and radar distances are comparable. + if np.array_equal(dvar.variables['x'].values, 1000.0*rvar.variables['x'].values): + dvar['x'] = rvar['x'] + dvar['y'] = rvar['y'] + dvar['z'] = rvar['z'] + xsubmin = np.where(rvar.variables['x']==np.min(dvar.variables['x']))[0][0] xsubmax = np.where(rvar.variables['x']==np.max(dvar.variables['x']))[0][0] @@ -292,45 +333,49 @@ def polarris_driver(configfile): ysubmax = np.where(rvar.variables['y']==np.max(dvar.variables['y']))[0][0] zsubmin = np.where(rvar.variables['z']==np.min(dvar.variables['z']))[0][0] - zsubmax = np.where(rvar.variables['z']==np.max(dvar.variables['z']))[0][0] - + zsubmax = np.where(rvar.variables['z']==np.max(dvar.variables['z']))[0][0] + for q,d in enumerate(dmatch.keys()): #print(q,'i outer',dmatch[d]) if dmatch[d] is not None: #print('good, dmatch is not none') dfile = dmatch[d] if dfile in dfiles1: + #print(dfile) i = dfiles1.index(dfile) #print(i,'i inner') - wvar[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['wname']].sel(d=i) - unew[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['uname']].sel(d=i) - vnew[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['vname']].sel(d=i) - conv[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['convname']].sel(d=i) - + #wvar[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['wname']].sel(d=i) + wvar[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['wname']][0,i,:,:,:] + unew[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['uname']][0,i,:,:,:] + vnew[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['vname']][0,i,:,:,:] + #conv[q,zsubmin:zsubmax+1,ysubmin:ysubmax+1,xsubmin:xsubmax+1] = dvar[config['convname']][i,:,:,:] rvar[config['wname']] = (['d','z','y','x'],wvar) rvar[config['uname']] = (['d','z','y','x'],unew) rvar[config['vname']] = (['d','z','y','x'],vnew) - rvar[config['convname']] = (['d','z','y','x'],conv) - - print('sending data to RadarData!') - rdata = RadarData.RadarData(rvar,tm,ddata = None,dz =config['dz_name'],zdr=config['dr_name'], - kdp=config['kd_name'],rho=config['rh_name'],temp=config['t_name'], - u=config['uname'],v=config['vname'],w=config['wname'],conv=config['convname'],x=config['xname'], - rr=config['rr_name'],band = config['band'],vr = config['vr_name'],lat_r=lat_r,lon_r=lon_r, - y=config['yname'],z=config['zname'],lat=config['latname'], lon=config['lonname'],lat_0=lat_0,lon_0=lon_0, - exper=config['exper'],mphys=config['mphys'],radar_name =config['radarname'], - z_thresh=0,conv_types = config['conv_types'], - strat_types = config['strat_types']) - + #rvar[config['convname']] = (['d','z','y','x'],conv) + + print('\nSending data to RadarData...') + + rdata = RadarData.RadarData(rvar,tm,ddata = None,dz=config['dz_name'],zdr=config['dr_name'],kdp=config['kd_name'],rho=config['rh_name'],temp=config['t_name'],u=config['uname'],v=config['vname'],w=config['wname'],conv=config['convname'],x=config['xname'],rr=config['rr_name'],band = config['band'],vr = config['vr_name'],lat_r=config['lat'],lon_r=config['lon'],y=config['yname'],z=config['zname'],lat=config['latname'], lon=config['lonname'],lat_0=config['lat'],lon_0=config['lon'],exper=config['exper'],mphys=config['mphys'],radar_name =config['radarname'],z_thresh=0,conv_types=config['conv_types'],strat_types=config['strat_types'],color_blind=config['cb_friendly']) + + if config['snd_on'] == True: + print('In your config file, snd_on is set to True.') + time.sleep(3) + smatch = find_snd_match(config) + #print("rfiles",rfiles[0]) + sfile = smatch[rfiles[0]] + print('Matching Sounding') + else: + smatch = None + if smatch is not None: - print ('Smatch',sfile) + print ('Found sounding match!',sfile,'\n') snd = SkewT.Sounding(sfile) rdata.add_sounding_object(snd) # this will add the sounding object to the radar object # and then will take the heights and temps rdata.interp_sounding() - if config['convert_Tk_Tc'] == True: print('converting T') rdata.convert_t() @@ -348,8 +393,8 @@ def polarris_driver(configfile): if config['comb_vicr'] == True: whvi = np.where(rdata.hid == 6) rdata.hid[whvi] = 3 - - + + #Do some quick masking of the data#### # mask = np.zeros([rdata.data.dims['d'],rdata.data.dims['z'],rdata.data.dims['y'],rdata.data.dims['x']]) # whbad = np.logical_or(np.logical_or(np.logical_or(np.logical_or(rdata.data[rdata.dz_name].values>-20.,rdata.data[rdata.zdr_name].values>-2.),rdata.data[rdata.kdp_name].values<10.),rdata.data[rdata.zdr_name].values<10.),rdata.data[rdata.dz_name].values<70.) @@ -367,5 +412,8 @@ def polarris_driver(configfile): # rdata.data[rdata.rho_name].values[whbad2] = np.nan # rdata.data[rdata.w_name].values[whbad2] = np.nan + print('\n\n#################################################') + print('####### Returning to run_ipolarris_new.py #######') + print('#################################################\n') return rdata, config diff --git a/run_ipolarris_new.py b/run_ipolarris_new.py index 40d8484..9372f26 100644 --- a/run_ipolarris_new.py +++ b/run_ipolarris_new.py @@ -1,40 +1,47 @@ -import numpy as np -import os +#=================================================== +#============== RUN_IPOLARRIS_NEW.PY =============== +#=================================================== + +# Import core Python packages +from collections import OrderedDict +import csv +import datetime import glob -from netCDF4 import Dataset import matplotlib matplotlib.use('Agg') +from netCDF4 import Dataset +import numpy as np +import os import matplotlib.pyplot as plt import pandas as pd +import sys +import warnings +warnings.filterwarnings('ignore') import xarray as xr -import numpy as np - +# Import iPOLARRIS functions +import GeneralFunctions as GF +from polarris_driver_new import polarris_driver +import plot_driver import RadarData -import datetime - import RadarConfig -import plot_driver -#from polarris_config import run_exper -#from polarris_config import get_data -import warnings -warnings.filterwarnings('ignore') -import GeneralFunctions as GF from skewPy import SkewT -from collections import OrderedDict -from polarris_driver_new import polarris_driver -import os -import sys +#--------------- Main Program ---------------- + +print('\n#############################################') +print('####### Starting run_ipolarris_new.py #######') +print('#############################################') -configfile = sys.argv[1:] +configfile = sys.argv[1:] # Feed config file name as arg #print sys.argv[1:] rdata, config = polarris_driver(configfile) #config['image_dir'] ='./' -print(config['extrax'],'EXTRA 1 is') -######################################### +#print(config['extrax'],'EXTRA 1 is') +# If a second argument is passed for WRF config file, produce a bunch of comparison plots! +# More comments in this section TBD! if sys.argv[2:]: configfile1 = sys.argv[2:] rdata2, config2 = polarris_driver(configfile1) @@ -75,7 +82,7 @@ fig,ax = plot_driver.plot_hid_comparison_cfad(rdata,rdata2,config=config,cscfad='convective',savefig=True) ################################################################################ -##################Now you can just start plotting!############################## +################## Now you can just start plotting! ############################ ################################################################################ ### To see the variables that are available to plot, type: @@ -86,57 +93,92 @@ ################################################################################ ##Plot a composite reflectivity at a given time. - - #tdate = datetime.datetime(2011,5,23,22,00) + # tdate = datetime.datetime(2011,5,23,22,00) # tdate = datetime.datetime(2006,1,23,18,0,0) # whdate = np.where(np.abs(tdate-np.array(rdata.date)) == np.min(np.abs(tdate-np.array(rdata.date)))) - print('In run_ipolarris...running the COMPOSITE figs.') - for i,d in enumerate(np.array(rdata.date)): - print('plotting composites by time....') - fig, ax = plot_driver.plot_composite(rdata,rdata.dz_name,i,cs_over=True) - print('made composite') - rtimematch = d - ax.set_title('{e} {r} composite {d:%Y%m%d %H%M}'.format(d=rtimematch,e=rdata.exper,r=rdata.radar_name)) - minlat = config['ylim'][0] - maxlat = config['ylim'][1] - minlon = config['xlim'][0] - maxlon = config['xlim'][1] - ax.set_extent([minlon, maxlon, minlat,maxlat]) - - plt.tight_layout() - plt.savefig('{i}Composite_{v}_{t:%Y%m%d%H%M}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],v=rdata.dz_name,t=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() - print('plotting cappis at 1 km by time...') - fig, ax = plt.subplots(1,1,figsize=(8,8)) - if 'd' in rdata.data[rdata.z_name].dims: - try: - whz = np.where(rdata.data[rdata.z_name].sel(d=i).values==config['z'])[0][0] - except IndexError as ie: - #print('checking z...',rdata.data[rdata.z_name].sel(d=i).values) - zdiffs = np.median(np.diff(rdata.data[rdata.z_name].values)) - whz = np.where(np.isclose(rdata.data[rdata.z_name].sel(d=i).values,config['z'],rtol=zdiffs))[0][0] - - else: - whz = np.where(rdata.data[rdata.z_name].values==config['z'])[0][0] - print('whz in run 122',whz) - rdata.cappi(rdata.dz_name,z=whz,ts=d,contour='CS',ax=ax) - ax.set_title('CAPPI DZ {t:%Y%m%d_%M%H%S} {h} km'.format(t=d,h=rdata.data['z'].values[whz])) - ax.set_xlim(config['xlim'][0],config['xlim'][1]) - ax.set_ylim(config['ylim'][0],config['ylim'][1]) -# ax.set_extent([minlon, maxlon, minlat,maxlat]) - plt.savefig('{i}DZ_CAPPI_{h}_{v}_{t:%Y%m%d%H%M}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],h=config['z'],v=rdata.dz_name,t=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() + if config['compo_ref']: + + print('IN RUN_IPOLARRIS_NEW... creating COMPOSITE figures.') + print('\nPlotting composites by time for variable '+rdata.dz_name+'...') + + for i,rtimematch in enumerate(np.array(rdata.date)): + + fig, ax = plot_driver.plot_composite(rdata,rdata.dz_name,i,cs_over=False,statpt=True) + #ax.set_title('{e} {r} Composite {d:%Y-%m-%d %H%M} UTC'.format(d=rtimematch,e=rdata.exper,r=rdata.radar_name), fontsize=18) + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + ax.text(1, 1, '{d:%Y-%m-%d %H:%M:%S} UTC'.format(d=rtimematch), horizontalalignment='right', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + + #minlat = config['ylim'][0] + #maxlat = config['ylim'][1] + #minlon = config['xlim'][0] + #maxlon = config['xlim'][1] + #ax.set_extent([minlon, maxlon, minlat,maxlat]) + + #plt.tight_layout() + os.makedirs(config['image_dir']+'composite_'+rdata.dz_name+'/',exist_ok=True) + plt.savefig('{i}composite_{v}_{d:%Y-%m-%d_%H%M%S}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir']+'composite_'+rdata.dz_name+'/',v=rdata.dz_name,d=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400, bbox_inches='tight') + plt.close() + print(rtimematch) + + print('\nDone! Saved to '+config['image_dir']+'composite_'+rdata.dz_name+'/') + print('Moving on.\n') + + if config['cappi_ref']: + + print('\nIN RUN_IPOLARRIS_NEW... creating CAPPI figures.') + print('\nPlotting CAPPIs at height z = '+str(config['z'])+'km by time for variable '+rdata.dz_name+'...') + + for i,rtimematch in enumerate(np.array(rdata.date)): + + fig, ax = plot_driver.plot_cappi(rdata,rdata.dz_name,rdata.z_name,config['z'],i,rtimematch,cs_over=False,statpt=True) + + #ax.set_xlim(config['xlim'][0],config['xlim'][1]) + #ax.set_ylim(config['ylim'][0],config['ylim'][1]) + #ax.set_extent([minlon, maxlon, minlat,maxlat]) + + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + ax.text(1, 1, '{d:%Y-%m-%d %H:%M:%S} UTC'.format(d=rtimematch), horizontalalignment='right', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + ax.text(0.99, 0.99, 'z = {a} km'.format(a=config['z']), horizontalalignment='right',verticalalignment='top', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) + + os.makedirs(config['image_dir']+'cappi_'+rdata.dz_name+'/',exist_ok=True) + plt.savefig('{i}dz_cappi_{h}_{v}_{t:%Y-%m-%d_%H%M%S}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir']+'cappi_'+rdata.dz_name+'/',h=config['z'],v=rdata.dz_name,t=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') + plt.close() + + print(rtimematch) + + print('\nDone! Saved to '+config['image_dir']+'cappi_'+rdata.dz_name+'/') + print('Moving on.\n') + + if config['cappi_rr']: - fig, ax = plt.subplots(1,1,figsize=(8,8)) - # whz = np.where(rdata.data[rdata.z_name].values==config['z'])[0][0] - rdata.cappi(rdata.rr_name,z=whz,ts=d,contour='CS',ax=ax) - ax.set_xlim(config['xlim'][0],config['xlim'][1]) - ax.set_ylim(config['ylim'][0],config['ylim'][1]) - ax.set_title('CAPPI RR {t:%Y%m%d_%M%D%S} {h} km'.format(t=d,h=rdata.data['z'].values[2])) - plt.savefig('{i}RR_CAPPI_{h}_{v}_{t:%Y%m%d%H%M}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],h=config['z'],v=rdata.dz_name,t=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() - + print('\nIN RUN_IPOLARRIS_NEW... creating CAPPI figures.') + print('\nPlotting CAPPIs at height z = '+str(config['z'])+'km by time for variable '+rdata.dz_name+'...') + + print('Plotting CAPPIs at height z = '+str(config['z'])+'km by time for variable '+rdata.rr_name+'...') + + for i,rtimematch in enumerate(np.array(rdata.date)): + + fig, ax = plot_driver.plot_cappi(rdata,rdata.rr_name,rdata.z_name,config['z'],i,rtimematch,cs_over=False,statpt=True) + #fig, ax = plt.subplots(1,1,figsize=(10,8)) + # whz = np.where(rdata.data[rdata.z_name].values==config['z'])[0][0] + #rdata.cappi(rdata.rr_name,z=whz,ts=rtimematch,contour='CS',ax=ax) + #ax.set_xlim(config['xlim'][0],config['xlim'][1]) + #ax.set_ylim(config['ylim'][0],config['ylim'][1]) + #ax.set_title('CAPPI RR {t:%Y%m%d_%M%D%S} {h} km'.format(t=rtimematch,h=rdata.data['z'].values[2])) + + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + ax.text(1, 1, '{d:%Y-%m-%d %H:%M:%S} UTC'.format(d=rtimematch), horizontalalignment='right', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + ax.text(0.99, 0.99, 'z = {a} km'.format(a=config['z']), horizontalalignment='right',verticalalignment='top', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) + + os.makedirs(config['image_dir']+'cappi_'+rdata.rr_name+'/',exist_ok=True) + plt.savefig('{i}rr_cappi_{h}_{v}_{t:%Y-%m-%d_%H%M%S}_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir']+'cappi_'+rdata.rr_name+'/',h=config['z'],v=rdata.dz_name,t=rtimematch,e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') + plt.close() + + print(rtimematch) + + print('\nDone! Saved to '+config['image_dir']+'cappi_'+rdata.rr_name+'/') + print('Moving on.\n') # tdate = datetime.datetime(2006,1,23,18,00) # whdate = np.where(np.abs(tdate-np.array(rdata.date)) == np.min(np.abs(tdate-np.array(rdata.date)))) @@ -148,165 +190,258 @@ # plt.clf() ################################################################################ - ##Calculate a timeseries for writing out - rrstratu,rrconvu,rrallu = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=-0.1) - rrstrat,rrconv,rrall = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=0.0) - - import csv - tformat = '%Y%m%d-%H%M%S' - with open('{i}{e}_rr_uncondmean_stats.txt'.format(i=config['image_dir'],e=config['exper']), mode='w') as csv_file: - v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) - v_writer.writerow(['Date', 'Unc_Conv_RR', 'Unc_Strat_RR', 'Unc_Tot_RR']) - for i,v in enumerate(rdata.date): - print( v) - tim = v.strftime(tformat) - dum =[tim,rrconvu[i].values,rrstratu[i].values,rrallu[i].values] - v_writer.writerow(dum) - - tformat = '%Y%m%d-%H%M%S' - with open('{i}{e}_rr_condmean_stats.txt'.format(i=config['image_dir'],e=config['exper']), mode='w') as csv_file: - v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) - v_writer.writerow(['Date', 'Conv_RR', 'Strat_RR', 'Tot_RR']) - for i,v in enumerate(rdata.date): - print (v) - tim = v.strftime(tformat) - dum =[tim,rrconv[i].values,rrstrat[i].values,rrall[i].values] - v_writer.writerow(dum) - - conv = np.where(rdata.data[rdata.cs_name].values == 2) - strat = np.where(rdata.data[rdata.cs_name].values == 1) - hist, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values)),bins=np.logspace(-1,2.4,40)) - histc, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values[conv])),bins=np.logspace(-1,2.4,40)) - hists, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values[strat])),bins=np.logspace(-1,2.4,40)) - - - tformat = '%Y%m%d-%H%M%S' - with open('{i}{e}_rr_histgram_{m}.txt'.format(i=config['image_dir'],e=config['exper'],m=config['mphys']), mode='w') as csv_file: - v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) - v_writer.writerow(['Date', 'Con', 'Strat', 'Tot']) - for i,v in enumerate(eg[:-1]): - dum =[v,histc[i],hists[i],hist[i]] - v_writer.writerow(dum) - - - ###Areas - - rrstratu_area,rrconvu_area,rrallu_area = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=-0.1,areas=True) - rrstrat_area,rrconv_area,rrall_area = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=0.0,areas=True) - - #grid_area=rdata.radar_area() - rain_area = rdata.radar_area - import csv - tformat = '%Y%m%d-%H%M%S' - with open('{i}{e}_domain_area_stats.txt'.format(i=config['image_dir'],e=config['exper']), mode='w') as csv_file: - v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) - v_writer.writerow(['Date', 'Unc_Con', 'Unc_Strat', 'Unc_Tot']) - for i,v in enumerate(rdata.date): - print (v) - tim = v.strftime(tformat) - dum =[tim,rrconvu_area[i].values.astype(float)*rdata.dx*rdata.dy,rrstratu_area[i].values.astype(float)*rdata.dx*rdata.dy,rrallu_area[i].values.astype(float)*rdata.dx*rdata.dy] - v_writer.writerow(dum) - - tformat = '%Y%m%d-%H%M%S' - with open('{i}{e}_rel_frequency_stats.txt'.format(i=config['image_dir'],e=config['exper']), mode='w') as csv_file: - v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) - v_writer.writerow(['Date', 'Conv', 'Strat', 'Tot']) - for i,v in enumerate(rdata.date): - print( v) - tim = v.strftime(tformat) - dum =[tim,rrconv[i].values*rdata.dx*rdata.dy/rain_area*100.,rrstrat[i].values*rdata.dx*rdata.dy/rain_area*100.,rrall[i].values*rdata.dx*rdata.dy/rain_area*100.] - v_writer.writerow(dum) - - - - - ################################################################################ - ##First make a timeseries of rain rate, unconditional and conditional. This puts strat, conv, and total on the same plot but you can split the out by putting cs==False. - ## The conditional rain rate is achieved by sending threshold = 0. - fig,ax = plt.subplots(1,1,figsize=(10,10)) - ax = plot_driver.plot_timeseries(rdata.data[rdata.rr_name],rdata.date,ax,cs=True,rdata=rdata,thresh=0,zlev=1,make_zeros=False) - ax = plot_driver.plot_timeseries(rdata.data[rdata.rr_name],rdata.date,ax,cs=True,rdata=rdata,thresh=0,zlev=1,ls='--',typ='uncond',make_zeros=True)#,zlev=0) - - ax.set_ylabel('Rain Rate (mm/hr)') - ax.set_title('Precipitation Timeseries ') - plt.tight_layout() - plt.savefig('{i}Precip_timeseries_convstrat_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() - + + if config['rrstats_txt']: + + print('\nIN RUN_IPOLARRIS_NEW... creating text files.') + print('Printing unconditional-mean statistics for variable '+rdata.rr_name+'...') + + ##Calculate a timeseries for writing out + rrstratu,rrconvu,rrallu = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2.5,cs_flag=True,thresh=-0.1) + rrstrat,rrconv,rrall = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2.5,cs_flag=True,thresh=0.0) + + tformat = '%Y%m%d-%H%M%S' + os.makedirs(config['image_dir']+'txtfiles/',exist_ok=True) + with open('{i}{e}_rr_uncondmean_stats.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Unc_Conv_RR', 'Unc_Strat_RR', 'Unc_Tot_RR']) + for i,v in enumerate(rdata.date): + #print( v) + tim = v.strftime(tformat) + dum =[tim,rrconvu[i].values,rrstratu[i].values,rrallu[i].values] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Printing conditional-mean statistics for variable '+rdata.rr_name+'...') + + with open('{i}{e}_rr_condmean_stats.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Conv_RR', 'Strat_RR', 'Tot_RR']) + for i,v in enumerate(rdata.date): + #print (v) + tim = v.strftime(tformat) + dum =[tim,rrconv[i].values,rrstrat[i].values,rrall[i].values] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Printing relative frequency statistics for variable '+rdata.rr_name+'...') + + rain_area = rdata.radar_area + with open('{i}{e}_rel_frequency_stats.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Conv', 'Strat', 'Tot']) + for i,v in enumerate(rdata.date): + #print( v) + tim = v.strftime(tformat) + dum =[tim,rrconv[i].values*rdata.dx*rdata.dy/rain_area*100.,rrstrat[i].values*rdata.dx*rdata.dy/rain_area*100.,rrall[i].values*rdata.dx*rdata.dy/rain_area*100.] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Moving on.\n') + + if config['rrhist_txt']: + + print('\nIN RUN_IPOLARRIS_NEW... creating text files.') + print('Printing histogram data for variable '+rdata.rr_name+'...') + + conv = np.where(rdata.data[rdata.cs_name].values == 2) + strat = np.where(rdata.data[rdata.cs_name].values == 1) + hist, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values)),bins=np.logspace(-1,2.4,40)) + histc, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values[conv])),bins=np.logspace(-1,2.4,40)) + hists, eg = np.histogram(np.ravel((rdata.data[rdata.rr_name].values[strat])),bins=np.logspace(-1,2.4,40)) + + #tformat = '%Y%m%d-%H%M%S' + os.makedirs(config['image_dir']+'txtfiles/',exist_ok=True) + with open('{i}{e}_rr_histgram_{m}.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper'],m=config['mphys']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Con', 'Strat', 'Tot']) + for i,v in enumerate(eg[:-1]): + dum =[v,histc[i],hists[i],hist[i]] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Moving on.\n') + + if config['rrstats_areas_txt']: + ###Areas + print('\nIN RUN_IPOLARRIS_NEW... creating text files.') + print('Printing domain area statistics for '+rdata.rr_name+'...') + + rrstratu_area,rrconvu_area,rrallu_area = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=-0.1,areas=True) + rrstrat_area,rrconv_area,rrall_area = rdata.calc_timeseries_stats(rdata.rr_name,ht_lev=2,cs_flag=True,thresh=0.0,areas=True) + + #grid_area=rdata.radar_area() + rain_area = rdata.radar_area + tformat = '%Y%m%d-%H%M%S' + os.makedirs(config['image_dir']+'txtfiles/',exist_ok=True) + with open('{i}{e}_domain_area_stats.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Unc_Con', 'Unc_Strat', 'Unc_Tot']) + for i,v in enumerate(rdata.date): + print (v) + tim = v.strftime(tformat) + dum =[tim,rrconvu_area[i].values.astype(float)*rdata.dx*rdata.dy,rrstratu_area[i].values.astype(float)*rdata.dx*rdata.dy,rrallu_area[i].values.astype(float)*rdata.dx*rdata.dy] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Moving on.\n') + + ''' + if config['rrstats_txt']: + + tformat = '%Y%m%d-%H%M%S' + os.makedirs(config['image_dir']+'txtfiles/',exist_ok=True) + with open('{i}{e}_rel_frequency_stats.txt'.format(i=config['image_dir']+'txtfiles/',e=config['exper']), mode='w') as csv_file: + v_writer = csv.writer(csv_file, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_NONNUMERIC) + v_writer.writerow(['Date', 'Conv', 'Strat', 'Tot']) + for i,v in enumerate(rdata.date): + print( v) + tim = v.strftime(tformat) + dum =[tim,rrconv[i].values*rdata.dx*rdata.dy/rain_area*100.,rrstrat[i].values*rdata.dx*rdata.dy/rain_area*100.,rrall[i].values*rdata.dx*rdata.dy/rain_area*100.] + v_writer.writerow(dum) + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Moving on.\n') + ''' + if config['rr_timeseries']: + + print('\nIN RUN_IPOLARRIS_NEW... creating timeseries.') + print('Plotting timeseries for variable '+rdata.rr_name+'...') + + ################################################################################ + ##First make a timeseries of rain rate, unconditional and conditional. This puts strat, conv, and total on the same plot but you can split the out by putting cs==False. + ## The conditional rain rate is achieved by sending threshold = 0. + fig,ax = plt.subplots(1,1,figsize=(12,8)) + ax = plot_driver.plot_timeseries(rdata.data[rdata.rr_name],rdata.date,ax,cs=True,rdata=rdata,thresh=0,zlev=1,make_zeros=False) + ax = plot_driver.plot_timeseries(rdata.data[rdata.rr_name],rdata.date,ax,cs=True,rdata=rdata,thresh=0,zlev=1,ls='--',typ='uncond',make_zeros=True)#,zlev=0) + + ax.set_ylabel('Rain Rate (mm/hr)',fontsize=16) + #ax.set_title('Precipitation Timeseries ') + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + #plt.tight_layout() + plt.savefig('{i}precip_timeseries_convstrat_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') + plt.close() + + print('\nDone! Saved to '+config['image_dir']) + print('Moving on.\n') + ############################################################################ ################################################################################ ##Next let's make quantile (50,90,99) plots of the vertical velocity. This splits it by up and down, but you can turn split_updn == False if rdata.w_name is not None: - fig,ax = plt.subplots(1,1,figsize=(10,10)) - ax = plot_driver.plot_quartiles(rdata.data[rdata.w_name],0.9,0.5,0.99,rdata.data[rdata.z_name],ax,split_updn=True) - ax = plot_driver.plot_quartiles(rdata.data[rdata.w_name],0.9,0.5,0.99,rdata.data[rdata.z_name],ax,split_updn=False) - ax.set_xlabel('Vertical velocity m/s') - ax.set_title('Vertical velocity profiles') - plt.tight_layout() - plt.savefig('{i}Quantile_vvel_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() - - p99u,p90u,p50u,ht = rdata.percentile(wup=True) - p99d,p90d,p50d,ht = rdata.percentile(wdown=True) - p99a,p90a,p50a,ht = rdata.percentile(wdown=False) - file = open('{i}{e}_{m}_updown_percentiles.txt'.format(i=config['image_dir'],e=rdata.exper,m=rdata.mphys),'w') + if config['vv_profiles']: - file.write("Updraft\n") - file.write("Height (km). P99. P90. P50\n") - for i,h in enumerate(ht): - file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99u[i],p2=p90u[i],p3=p50u[i])) - - file.write("Downdraft\n") - file.write("Height (km). P99. P90. P50\n") - for i,h in enumerate(ht): - file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99d[i],p2=p90d[i],p3=p50d[i])) - - file.write("ALL\n") - file.write("Height (km). P99. P90. P50\n") - for i,h in enumerate(ht): - file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99a[i],p2=p90a[i],p3=p50a[i])) - + print('\nIN RUN_IPOLARRIS_NEW... creating vertical profile figure.') + print('Plotting vertical profile for variable '+rdata.w_name+'...') + + fig,ax = plt.subplots(1,1,figsize=(12,8)) + ax = plot_driver.plot_quartiles(rdata.data[rdata.w_name],0.9,0.5,0.99,rdata.data[rdata.z_name],ax,split_updn=True) + ax = plot_driver.plot_quartiles(rdata.data[rdata.w_name],0.9,0.5,0.99,rdata.data[rdata.z_name],ax,split_updn=False) + ax.set_xlabel('Vertical Velocity (m/s)',fontsize=16) + #ax.set_title('Vertical velocity profiles') + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + #plt.tight_layout() + plt.savefig('{i}quantile_vvel_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') + plt.close() + + print('\nDone! Saved to '+config['image_dir']) + print('Moving on.\n') + + if config['percentiles_txt']: - file.close() + print('\nIN RUN_IPOLARRIS_NEW... creating percentile text file.') + print('Printing percentile data for variable '+rdata.w_name+'...') + + p99u,p90u,p50u,ht = rdata.percentile(wup=True) + p99d,p90d,p50d,ht = rdata.percentile(wdown=True) + p99a,p90a,p50a,ht = rdata.percentile(wdown=False) + + os.makedirs(config['image_dir']+'txtfiles/',exist_ok=True) + file = open('{i}{e}_{m}_updown_percentiles.txt'.format(i=config['image_dir']+'txtfiles/',e=rdata.exper,m=rdata.mphys),'w') + + file.write("Updraft\n") + file.write("Height (km). P99. P90. P50\n") + for i,h in enumerate(ht): + file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99u[i],p2=p90u[i],p3=p50u[i])) + + file.write("Downdraft\n") + file.write("Height (km). P99. P90. P50\n") + for i,h in enumerate(ht): + file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99d[i],p2=p90d[i],p3=p50d[i])) + + file.write("ALL\n") + file.write("Height (km). P99. P90. P50\n") + for i,h in enumerate(ht): + file.write("{h} {p1} {p2} {p3}\n".format(h=h,p1=p99a[i],p2=p90a[i],p3=p50a[i])) + + file.close() + + print('\nDone! Saved to '+config['image_dir']+'txtfiles/') + print('Moving on.\n') + else: - print("No vertical velocity data.") + print("\nNo vertical velocity data.") + print('Moving on.\n') + ################################################################################ ################################################################################ - ##Next let's make mean vertical profile of reflectivity - fig,ax = plt.subplots(1,1,figsize=(10,10)) - ax = plot_driver.plot_verprof(rdata.data[rdata.dz_name],rdata.data[rdata.z_name],ax,split_updn=False,lab='dz',thresh=-50) - ax.set_title('Vertical profile of reflectivity') - ax.set_xlabel('Reflectivity') - plt.tight_layout() - plt.savefig('{i}MeanProfile_refl_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) + + if config['vert_ref']: + + print('\nIN RUN_IPOLARRIS_NEW... creating vertical profile figure.') + print('Plotting vertical profile for variable '+rdata.dz_name+'...') + + ##Next let's make mean vertical profile of reflectivity + fig,ax = plt.subplots(1,1,figsize=(12,8)) + ax = plot_driver.plot_verprof(rdata.data[rdata.dz_name],rdata.data[rdata.z_name],ax,split_updn=False,lab='dz',thresh=-50) + #ax.set_title('Vertical profile of reflectivity') + ax.set_xlabel('Reflectivity (dBZ)',fontsize=16) + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + #plt.tight_layout() + plt.savefig('{i}meanprofile_refl_{e}_{m}_{x}.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') - plt.close() + plt.close() + + print('\nDone! Saved to '+config['image_dir']) + print('Moving on.\n') + + if config['refcfad']: + + print('\nIN RUN_IPOLARRIS_NEW... creating CFAD figure.') + print('Plotting CFAD for variable '+rdata.dz_name+'...') + ################################################################################ ##Next let's make a reflectivity CFAD # cfaddat,vbins = plot_driver.cfad(rdata.data[rdata.dz_name],rdata,rdata.data[rdata.z_name],var=rdata.dz_name,nbins=40) - cfaddat,vbins,r1ht = rdata.cfad(rdata.dz_name,ret_z=1,z_resolution=1.0,value_bins=np.arange(0,82,2),cscfad=False) + cfaddat,vbins,r1ht = rdata.cfad(rdata.dz_name,ret_z=1,z_resolution=1.0,value_bins=np.arange(0,82,2),cscfad=False) - fig,ax = plt.subplots(1,1,figsize=(10,10)) - ax = plot_driver.plot_cfad(cfaddat, hts = r1ht, vbins = vbins,ax=ax,cfad_on = 0,tspan = config['date'],maxval=20,cont=True,levels = True) + fig,ax = plt.subplots(1,1,figsize=(12,8)) + ax = plot_driver.plot_cfad(fig,cfaddat, hts = r1ht, vbins = vbins,ax=ax,cfad_on = 0,tspan = config['date'],maxval=20,cont=True,levels = True) + + ax.set_xlabel('Reflectivity (dBZ)',fontsize=16) + #ax.set_title('{c} CFAD'.format(c=rdata.exper)) + ax.text(0, 1, '{e} {r}'.format(e=rdata.exper,r=rdata.radar_name), horizontalalignment='left', verticalalignment='bottom', size=16, color='k', zorder=10, weight='bold', transform=ax.transAxes) # (a) Top-left + #plt.tight_layout() + plt.savefig('{i}CFAD_refl_{e}_{m}_{x}_new.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400,bbox_inches='tight') + plt.close() + + print('\nDone! Saved to '+config['image_dir']) + print('Moving on.\n') + + ################################################################################# - ax.set_xlabel('Reflectivity') - ax.set_ylabel('Height (km)') - ax.set_title('{c} CFAD'.format(c=rdata.exper)) - plt.tight_layout() - plt.savefig('{i}CFAD_refl_{e}_{m}_{x}_new.{p}'.format(p=config['ptype'],i=config['image_dir'],e=rdata.exper,m=rdata.mphys,x=config['extrax']),dpi=400) - plt.close() - - flags = {} for k in eval(config['ks']): flags[k]=config[k] if any(flags.values()) == True: - plot_driver.make_single_pplots(rdata,flags,config) + print('\n############################################') + print('####### Exiting run_ipolarris_new.py #######') + print('############################################\n') - + plot_driver.make_single_pplots(rdata,flags,config) diff --git a/samples/KLGX_20151203_010200_V06.nc b/samples/KLGX_20151203_010200_V06.nc new file mode 100644 index 0000000..bf16cfc Binary files /dev/null and b/samples/KLGX_20151203_010200_V06.nc differ diff --git a/samples/KLGX_20151203_010607_V06.nc b/samples/KLGX_20151203_010607_V06.nc new file mode 100644 index 0000000..ea2b38a Binary files /dev/null and b/samples/KLGX_20151203_010607_V06.nc differ diff --git a/samples/KLGX_KATX_20151203_010200_V06_dopvel_gridded.nc b/samples/KLGX_KATX_20151203_010200_V06_dopvel_gridded.nc new file mode 100644 index 0000000..64dd694 Binary files /dev/null and b/samples/KLGX_KATX_20151203_010200_V06_dopvel_gridded.nc differ diff --git a/samples/KLGX_KATX_20151203_010607_V06_dopvel_gridded.nc b/samples/KLGX_KATX_20151203_010607_V06_dopvel_gridded.nc new file mode 100644 index 0000000..160d1ce Binary files /dev/null and b/samples/KLGX_KATX_20151203_010607_V06_dopvel_gridded.nc differ diff --git a/samples/UIL_20151203_000000.txt b/samples/UIL_20151203_000000.txt new file mode 100644 index 0000000..749620e --- /dev/null +++ b/samples/UIL_20151203_000000.txt @@ -0,0 +1,148 @@ + pressure height temperature dewpoint direction speed u_wind v_wind station station_number time latitude longitude elevation pw +0 1004.0 62 12.2 9.4 130.0 11.0 -8.42648887430876 7.070663706551931 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +1 1000.0 88 11.8 9.9 135.0 14.0 -9.899494936611665 9.899494936611664 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +2 974.3 305 10.3 9.1 145.0 28.0 -16.060140217829286 22.936257240091773 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +3 939.2 610 8.1 8.1 160.0 37.0 -12.654745303049749 34.76862696907861 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +4 937.0 630 8.0 8.0 161.0 37.0 -12.046021714914794 34.98418729717472 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +5 925.0 736 7.6 7.6 165.0 39.0 -10.09394275899832 37.67110722527366 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +6 905.1 914 6.7 6.7 170.0 40.0 -6.94592710667721 39.392310120488325 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +7 872.1 1219 5.2 5.2 180.0 36.0 -4.408728476930472e-15 36.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +8 860.0 1334 4.6 4.6 180.0 35.0 -4.2862637970157365e-15 35.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +9 850.0 1429 3.2 2.5 180.0 34.0 -4.1637991171010006e-15 34.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +10 845.0 1477 2.8 1.3 180.0 32.0 -3.91886975727153e-15 32.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +11 840.1 1524 2.8 0.8 180.0 31.0 -3.796405077356795e-15 31.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +12 834.0 1583 2.8 0.2 181.0 32.0 0.5584770059930764 31.99512624500452 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +13 828.0 1642 2.2 -2.8 182.0 32.0 1.1167838944800286 31.980506464611064 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +14 809.0 1829 1.2 -5.4 185.0 34.0 2.96329525342037 33.87061973511935 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +15 805.0 1869 1.0 -6.0 184.0 34.0 2.371720107300259 33.917177708834025 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +16 793.0 1989 0.6 -6.4 182.0 35.0 1.2214823845875313 34.97867894566835 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +17 778.8 2134 -0.2 -4.3 180.0 35.0 -4.2862637970157365e-15 35.0 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +18 769.0 2236 -0.7 -2.8 182.0 35.0 1.2214823845875313 34.97867894566835 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +19 756.0 2372 -1.5 -9.5 184.0 35.0 2.441476581044385 34.91474175909385 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +20 752.0 2414 -1.1 -7.1 185.0 35.0 3.050450996168028 34.8668144332111 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +21 749.8 2438 -1.2 -7.3 185.0 35.0 3.050450996168028 34.8668144332111 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +22 721.5 2743 -2.7 -9.5 200.0 45.0 15.39090644965509 42.28616793536588 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +23 714.0 2827 -3.1 -10.1 202.0 45.0 16.85729670371604 41.723273455505435 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +24 700.0 2983 -4.3 -11.3 205.0 46.0 19.44044004007217 41.690158203685904 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +25 696.0 3028 -4.7 -11.7 205.0 46.0 19.44044004007217 41.690158203685904 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +26 694.2 3048 -4.8 -13.0 205.0 46.0 19.44044004007217 41.690158203685904 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +27 684.0 3165 -5.5 -20.5 204.0 46.0 18.70988558148681 42.02309105155964 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +28 677.0 3246 -5.9 -21.9 203.0 47.0 18.364363038995858 43.2637281122647 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +29 660.0 3444 -7.9 -17.9 202.0 47.0 17.606509890547866 43.57764116463901 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +30 647.0 3599 -9.3 -19.3 200.0 48.0 16.416966879632096 45.105245797723605 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +31 642.0 3658 -9.9 -15.3 200.0 48.0 16.416966879632096 45.105245797723605 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +32 640.0 3683 -10.1 -13.6 200.0 48.0 16.416966879632096 45.105245797723605 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +33 634.0 3755 -10.5 -13.3 198.0 46.0 14.214781741247576 43.74859974957707 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +34 620.0 3927 -11.1 -13.9 196.0 44.0 12.128043655947955 42.29551462128603 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +35 617.2 3962 -10.7 -13.1 195.0 43.0 11.129218939408394 41.53481053042994 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +36 611.0 4040 -9.7 -11.4 200.0 46.0 15.732926592980759 43.22586055615179 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +37 593.2 4267 -11.1 -12.2 215.0 55.0 31.54670399930754 45.05336243589454 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +38 571.0 4560 -12.9 -13.2 225.0 53.0 37.476659402887016 37.47665940288703 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +39 559.0 4722 -13.7 -14.5 230.0 51.0 39.068266599067876 32.78216809401352 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +40 547.7 4877 -14.9 -15.8 235.0 50.0 40.95760221444958 28.67882181755232 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +41 544.0 4929 -15.3 -16.3 237.0 50.0 41.9335283972712 27.23195175075135 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +42 537.0 5027 -16.3 -18.4 240.0 49.0 42.43524478543748 24.50000000000002 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +43 532.0 5097 -17.3 -21.7 242.0 48.0 42.381484457228495 22.534635013722756 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +44 526.0 5182 -17.8 -22.2 245.0 47.0 42.59646599072256 19.86305830181286 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +45 502.0 5530 -19.7 -24.5 240.0 48.0 41.56921938165304 24.00000000000002 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +46 500.0 5560 -19.9 -24.5 240.0 48.0 41.56921938165304 24.00000000000002 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +47 488.0 5739 -21.1 -23.2 237.0 46.0 38.578846125489505 25.05339561069124 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +48 478.0 5891 -22.3 -23.9 234.0 44.0 35.596747752497684 25.862551100868824 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +49 464.7 6096 -23.9 -26.0 230.0 42.0 32.173866610997074 26.997079606834664 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +50 445.6 6401 -26.2 -29.2 230.0 42.0 32.173866610997074 26.997079606834664 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +51 438.0 6525 -27.1 -30.5 229.0 44.0 33.20722152980196 28.866597275582325 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +52 400.0 7170 -32.7 -36.5 225.0 54.0 38.18376618407356 38.183766184073576 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +53 399.0 7188 -32.7 -36.4 225.0 54.0 38.18376618407356 38.183766184073576 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +54 395.0 7259 -33.3 -38.0 225.0 55.0 38.89087296526011 38.89087296526012 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +55 375.0 7620 -36.5 -40.8 225.0 59.0 41.7193000900063 41.719300090006314 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +56 361.0 7886 -38.9 -42.9 227.0 60.0 43.88122209715023 40.91990160374991 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +57 349.0 8118 -40.7 -45.7 229.0 61.0 46.03728439358908 40.01960076842095 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +58 313.4 8839 -47.4 -51.1 235.0 63.0 51.60657879020647 36.13531549011592 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +59 302.0 9086 -49.7 -52.9 235.0 75.0 61.43640332167437 43.01823272632848 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +60 300.0 9130 -49.9 -52.9 235.0 77.0 63.07470741025235 44.16538559903057 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +61 299.3 9144 -50.0 -53.0 235.0 77.0 63.07470741025235 44.16538559903057 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +62 250.0 10290 -59.3 -63.8 245.0 104.0 94.25600985181161 43.95229922103271 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +63 247.0 10366 -59.9 -64.3 246.0 107.0 97.7493639677583 43.52082080911061 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +64 223.8 10973 -63.9 -68.7 250.0 128.0 120.28065546059628 43.77857834568557 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +65 223.0 10994 -64.1 -68.8 250.0 128.0 120.28065546059628 43.77857834568557 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +66 206.0 11481 -67.3 -72.3 255.0 107.0 103.35406341293032 27.693637825969706 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +67 202.0 11600 -66.5 -71.5 255.0 102.0 98.52443428148497 26.399542600457103 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +68 200.0 11660 -66.7 -71.5 255.0 99.0 95.62665680261776 25.623085465149543 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +69 196.0 11782 -66.5 -71.3 255.0 98.0 94.6607309763287 25.36426642004702 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +70 192.7 11887 -65.2 -72.5 255.0 98.0 94.6607309763287 25.36426642004702 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +71 179.0 12342 -59.5 -77.5 250.0 69.0 64.83879083422768 23.599389889471126 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +72 171.0 12629 -58.1 -81.1 247.0 51.0 46.945747526074456 19.927287552952965 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +73 166.3 12802 -58.8 -82.1 245.0 40.0 36.25231148146601 16.904730469627964 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +74 158.5 13106 -60.1 -83.9 225.0 54.0 38.18376618407356 38.183766184073576 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +75 155.0 13245 -60.7 -84.7 227.0 62.0 45.34392950038857 42.283898323874915 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +76 150.9 13411 -59.6 -84.4 230.0 71.0 54.38915546144743 45.63792028774431 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +77 150.0 13450 -59.3 -84.3 230.0 71.0 54.38915546144743 45.63792028774431 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +78 140.0 13885 -56.3 -84.3 237.0 55.0 46.12688123699832 29.955146925826483 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +79 130.5 14326 -58.1 -85.6 245.0 38.0 34.439695907392704 16.059493946146567 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +80 123.0 14700 -59.7 -86.7 239.0 50.0 42.85836503510561 25.751903745502723 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +81 120.0 14855 -57.5 -85.5 236.0 55.0 45.597066490527304 30.755609690891063 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +82 118.5 14935 -57.1 -85.4 235.0 58.0 47.510818568761515 33.26743330836069 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +83 115.0 15124 -56.1 -85.1 244.0 53.0 47.636084453855844 23.23367077982112 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +84 113.0 15236 -56.5 -85.5 250.0 50.0 46.98463103929542 17.101007166283424 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +85 112.9 15240 -56.5 -85.5 250.0 50.0 46.98463103929542 17.101007166283424 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +86 111.0 15349 -56.1 -86.1 248.0 45.0 41.72327345550543 16.857296703716052 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +87 107.6 15545 -56.9 -86.2 245.0 37.0 33.53338812035606 15.636875684405867 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +88 106.0 15641 -57.3 -86.3 244.0 37.0 33.255379713069175 16.219732431195876 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +89 103.0 15823 -56.5 -86.5 242.0 37.0 32.6690609357803 17.370447823077956 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +90 100.0 16010 -56.9 -85.9 240.0 37.0 32.04293994002422 18.500000000000018 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +91 93.5 16436 -56.3 -86.3 240.0 45.0 38.971143170299726 22.50000000000002 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +92 93.2 16459 -56.1 -86.2 240.0 45.0 38.971143170299726 22.50000000000002 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +93 88.3 16802 -53.5 -85.5 249.0 36.0 33.60889535389926 12.901246183630825 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +94 80.7 17374 -56.0 -86.7 265.0 20.0 19.92389396183491 1.743114854953165 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +95 77.1 17667 -57.3 -87.3 251.0 15.0 14.182778633989752 4.883522316857349 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +96 77.0 17678 -57.2 -87.3 250.0 15.0 14.095389311788626 5.130302149885027 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +97 75.3 17817 -55.9 -86.9 250.0 22.0 20.673237657289985 7.5244431531647065 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +98 73.3 17983 -56.2 -87.2 250.0 31.0 29.13047124436316 10.602624443095722 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +99 70.0 18280 -56.7 -87.7 250.0 21.0 19.733545036504076 7.182423009839038 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +100 69.9 18288 -56.7 -87.7 260.0 25.0 24.6201938253052 4.341204441673258 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +101 66.6 18593 -57.9 -87.9 215.0 29.0 16.63371665418034 23.755409284380757 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +102 66.5 18603 -57.9 -87.9 216.0 29.0 17.045772316481717 23.461492836873475 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +103 63.5 18898 -56.3 -87.4 230.0 36.0 27.577599952283204 23.140353948715426 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +104 61.3 19118 -55.1 -87.1 252.0 32.0 30.433808521444913 9.888543819998322 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +105 60.5 19202 -55.2 -87.0 260.0 31.0 30.529040343378448 5.3830935076748405 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +106 57.8 19492 -55.5 -86.5 241.0 31.0 27.113210921321272 15.029098227636442 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +107 57.7 19507 -55.4 -86.5 240.0 31.0 26.84678751731759 15.500000000000014 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +108 55.0 19812 -54.2 -86.0 250.0 33.0 31.00985648593498 11.286664729747061 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +109 52.4 20117 -52.9 -85.5 240.0 39.0 33.774990747593094 19.500000000000018 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +110 51.2 20267 -52.3 -85.3 245.0 33.0 29.908156972209454 13.94640263744307 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +111 50.0 20422 -52.9 -84.9 240.0 31.0 26.84678751731759 15.500000000000014 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +112 50.0 20420 -52.9 -84.9 250.0 27.0 25.37170076121953 9.234543869793049 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +113 47.7 20726 -53.0 -85.0 250.0 40.0 37.58770483143634 13.68080573302674 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +114 45.5 21031 -53.0 -85.0 275.0 24.0 23.908672754201895 -2.091737825943789 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +115 43.4 21336 -53.1 -85.1 255.0 19.0 18.352590699492296 4.917561856947892 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +116 42.0 21547 -53.1 -85.1 253.0 19.0 18.169790363297672 5.5550623897320035 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +117 39.5 21946 -50.9 -84.5 250.0 20.0 18.79385241571817 6.84040286651337 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +118 38.9 22045 -50.3 -84.3 258.0 20.0 19.56295201467611 4.158233816355196 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +119 37.7 22250 -50.8 -84.5 275.0 21.0 20.920088659926655 -1.8302705977008156 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +120 34.3 22860 -52.3 -85.2 250.0 16.0 15.035081932574535 5.472322293210696 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +121 32.7 23165 -53.0 -85.5 245.0 21.0 19.032463527769654 8.874983496554682 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +122 31.0 23518 -53.9 -85.9 267.0 25.0 24.965738368864347 1.3083989060736076 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +123 30.0 23730 -52.3 -85.3 280.0 27.0 26.58980933132962 -4.6885007970071095 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +124 29.8 23774 -52.2 -85.2 280.0 27.0 26.58980933132962 -4.6885007970071095 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +125 28.4 24079 -51.2 -84.7 245.0 26.0 23.564002462952903 10.988074805258178 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +126 27.3 24343 -50.3 -84.3 254.0 32.0 30.760374270026208 8.820395386143964 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +127 27.1 24384 -50.4 -84.4 255.0 33.0 31.875552267539256 8.54102848838318 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +128 25.9 24689 -51.2 -84.8 250.0 29.0 27.251086002791343 9.918584156444386 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +129 24.4 25072 -52.3 -85.3 269.0 29.0 28.995583159535347 0.5061197866812215 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +130 21.4 25908 -51.6 -84.6 310.0 29.0 22.215288850450367 -18.640840680909637 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +131 21.2 25983 -51.5 -84.5 308.0 28.0 22.06430110098821 -17.23852130911844 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +132 20.0 26360 -53.1 -85.1 300.0 23.0 19.918584287042087 -11.5 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +133 19.5 26518 -53.9 -85.4 300.0 23.0 19.918584287042087 -11.5 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +134 19.0 26689 -54.7 -85.7 303.0 31.0 25.998787606308152 -16.883810085465825 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +135 18.6 26822 -54.4 -85.6 305.0 38.0 31.127777682981687 -21.795904581339748 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +136 16.9 27432 -52.9 -85.3 305.0 36.0 29.489473594403705 -20.648751708637658 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +137 16.1 27755 -52.1 -85.1 305.0 37.0 30.308625638692696 -21.222328144988705 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +138 15.0 28211 -53.9 -85.9 305.0 38.0 31.127777682981687 -21.795904581339748 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +139 12.8 29234 -51.9 -84.9 305.0 41.0 33.58523381584866 -23.51663389039289 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +140 11.0 30206 -56.1 -87.1 305.0 43.0 35.22353790442665 -24.66378676309498 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +141 10.8 30323 -55.3 -86.3 305.0 44.0 36.04268994871564 -25.237363199446026 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +142 10.5 30480 -56.6 -87.0 305.0 44.0 36.04268994871564 -25.237363199446026 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +143 10.3 30623 -57.7 -87.7 305.0 45.0 36.86184199300463 -25.810939635797073 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +144 10.0 30810 -57.3 -88.3 305.0 46.0 37.68099403729362 -26.384516072148116 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +145 10.0 30785 -57.4 -88.2 305.0 46.0 37.68099403729362 -26.384516072148116 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 +146 9.4 31200 -58.7 -88.7 UIL 72797 2015-12-03 47.95 -124.55 62.0 20.93 diff --git a/samples/UIL_20151203_120000.txt b/samples/UIL_20151203_120000.txt new file mode 100644 index 0000000..13e7706 --- /dev/null +++ b/samples/UIL_20151203_120000.txt @@ -0,0 +1,111 @@ + pressure height temperature dewpoint direction speed u_wind v_wind station station_number time latitude longitude elevation pw +0 991.0 62 12.8 10.6 145.0 23.0 -13.192258036074056 18.840497018646815 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +1 973.0 215 14.0 12.1 154.0 32.0 -14.027876697250473 28.761409481573345 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +2 962.4 305 13.3 11.8 160.0 37.0 -12.654745303049749 34.76862696907861 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +3 927.7 610 11.0 10.7 165.0 45.0 -11.646857029613447 43.46666218300807 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +4 925.0 634 10.8 10.6 165.0 45.0 -11.646857029613447 43.46666218300807 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +5 910.0 771 10.0 9.9 175.0 51.0 -4.444942880130568 50.80592960267902 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +6 894.4 914 9.5 9.4 185.0 58.0 5.0550330793641605 57.77929248932124 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +7 862.2 1219 8.4 8.3 190.0 58.0 10.071594304681968 57.118849674708066 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +8 850.0 1337 8.0 7.9 190.0 56.0 9.724297949348106 55.14923416868365 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +9 829.0 1544 7.4 7.2 192.0 49.0 10.187672850070216 47.92923243595647 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +10 800.6 1829 5.6 5.4 195.0 40.0 10.352761804100831 38.63703305156273 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +11 771.3 2134 3.6 3.4 195.0 40.0 10.352761804100831 38.63703305156273 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +12 754.0 2319 2.4 2.2 195.0 40.0 10.352761804100831 38.63703305156273 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +13 743.0 2438 1.9 1.6 195.0 40.0 10.352761804100831 38.63703305156273 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +14 715.4 2743 0.5 0.2 195.0 38.0 9.83512371389579 36.70518139898459 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +15 700.0 2918 -0.3 -0.6 195.0 37.0 9.57630466879327 35.73925557269553 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +16 686.0 3080 -0.9 -1.3 197.0 37.0 10.817753074741262 35.38327597063231 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +17 637.8 3658 -4.6 -5.1 205.0 36.0 15.214257422665174 32.6270803333194 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +18 626.0 3805 -5.5 -6.1 203.0 38.0 14.847782882592394 34.979184431192735 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +19 601.0 4125 -6.5 -7.2 197.0 41.0 11.98723989363221 39.208494994484454 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +20 590.0 4267 -7.4 -8.2 195.0 43.0 11.129218939408394 41.53481053042994 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +21 545.1 4877 -11.3 -12.2 200.0 45.0 15.39090644965509 42.28616793536588 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +22 507.0 5435 -14.9 -16.0 204.0 48.0 19.52335886763841 43.85018196684484 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +23 500.0 5540 -15.7 -16.6 205.0 49.0 20.708294825294264 44.409081564795855 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +24 468.0 6035 -19.7 -20.7 205.0 56.0 23.66662265747916 50.7532360740524 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +25 464.1 6096 -20.2 -21.6 205.0 57.0 24.089240919219858 51.659543861089055 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +26 447.0 6374 -22.3 -25.4 204.0 62.0 25.21767187069961 56.63981837384125 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +27 401.0 7162 -27.9 -30.1 200.0 76.0 25.99353089275082 71.41663917972905 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +28 400.0 7180 -28.1 -30.4 200.0 76.0 25.99353089275082 71.41663917972905 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +29 394.0 7288 -28.9 -32.9 201.0 76.0 27.235964165442834 70.95211241378732 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +30 387.0 7416 -29.7 -33.7 203.0 77.0 30.086296893674064 70.87887371583791 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +31 376.0 7620 -31.6 -37.2 205.0 77.0 32.54160615403384 69.78569960182206 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +32 374.0 7658 -31.9 -37.9 205.0 77.0 32.54160615403384 69.78569960182206 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +33 340.0 8323 -37.9 -42.5 211.0 85.0 43.7782363673546 72.85922055967954 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +34 335.0 8424 -38.9 -44.9 212.0 86.0 45.57305672405561 72.93213626945264 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +35 325.0 8631 -40.5 -45.5 213.0 89.0 48.47287411633741 74.64168054714274 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +36 321.0 8715 -41.3 -48.3 214.0 90.0 50.3273613123672 74.61338152995377 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +37 317.0 8800 -42.3 -50.3 215.0 91.0 52.1954557079452 74.54283603029825 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +38 315.2 8839 -42.6 -50.6 215.0 91.0 52.1954557079452 74.54283603029825 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +39 301.2 9144 -45.1 -53.1 215.0 88.0 50.47472639889206 72.08537989743127 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +40 300.0 9170 -45.3 -53.3 215.0 88.0 50.47472639889206 72.08537989743127 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +41 269.0 9888 -51.9 -58.9 212.0 83.0 43.983298931355996 70.38799198098336 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +42 250.0 10360 -55.5 -61.5 210.0 79.0 39.50000000000001 68.41600689897065 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +43 228.0 10939 -60.7 -64.6 210.0 77.0 38.50000000000001 66.68395609140177 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +44 205.2 11582 -65.7 -69.7 210.0 75.0 37.50000000000001 64.9519052838329 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +45 200.0 11740 -66.9 -70.9 210.0 79.0 39.50000000000001 68.41600689897065 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +46 198.0 11801 -67.7 -71.8 210.0 82.0 41.00000000000001 71.01408311032397 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +47 195.2 11887 -67.3 -71.3 215.0 93.0 53.342608580647294 76.18114011887623 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +48 191.0 12019 -66.6 -70.6 215.0 93.0 53.342608580647294 76.18114011887623 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +49 189.0 12083 -66.3 -70.3 215.0 91.0 52.1954557079452 74.54283603029825 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +50 185.0 12214 -63.5 -68.5 215.0 86.0 49.32757352618997 70.44707580885328 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +51 180.0 12382 -63.3 -70.3 215.0 80.0 45.88611490808369 65.53216354311934 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +52 174.0 12593 -61.1 -71.1 215.0 73.0 41.87107985362637 59.798099233096394 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +53 171.0 12702 -58.3 -69.3 215.0 69.0 39.57677410822218 56.52149105594042 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +54 152.8 13411 -57.8 -71.4 215.0 45.0 25.810939635797077 36.861841993004624 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +55 150.0 13530 -57.7 -71.7 210.0 51.0 25.500000000000007 44.167295593006365 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +56 149.0 13572 -57.5 -71.5 210.0 53.0 26.500000000000007 45.89934640057525 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +57 143.0 13831 -58.7 -72.7 213.0 65.0 35.40153727597676 54.51358691645256 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +58 140.0 13965 -56.9 -70.9 215.0 71.0 40.723926980924276 58.15979514451841 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +59 139.0 14010 -56.6 -70.6 215.0 73.0 41.87107985362637 59.798099233096394 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +60 138.8 14021 -56.6 -70.6 215.0 73.0 41.87107985362637 59.798099233096394 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +61 137.0 14102 -56.1 -70.1 219.0 71.0 44.68174776453847 55.177363263444924 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +62 135.0 14196 -53.7 -68.7 224.0 70.0 48.62608593212982 50.353786023705574 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +63 126.1 14630 -55.0 -70.0 245.0 61.0 55.28477500923566 25.779713966182644 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +64 120.2 14935 -55.9 -70.9 240.0 38.0 32.90896534380866 19.000000000000018 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +65 120.0 14949 -55.9 -70.9 240.0 38.0 32.90896534380866 19.000000000000018 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +66 118.0 15055 -55.9 -69.9 236.0 39.0 32.33246532964663 21.808523235359118 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +67 115.0 15220 -54.5 -68.5 231.0 41.0 31.862984419735792 25.80213603304335 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +68 114.6 15240 -54.6 -68.6 230.0 41.0 31.407822167878095 26.354291997148124 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +69 111.0 15446 -55.5 -69.5 228.0 41.0 30.468937844573166 27.43435486071318 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +70 110.0 15503 -55.5 -69.5 227.0 41.0 29.98550176638599 27.96193276256244 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +71 100.0 16110 -55.9 -69.9 220.0 40.0 25.71150438746157 30.641777724759123 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +72 86.0 17069 -57.1 -71.1 225.0 25.0 17.677669529663685 17.677669529663692 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +73 83.9 17224 -57.3 -71.3 228.0 27.0 20.064910287889646 18.066526371689168 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +74 77.5 17729 -54.5 -68.5 236.0 35.0 29.016315039426466 19.571751621476132 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +75 74.1 18016 -55.3 -69.3 241.0 40.0 34.98478828557583 19.392384809853475 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +76 71.0 18288 -54.2 -67.5 245.0 44.0 39.877542629612606 18.59520351659076 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +77 70.0 18380 -53.9 -66.9 255.0 39.0 37.671107225273666 10.093942758998304 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +78 69.9 18389 -53.9 -66.9 255.0 38.0 36.70518139898459 9.835123713895785 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +79 67.7 18593 -53.0 -66.8 260.0 26.0 25.60500157831741 4.514852619340188 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +80 67.5 18614 -52.9 -66.8 258.0 25.0 24.45369001834514 5.197792270443995 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +81 64.6 18898 -53.5 -66.8 235.0 16.0 13.106432708623865 9.177222981616742 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +82 61.6 19202 -54.0 -66.7 220.0 19.0 12.212964584044245 14.554844419260585 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +83 58.7 19507 -54.6 -66.6 250.0 20.0 18.79385241571817 6.84040286651337 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +84 54.7 19964 -55.5 -66.5 248.0 30.0 27.81551563700362 11.23819780247737 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +85 50.8 20438 -53.7 -66.3 245.0 40.0 36.25231148146601 16.904730469627964 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +86 50.0 20540 -53.7 -66.3 245.0 42.0 38.06492705553931 17.749966993109364 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +87 48.6 20726 -53.7 -66.3 255.0 43.0 41.53481053042994 11.129218939408387 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +88 46.3 21031 -53.7 -66.2 260.0 38.0 37.42269461446391 6.598630751343353 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +89 44.2 21336 -53.7 -66.1 240.0 36.0 31.17691453623978 18.000000000000014 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +90 42.6 21569 -53.7 -66.0 236.0 40.0 33.161502902201676 22.367716138829863 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +91 42.1 21641 -53.4 -66.0 235.0 41.0 33.58523381584865 23.516633890392903 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +92 40.2 21946 -52.0 -65.9 270.0 30.0 30.0 5.510910596163089e-15 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +93 39.5 22057 -51.5 -65.9 268.0 28.0 27.98294315653468 0.9771859076700211 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +94 38.3 22250 -51.6 -65.8 265.0 24.0 23.908672754201895 2.091737825943798 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +95 34.9 22860 -51.7 -65.7 260.0 34.0 33.48346360241507 5.904038040675632 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +96 32.0 23423 -51.9 -65.5 258.0 30.0 29.344428022014167 6.237350724532795 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +97 28.2 24247 -49.5 -65.3 255.0 25.0 24.148145657226706 6.470476127563016 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +98 27.6 24384 -49.8 -65.3 255.0 24.0 23.18221983093764 6.211657082460495 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +99 26.3 24689 -50.6 -65.2 235.0 19.0 15.56388884149084 10.897952290669881 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +100 24.0 25298 -52.0 -65.0 265.0 24.0 23.908672754201895 2.091737825943798 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +101 22.9 25603 -52.8 -64.9 255.0 18.0 17.38666487320323 4.6587428118453715 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +102 20.8 26213 -54.3 -64.8 250.0 32.0 30.07016386514907 10.944644586421392 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +103 20.0 26470 -54.9 -64.7 260.0 30.0 29.544232590366242 5.2094453300079095 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +104 19.4 26664 -55.5 -64.6 266.0 33.0 32.9196136585742 2.301963633556144 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +105 17.2 27432 -52.5 -64.4 290.0 45.0 42.286167935365874 -15.390906449655107 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +106 16.9 27550 -52.1 -64.4 293.0 46.0 42.34322325881226 -17.97363191050658 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +107 15.7 28042 -52.6 -64.3 305.0 50.0 40.95760221444959 -28.678821817552304 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +108 14.9 28346 -52.8 -64.2 295.0 29.0 26.282925824062843 -12.255929590480289 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 +109 14.3 28629 -53.1 -64.1 UIL 72797 "2015-12-03 12:00:00" 47.95 -124.55 62.0 31.57 diff --git a/skewPy/SkewT.py b/skewPy/SkewT.py index 72d0a2d..9bc0639 100644 --- a/skewPy/SkewT.py +++ b/skewPy/SkewT.py @@ -18,7 +18,8 @@ from .thermodynamics import Rs_da, Cp_da, Epsilon -from collections import UserDict +#from collections import UserDict +from UserDict import UserDict from datetime import datetime import os,sys import scipy.interpolate as si @@ -465,7 +466,6 @@ def __init__(self, filename=None, data=None, fmt='UWYO', station_name=None, parc if data is None: self.data={} self.readfile(filename) - else: self.data=data self['SoundingDate']="" @@ -517,8 +517,9 @@ def interp_parcel(self): # let's try to interpolate the parcel stuff onto the environmental stuff - f_pt = si.interp1d(all_parcel_p, all_parcel_t, bounds_error=True) - + #print('Hello') + #f_pt = si.interp1d(all_parcel_p, all_parcel_t, bounds_error=True) + f_pt = si.interp1d(all_parcel_p, all_parcel_t, fill_value="extrapolate") self.parcel_p = self.data['pres'].copy() ###BD added .compressed here because I was getting an error about being a masked array in si.interp1d 3/2017 self.parcel_t = f_pt(self.parcel_p.compressed()) @@ -920,11 +921,26 @@ def readfile(self, fname): ndata = nlines-skip output = {} - fields = lines[3].split() - units = lines[4].split() + txtfields = lines[0].split() + fields = deepcopy(txtfields) + for ii, var in enumerate(fields): + if var == 'pressure': + fields[ii] = 'pres' + elif var == 'height': + fields[ii] = 'hght' + elif var == 'temperature': + fields[ii] = 'temp' + elif var == 'dewpoint': + fields[ii] = 'dwpt' + + #print(fields[0].lower()) + #fields = lines[3].split() + #units = lines[4].split() + #print(fields) + #print(units) # First line for WRF profiles differs from the UWYO soundings - header = lines[0] + header = lines[1] if header[:5] == '00000': # WRF profile self.station = '-99999' @@ -932,36 +948,75 @@ def readfile(self, fname): self['Latitude'] = float(header.split()[6]) self.sounding_date = header.split()[-1] else: - self.station = header[:5] - dstr = (' ').join(header.split()[-4:]) - self.sounding_date = datetime.strptime(dstr, "%HZ %d %b %Y").strftime("%Y-%m-%d_%H:%M:%S") - + #self.station = header[:5] + findstat = np.where([x.isalpha() for x in str(header)])[0] + self.station = header[min(findstat):max(findstat)+1] + #self.station = header[header.find(header.isalpha())] + #dstr = (' ').join(header.split()[-4:]) + dstr = (' ').join(header.split()[-5:-4]) + #self.sounding_date = datetime.strptime(dstr, "%HZ %d %b %Y").strftime("%Y-%m-%d_%H:%M:%S") + self.sounding_date = datetime.strptime(dstr, "%Y-%m-%d").strftime("%Y-%m-%d_%H:%M:%S") + if self.station_name is not None: self.station = self.station_name - for ff in fields: - output[ff.lower()]=zeros((nlines-skip)) - 999. - + + #print(output) #print 'output keys: {}'.format(output.keys()) - + lhi=[1, 9,16,23,30,37,46,53,58,65,72] rhi=[7,14,21,28,35,42,49,56,63,70,77] - - lcounter = 5 - for line,idx in zip(lines[6:],range(ndata)): - lcounter += 1 - + + ''' + #lcounter = 1 + #for line,idx in zip(lines[6:],range(ndata)): + for line,idx in zip(lines[1:],range(ndata)): + #lcounter += 1 + + #print(line) + #print(line[0],line[1],line[2]) + #print(line[lhi[0]:rhi[0]]) + #print('') try: output[fields[0].lower()][idx] = float(line[lhi[0]:rhi[0]]) except ValueError: + print('BREAK') break + #print(line) + #print(idx) + #print(output[fields[0].lower()][idx]) + + #input() for ii in range(1, len(rhi)): try: # Debug only: # print fields[ii].lower(), float(line[lhi[ii]:rhi[ii]].strip()) output[fields[ii].lower()][idx]=float(line[lhi[ii]:rhi[ii]].strip()) + #print(output[fields[ii].lower()][idx]) except ValueError: pass - + + #print(output) + #input() + ''' + # NEW! For loop to read in output data from sounding file + # First, determine number of columns, and find lines in the text file where data is missing ==> omit. May not be necessary if masking handles missing values properly... + numcols = len(lines[0].split())+1 + for line in lines[1:]: + if len(line.split()) != numcols: + lines.remove(line) + continue + + for ff in fields: + output[ff.lower()] = zeros((len(lines[1:])-1)) - 999. + + # Next, read each column into a dictionary. + for line, ii in zip(lines[1:],range(1,len(lines[1:]))): + for jj in range(0,numcols-1): + try: + output[fields[jj].lower()][ii-1] = float(line.split()[jj+1]) + except ValueError: + continue + for field in fields: #print field ff=field.lower() @@ -969,7 +1024,6 @@ def readfile(self, fname): self.data[ff]=ma.masked_values(output[ff], -999.) - elif self.fmt == 'EDT': # THIS IS FOR READING IN VAISALA EDITED DATA FILE FORMATS, FOR SPECIAL FIELD OPS SOUNDINGS rawfile = open(fname).readlines() serial_number = rawfile[1].split('#')[1].strip() @@ -1046,7 +1100,7 @@ def lift_parcel(self,startp,startt,startdp, plotkey = False): """ from numpy import interp - # print startp, startt, startdp + #print(startp, startt, startdp) assert startt >startdp, "Not a valid parcel. Check Td