from xml.dom import minidom from sys import stdin from subprocess import call import numpy as np from pyproj import Geod import matplotlib.patheffects as PathEffects from mpl_toolkits.basemap import Basemap from urllib.request import urlopen import matplotlib.colors as col import matplotlib.cm as cm import matplotlib.pyplot as plt import netCDF4 import sys sys.path.insert(0, '/home/djv/python/djvlibs/') #DEFINES DJV LIB PATH from dvpylib_shape import * import shapefile from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection from matplotlib.patches import PathPatch ''' #################################################################### RADAR2_REMOTE_SAVE_PLT.PY This script will retrieve a Nexrad level II file from NCEI NexRad BigData site and save locally in NetCDF or Shapefile format. User will define a STA DATE TIME OUTTYP PLTTYP and VAR as input parameters Script will retrieve a list of remote files for a given STA+DATE and then find the file closest to desired time. Will then retrieve and save as OUTTYP format in a user-defined location. Will then create plot of desired VAR and save as a PLTTYP file. Uses Weather Climate Toolkit(WCT) for data conversion http://www.ncdc.noaa.gov/wct/install.php Shapefile will generate several files(shp, sbf and shx files) Run: pyrun 6 radar2_remote_save_plt.py sta date time outtyp plttyp var ie pyrun 6 radar2_remote_save_plt.py KBRO 2017/08/26 003000 NC png Reflectivity outtyp: SHP or NC plttyp: png,ps,eps,jpg,screen,none=noplot var: "Reflectivity","Reflectivity_HI","RadialVelocity","RadialVelocity_HI" FUNCTIONS: pltmapE,pltmap, maximum, find_nearest_below, get_lonlat, getText #################################################################### ''' def pltmapE(wlon,elon,blat,ulat,inc,lalo): # Plot Map in cylindrical projection # Instead of using blank background, get a high resolution image from ESRI # lalo- 4dig tuple label plotting flag: (LRTB) # inc- increment for gridlines # wlon,elon,blat,ulat: map bounds(LR/UL corners) mp = Basemap(resolution='i',projection='cyl',\ llcrnrlon=wlon,llcrnrlat=blat,urcrnrlon=elon,urcrnrlat=ulat) maps = [ 'ESRI_Imagery_World_2D', # 0 # 'ESRI_StreetMap_World_2D', #1 # 'NatGeo_World_Map', # 2 # 'NGS_Topo_US_2D', # 3 'Ocean_Basemap', # 4 # 'USA_Topo_Maps', # 5 'World_Imagery', # 6 'World_Physical_Map', # 7 'World_Shaded_Relief', # 8 # 'World_Street_Map', # 9 'World_Terrain_Base', # 10 'World_Topo_Map' # 11 ] mp.arcgisimage(service=maps[2], xpixels = 1000, verbose= True) mp.drawcoastlines(zorder=1000) mp.drawstates(zorder=1000) mp.drawcounties(linewidth=.5,zorder=1000) # Set font for lalo labels font = {'family' : 'normal', 'size' : 8} plt.rc('font', **font) mp.drawparallels(np.arange(int(blat),int(ulat),inc),labels=lalo) mp.drawmeridians(np.arange(int(wlon),int(elon),inc),labels=lalo) plt.rcdefaults() #resets to default font sizes return mp def pltmap(wlon,elon,blat,ulat,inc,lalo): # Plot Map in cylindrical projection # lalo- 4dig tuple label plotting flag: (LRTB) # inc- increment for gridlines # wlon,elon,blat,ulat: map bounds(LR/UL corners) mp=Basemap(resolution='i',projection='cyl',\ llcrnrlon=wlon,llcrnrlat=blat,urcrnrlon=elon,urcrnrlat=ulat) # map.drawmapboundary(fill_color='aqua') # map.fillcontinents(color='#ddaa66',lake_color='aqua') mp.drawrivers(color='blue') mp.drawcoastlines(color='black',linewidth=1.5) mp.drawstates(color='black',linewidth=1.0) mp.drawcountries(color='black',linewidth=1.5) # Set font for lalo labels font = {'family' : 'normal', 'size' : 8} plt.rc('font', **font) mp.drawparallels(np.arange(int(blat),int(ulat),inc),labels=lalo) mp.drawmeridians(np.arange(int(wlon),int(elon),inc),labels=lalo) plt.rcdefaults() #resets to default font sizes return mp def maximum(items, default=None): # Finds maximum value in an array a=list(items.flat) mx=-999 mn=999 print("ll",len(a)) for x in range(len(a)): if (a[x] > mx and a[x]!= np.nan):mx = a[x] if (a[x] < mn and a[x]!= np.nan):mn = a[x] return mx,mn def find_nearest_below(array, value): # Will return index of array_val closest to but less than value # Return: index of nearest value below array = np.asarray(array) idx = (np.abs(array - value)).argmin() if value < array[idx]: idx=idx-1 return idx def get_lonlat(x,y): # Script to convert x/y values to proper latlon values for NexRad file # Inupt are x/y arrays, output lat/lon arrays # ie lon, lat = get_lon_lat(xlocs,ylocs) from pyproj import Proj p = Proj(proj='aeqd', ellps='sphere', lon_0=nexrad_lon,lat_0=nexrad_lat) return p(x,y,inverse=True) def getText(nodelist): rc = [] for node in nodelist: if node.nodeType == node.TEXT_NODE: rc.append(node.data) return ''.join(rc) #################################################################### # MAIN #################################################################### site = sys.argv[1] #station ID date = sys.argv[2] #date YYYY/MM/DD itime = sys.argv[3] #time HHMMSS out_typ = sys.argv[4] #Output type: NC or SHP plttyp = sys.argv[5] #Plot output typ(png,ps,eps,jpg,screen,none=noplot pltvar = sys.argv[6] #Var to plot: Reflectivity RadialVelocity ###################### USER DEFINED VARIABLES ####################### outpth = "./tmp/" #path to dump output varnames = ["Reflectivity","Reflectivity_HI","RadialVelocity","RadialVelocity_HI"] varsuff = ["R","R_HI","V","V_HI"] nrow=3 #nrows in panel plot ncol=3 #ncols in panel plot plotmap = "True" #plot w/map backgnd (True/False) cbar='True' #Plot color bar (True/False) cb_pth="/home11/staff/djv/python/" day=date.split('/') ofile0 = "NexRadII_"+site+day[0]+day[1]+day[2]+itime[0:4] ofile=ofile0+"."+plttyp print(ofile) ### Finds location of desired var in varname array #### loc=[] for i,v in enumerate(varnames): if v == pltvar: loc=i v0=-40 v1=75 cbmin=-35 cbmax=76 cbinc=5 cblab = "DbZ" cfil=cb_pth+"djv_radar_bin5.rgb" #change color table if velocity if loc >1: v0=-65. v1=65. cbmin=-50 cbmax=51 cbinc=10 cblab = "m/s" cfil=cb_pth+"djv_vad_bin5.rgb" ################## REMOTE LISTING OF DESIRED STA+DATE ############### wctpth = "/jm11/djv/radar/wct-4.3.1/" bucketURL = "http://noaa-nexrad-level2.s3.amazonaws.com" dirListURL = bucketURL+ "/?prefix=" + date + "/" + site xmldoc = minidom.parse(urlopen(dirListURL)) itemlist = xmldoc.getElementsByTagName('Key') print ("listing files from {} {} keys found".format(dirListURL,len(itemlist))) filelist=[] times=[] ############## IDENTIFY NAME OF REMOTE FILE TO RETRIEVE ############## for i,x in enumerate( itemlist): ifile = getText(x.childNodes) words=ifile.split("_") times.append(int(words[1])) filelist.append(ifile) print (i,"Found %s " % ifile) itim=find_nearest_below(times,int(itime)) print("Desired time=",date,itime) print("Closest file=",filelist[itim],itim) ############## Creates IR colormap ############### chan="RADAR" dvrad,ncolor=dvColormap(cfil,chan,"") #RGB colorfile,channel,path cmap1=col.LinearSegmentedColormap('my_colormap',dvrad,N=ncolor,gamma=1.) ################## FILE EXTRACTION AND CONVERSION ############### # Uses WCT to convert binary file to different format file = filelist[itim] name=file.split("/") filename="tmp/"+name[-1] ############ Converting to NetCDF +create plot ################# if out_typ == "NC": print("saving as netcdf file: ",filename) call(["sh", wctpth+"wct-export", "%s/%s"%(bucketURL,file), outpth, "rnc", wctpth+"wctBatchConfig.xml"]) if not plttyp in ["None","none"]: f = netCDF4.Dataset(filename+".nc",'r') data= np.array(f.variables[pltvar][:]) rng= np.array(f.variables['distance'+varsuff[loc]][:]) az= np.array(f.variables['azimuth'+varsuff[loc]][:]) elang= np.array(f.variables['elevation'+varsuff[loc]][:]) data = np.ma.array(data) #creates masked array omits vals=0 data = np.ma.masked_equal(data,0) # print(data[0,:,300]) print('elev= ', elang[:,100]) print("shapes",np.shape(data),np.shape(az),np.shape(rng),loc,varsuff[loc]) # Get map coordinates of station from Global file attributes btime = f.getncattr('time_coverage_start') clat = f.getncattr('StationLatitude') clon = f.getncattr('StationLongitude') wlon = int(clon-1.5) elon = int(clon+1.5) blat = int(clat-1.5) ulat = int(clat+1.5) blat = int(f.getncattr('geospatial_lat_min')) ulat = int(f.getncattr('geospatial_lat_max')+1.0) wlon = int(f.getncattr('geospatial_lon_min')-1.0) elon = int(f.getncattr('geospatial_lon_max')) VCP = f.getncattr('VolumeCoveragePattern') # Finds proper increment of lalo labels (~3) tinc=(max(abs(ulat-blat),abs(wlon-elon))//3) print(wlon,elon,blat,ulat,tinc,VCP) n_az0 =data.shape[1] n_rng0=data.shape[2] fig = plt.figure() cmap1 = col.LinearSegmentedColormap('my_colormap',dvrad,N=ncolor,gamma=1.) npan = min(data.shape[0],nrow*ncol) for sc in range(npan): # for sc in range(nrow*ncol): #LRTB ax = fig.add_subplot(nrow,ncol,sc+1) lalo_lab=[1,0,0,1] if sc>= (nrow-1)*ncol: lalo_lab[3]=1 else: lalo_lab[3]=0 ax.set_xticklabels(" ") if sc%ncol == 0 or ncol == 1: lalo_lab[0]=1 else: lalo_lab[0]=0 ax.set_yticklabels(" ") if plotmap == "False": #Plot data no map bckgnd print("SWEEP",sc,data.shape[0]) xlocs = rng/1000. * np.sin(np.deg2rad(az[sc,:, np.newaxis])) ylocs = rng/1000. * np.cos(np.deg2rad(az[sc,:, np.newaxis])) rad=ax.pcolormesh(xlocs, ylocs, data[sc,:,:], cmap=cmap1,vmin=v0 ,vmax=v1) ax.set_aspect('equal','box') #'datalim') ax.set_xlim(-400, 400) #creates plot +/- Nkm ax.set_ylim(-400, 400) plt.title('elev= '+str('%6.3f'% elang[sc,100])) print("LL",xlocs.shape,rng.shape,az.shape) else: # Convert az, range to a lat/lon on elliptical earth g = Geod(ellps='clrk66') # Type of ellipse the earth is projected on # There are other types of ellipses you can use, # but the differences will be small center_lat = np.ones([n_az0,n_rng0])*clat center_lon = np.ones([n_az0,n_rng0])*clon az2D = np.ones_like(center_lat)*az[sc,:,None] rng2D = np.ones_like(center_lat)*np.transpose(rng[:,None]) lon,lat,back=g.fwd(center_lon,center_lat,az2D,rng2D) mp=pltmap(wlon,elon,blat,ulat,tinc,lalo_lab) rad=ax.pcolormesh(lon,lat,data[sc,:,:],cmap=cmap1,vmin=v0 ,vmax=v1) plt.imshow(data[sc,:,:], interpolation='bilinear') # Plot station loc and id-label mp.scatter(clon,clat,color='darkblue',edgecolor='k',linewidth='1',s=20,zorder=500) plt.text(clon,clat,site,zorder=500,fontdict={'weight':'bold','color':'oldlace','size':4},path_effects=[PathEffects.withStroke(linewidth=1,foreground="k")]) plt.title('elev= '+str('%6.3f'% elang[sc,100])) if (cbar == 'True'): fig.subplots_adjust(bottom=0.1, top=0.8, left=0.1, right=0.8, wspace=0.1,hspace=0.2) #wsp=fraction axwidth cb_ax = fig.add_axes([0.93, 0.2, 0.02, 0.5]) cb=fig.colorbar(rad, cax=cb_ax)#orientation='vertical',cax=cb_ax) cbticlocs=[] cbticlabs=[] for it in range(cbmin,cbmax,cbinc): cbticlocs.append(it) cbticlabs.append(str(it)) cb.set_ticklabels(cbticlabs,update_ticks=True) cb.set_ticks(cbticlocs,update_ticks=True) cb.set_label(cblab) cb.ax.tick_params(labelsize=12) cb.ax.yaxis.set_ticks_position('left') plt.suptitle("Level II Radar all angles "+site+" "+btime+" VCP="+str(VCP)) if plttyp in ["Screen","screen"]: plt.show() else: fig.savefig(ofile,format=plttyp,orientation='portrait',quality=100) else: ############ Converting to Shapfile +create plot ############### print("saving as shape file: ",filename) # call(["sh", wctpth+"wct-export", "%s/%s"%(bucketURL,file), outpth, "shp", wctpth+"wctBatchConfig.xml"]) ###### Reads+lists shapefile info ####### sf = shapefile.Reader(filename) shape_info(sf) shapes = sf.shapes() recs=sf.records() #records are list of values corresponding to fields Nshp = len(recs) ############## Plot shape onto Map ######## if not plttyp in ["None","none"]: fig = plt.figure() ax = fig.add_subplot(111) area = sf.bbox # map = Basemap(resolution='i',projection='cyl',\ # llcrnrlon=area[0],llcrnrlat=area[1],urcrnrlon=area[2],urcrnrlat=area[3]) # map.drawmapboundary(fill_color='aqua') # map.fillcontinents(color='#ddaa66',lake_color='aqua') # map.drawcoastlines() # map.readshapefile(filename,'value', drawbounds = False) mp = Basemap(resolution='i',projection='cyl',\ llcrnrlon=area[0],llcrnrlat=area[1],urcrnrlon=area[2],urcrnrlat=area[3]) mp.drawmapboundary(fill_color='aqua') mp.fillcontinents(color='#ddaa66',lake_color='aqua') mp.drawcoastlines() mp.readshapefile(filename,'value', drawbounds = False) patchesa = [] patchesb = [] # for info, shape in zip(map.value_info, map.value): # if info['value'] >= 20: # patchesa.append( Polygon(np.array(shape), True) ) # else: # patchesb.append( Polygon(np.array(shape), True) ) iskp=1000 #plot every nth point for nshp in range(0,Nshp,iskp): # vals = [] pts = np.array(shapes[nshp].points) prt = shapes[nshp].parts par = list(prt) + [pts.shape[0]] tmp=recs[nshp][3] cindx=cmap1(int(tmp)+40) for pij in range(len(prt)): vals.append(Polygon(pts[par[pij]:par[pij+1]])) ax.add_collection(PatchCollection(vals, facecolor= cindx, edgecolor=cindx, linewidths=1., zorder=2)) ######## two color patches ##### # ax.add_collection(PatchCollection(patchesa, facecolor= 'r', edgecolor='r', linewidths=1., zorder=2)) # ax.add_collection(PatchCollection(patchesb, facecolor= 'g', edgecolor='g', linewidths=1., zorder=2)) plt.title(site+" "+recs[0][1]) if plttyp in ["Screen","screen"]: plt.show() else: fig.savefig(ofile,format=plttyp,orientation='portrait',quality=100) exit()