import numpy as np import os from netCDF4 import Dataset as NetCDFFile import drawing_tools as drw from optparse import OptionParser import sys from cStringIO import StringIO ## e.g. # drawing.py -f /media/data2/etudes/WRF_LMDZ/WL_HyMeX/IIphase/medic950116/wlmdza/wrfout/wrfout_d01_1995-01-13_00:00:00 -o create_movie -S 'draw_2D_shad#Time@WRFTimes@10@95@191@1#tas:East_West|-1,North_South|-1,Time|2:longitude:latitude:Summer:270.,300.:tas|at|t=0:pdf:None:None' -v T2 ## e.g. # drawing.py -f wrfout_d01_1980-03-01_00\:00\:00_Time_B0-E48-I1.nc -o draw_2D_shad -S 'tas:East_West|-1,North_South|-1,Time|2:longitude:latitude:Summer:270.,300.:tas|at|t=0:pdf:None:cyl,i' -v T2 ## e.g. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'landmask,height:Time|0:Time|0:XLONG_M:XLAT_M:rainbow:fixc,k:%.2f:0,1:0.,3000.,10:landmask & height:pdf:False:lcc,i' -v LANDMASK,HGT_M ## e.g. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'height,landmask:Time|0:Time|0:XLONG_M:XLAT_M:terrain:fixc,k:None:0.,3000.:0,1,10:MEDCORDEX height & landmask:pdf:False:lcc,i' -v HGT_M,LANDMASK ## e.g. # drawing.py -o draw_2D_shad_line -f 'mean_dtcon-pluc-pres_lat.nc,mean_dtcon-pluc-pres_lat.nc' -S 'dtcon,prc:bottom_top|-1,south_north|-1:latmean:presmean:seismic,k:-5.,5.:monthly|dtcon|&|prc:pdf:flip@y:None:True' -v 'dtconmean,prcmean' ## e.g. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i' ## e.g. # drawing.py -o draw_topo_geogrid_boxes -S '0.,3000.:None:FF domains:pdf:lcc,i:d01,d02,d03' -f geo_em.d01.nc,geo_em.d02.nc,geo_em.d03.nc ## e.g. # drawing.py -o draw_trajectories -f 'WRF/control/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_HighRes_C/geo_em.d03.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdza/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb_ii/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M' -S '$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$,$LMDZ_{NPv3.1b}$|None|medicane trajectories|pdf|cyl,i' -v obs/trajectory.dat,satellite,-1 ## e.g. # drawing.py -o draw_vals_trajectories -f WRF_LMDZ/wlmdza/tevolboxtraj_T2.nc,WRF_LMDZ/wlmdzb/tevolboxtraj_T2.nc,WRF/control/tevolboxtraj_T2.nc -S 'mean:-1:$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$@4:tas:time|($[DD]^[HH]$):exct,6,h:$%d^{%H}$:trajectory|following|mean:pdf' -v T2 ## e.g. # drawing.py -o draw_2D_shad_time -f 'netcdf_concatenated.nc' -S 'dtcon:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True' -v 'dtconmean' ## e.g. # drawing.py -o variable_values -S PSFC ## e.g. # drawing.py -o draw_timeSeries -f wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc -S 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf' -v 'LDQCON,time' ## e.g. # drawing.py -f wrfout_d01_1979-12-01_00:00:00 -o draw_Neighbourghood_evol -S 'q:Time|-1|Times,bottom_top|6|ZNU,south_north|3|XLAT,west_east|26|XLONG:south_north,west_east:5:auto:time|($[DD]^{[HH]}$),time|($[DD]^{[HH]}$):exct,2,h|exct,1,d:$%d^{%H}$,$%d^{%H}$:5|pts|neighbourghood|temporal|evolution:0.0,0.004:BuPu:pdf:True' -v QVAPOR ## e.g. # drawing.py -o draw_lines_time -f wrfout_d01_1980-03-01_00:00:00_Time_B0-E48-I1_south_north_B15-E15-I1_west_east_B15-E15-I1.nc -S 'time;y;time ([DD]${[HH]}$);file1;tas;evolution;time|hours!since!1949-12-01_00:00:00|exct,12,h|%d$^{%H}$;pdf' -v T2 main = 'drawing.py' errormsg = 'ERROR -- error -- ERROR -- error' warnmsg = 'WARNING -- waring -- WARNING -- warning' fillValue=1.e20 namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time', \ 'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line', \ 'draw_2D_shad_line_time', 'draw_timeSeries', 'draw_topo_geogrid', \ 'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \ 'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics', \ 'variable_values'] def draw_2D_shad(ncfile, values, varn): """ plotting a fields with shading draw_2D_shad(ncfile, values, varn) ncfile= file to use values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]: [kindfig]:[reverse]:[mapv]:[close] [vnamefs]: Name in the figure of the variable to be shaded [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the variable a given value is required (-1, all the length) [dimx/yvn]: name of the variables with the values of the final dimensions (x,y) [colorbar]: name of the color bar [smin/axv]: minimum and maximum value for the shading or: 'Srange': for full range 'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean) 'Saroundminmax@val': for min*val,max*val 'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val), percentile_(100-val)-median) 'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean) 'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median) 'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val), percentile_(100-val)-median) [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: Transformation of the values * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y [mapv]: map characteristics: [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lamvbert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None' varn= [varsn] name of the variable to plot with shading """ fname = 'draw_2D_shad' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad.__doc__ quit() expectargs = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \ '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[mapv]', '[close]'] drw.check_arguments(fname,len(expectargs),values,':',expectargs) vnamesfig = values.split(':')[0] dimvals= values.split(':')[1].replace('|',':') vdimxn = values.split(':')[2] vdimyn = values.split(':')[3] colbarn = values.split(':')[4] shadminmax = values.split(':')[5] figtitle = values.split(':')[6].replace('|',' ') figkind = values.split(':')[7] revals = values.split(':')[8] mapvalue = values.split(':')[9] # varn = values.split(':')[10] ncfiles = ncfile if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') varns = varn.split(',')[0] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) # Dimensions names ## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(',')) ## dimnamesv = [] ## for idd in range(len(objvars.dimensions)): ## cutdim = False ## for idc in range(len(dimvals.split(','))): ## dimcutn = dimvals.split(',')[idc].split(':')[0] ## print objvars.dimensions[idd], dimcutn ## if objvars.dimensions[idd] == dimcutn: ## cutdim = True ## break ## if not cutdim: dimnamesv.append(objvars.dimensions[idd]) dimnamesv = [vdimyn, vdimxn] if drw.searchInlist(objvars.ncattrs(),'units'): varunits = objvars.getncattr('units') else: print warnmsg print ' ' + fname + ": variable '" + varn + "' without units!!" varunits = '-' if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimyn + '" !!' quit(-1) objdimx = objsf.variables[vdimxn] objdimy = objsf.variables[vdimyn] if drw.searchInlist(objdimx.ncattrs(),'units'): odimxu = objdimx.getncattr('units') else: print warnmsg print ' ' + fname + ": variable dimension '" + vdimxn + "' without units!!" odimxu = '-' if drw.searchInlist(objdimy.ncattrs(),'units'): odimyu = objdimy.getncattr('units') else: print warnmsg print ' ' + fname + ": variable dimension '" + vdimyn + "' without units!!" odimyu = '-' if len(objdimx.shape) <= 2: # odimxv = objdimx[valshad.shape] # odimyv = objdimy[valshad.shape] odimxv = objdimx[:] odimyv = objdimy[:] elif len(objdimx.shape) == 3: # dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])] # odimxv = objdimx[tuple(dimcut)] # odimyv = objdimy[tuple(dimcut)] odimxv = objdimx[0,:] odimyv = objdimy[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \ ' not ready!!' quit(-1) shading_nx = [] if shadminmax.split(',')[0][0:1] != 'S': shading_nx.append(np.float(shadminmax.split(',')[0])) else: shading_nx.append(shadminmax.split(',')[0]) if shadminmax.split(',')[1][0:1] != 'S': shading_nx.append(np.float(shadminmax.split(',')[1])) else: shading_nx.append(shadminmax.split(',')[1]) if mapvalue == 'None': mapvalue = None drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\ colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True) return def draw_2D_shad_time(ncfile, values, varn): """ plotting a fields with shading with time values draw_2D_shad(ncfile, values, varn) ncfile= file to use values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]: [kindfig]:[reverse]:[timevals]:[close] [vnamefs]: Name in the figure of the variable to be shaded [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the variable a given value is required (-1, all the length) [dimx/yvn]: name of the variables with the values of the final dimensions (x,y) [colorbar]: name of the color bar [smin/axv]: minimum and maximum value for the shading or: 'Srange': for full range 'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean) 'Saroundminmax@val': for min*val,max*val 'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val), percentile_(100-val)-median) 'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean) 'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median) 'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val), percentile_(100-val)-median) [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: Transformation of the values * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics [timen]; name of the time variable [units]; units string according to CF conventions ([tunits] since [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces) [kind]; kind of output 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [tfmt]; desired format [label]; label at the graph ('!' for spaces) [close]: should figure be closed (finished) values='dtcon:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|' 'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True varn= [varsn] name of the variable to plot with shading """ fname = 'draw_2D_shad_time' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad_time.__doc__ quit() farguments = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \ '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', '[close]'] drw.check_arguments(fname,11,values,':',farguments) vnamesfig = values.split(':')[0] dimvals= values.split(':')[1].replace('|',':') vdimxn = values.split(':')[2] vdimyn = values.split(':')[3] colbarn = values.split(':')[4] shadminmax = values.split(':')[5] figtitle = values.split(':')[6].replace('|',' ') figkind = values.split(':')[7] revals = values.split(':')[8] timevals = values.split(':')[9] close = values.split(':')[10] ncfiles = ncfile if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') varns = varn.split(',')[0] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) dimnamesv = [vdimyn, vdimxn] varunits = objvars.getncattr('units') if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimensino variable "' + vdimyn + '" !!' quit(-1) objdimx = objsf.variables[vdimxn] objdimy = objsf.variables[vdimyn] odimxu = objdimx.getncattr('units') odimyu = objdimy.getncattr('units') if len(objdimx.shape) <= 2: odimxv = objdimx[:] odimyv = objdimy[:] elif len(objdimx.shape) == 3: odimxv = objdimx[0,:] odimyv = objdimy[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \ ' not ready!!' quit(-1) timename = timevals.split('|')[0] timeunit = timevals.split('|')[1].replace('!',' ') timekind = timevals.split('|')[2] timefmt = timevals.split('|')[3] timelabel = timevals.split('|')[4].replace('!',' ') if vdimxn == timename: odimxv = objsf.variables[vdimxn][:] odimxu = timelabel timeaxis = 'x' odimyv = objsf.variables[vdimyn] odimyu = odimyv.getncattr('units') timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt) elif vdimyn == timename: odimyv = objsf.variables[vdimyn][:] odimyu = timelabel timeaxis = 'y' odimxv = objsf.variables[vdimxn] odimxu = odimxv.getncattr('units') timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt) else: print errormsg print ' ' + fname + ": time variable '" + timename + "' not found!!" quit(-1) shading_nx = [] if shadminmax.split(',')[0][0:1] != 'S': shading_nx.append(np.float(shadminmax.split(',')[0])) else: shading_nx.append(shadminmax.split(',')[0]) if shadminmax.split(',')[1][0:1] != 'S': shading_nx.append(np.float(shadminmax.split(',')[1])) else: shading_nx.append(shadminmax.split(',')[1]) closeval = drw.Str_Bool(close) drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \ dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \ timepos, timelabels, closeval) return def draw_2D_shad_cont(ncfile, values, varn): """ plotting two fields, one with shading and the other with contour lines draw_2D_shad_cont(ncfile, values, varn) ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file) values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv] [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables [dimvals/c]: list of [dimname]|[value] telling at which dimension of the variable a given value is required (no dimension name, all the length) [dimx/yvn]: ',' list with the name of the variables with the values of the dimensions [colorbar]: name of the color bar [ckind]: kind of contours 'cmap': as it gets from colorbar 'fixc,[colname]': fixed color [colname], all stright lines 'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line [clabfmt]: format of the labels in the contour (None, also possible) [smin/axv]: minimum and maximum value for the shading [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: does the values be transposed? 'True/False', [mapv]: map characteristics: [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lamvbert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full valules= 'rh,ta:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal|average|of|rh|&|ta:pdf:flip@y:None' varn= [varsn],[varcn] name of the variable to plot with shading variable with contour """ fname = 'draw_2D_shad_cont' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad_cont.__doc__ quit() vnamesfig = values.split(':')[0].split(',') dimvals= values.split(':')[1].replace('|',':') dimvalc= values.split(':')[2].replace('|',':') vdimxn = values.split(':')[3] vdimyn = values.split(':')[4] colbarn = values.split(':')[5] countkind = values.split(':')[6] countlabelfmt = values.split(':')[7] shadminmax = values.split(':')[8] contlevels = values.split(':')[9] figtitle = values.split(':')[10].replace('|',' ') figkind = values.split(':')[11] revals = values.split(':')[12] mapvalue = values.split(':')[13] if2filenames = ncfile.find(',') if if2filenames != -1: ncfiles = ncfile.split(',')[0] ncfilec = ncfile.split(',')[1] else: ncfiles = ncfile ncfilec = ncfile if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) if not os.path.isfile(ncfilec): print errormsg print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') objcf = NetCDFFile(ncfilec, 'r') varns = varn.split(',')[0] varnc = varn.split(',')[1] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) if not objcf.variables.has_key(varnc): print errormsg print ' ' + fname + ': contour file "' + ncfilec + \ '" does not have variable "' + varnc + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] objvarc = objcf.variables[varnc] if len(objvars.shape) != len(objvarc.shape): print errormsg print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \ objvars.shape, 'different than contour variable "' + varnc + '": ', \ objvarc.shape,' !!!' quit(-1) for idim in range(len(objvars.shape)): if objvars.shape[idim] != objvarc.shape[idim]: print errormsg print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \ objvars.shape, 'different than contour variable "' + varnc + '": ', \ objvarc.shape,' !!!' quit(-1) valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|')) # Dimensions names ## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(',')) ## dimnamesv = [] ## for idd in range(len(objvars.dimensions)): ## cutdim = False ## for idc in range(len(dimvals.split(','))): ## dimcutn = dimvals.split(',')[idc].split(':')[0] ## print objvars.dimensions[idd], dimcutn ## if objvars.dimensions[idd] == dimcutn: ## cutdim = True ## break ## if not cutdim: dimnamesv.append(objvars.dimensions[idd]) dimnamesv = [vdimyn, vdimxn] varunits = [] varunits.append(objvars.getncattr('units')) varunits.append(objvarc.getncattr('units')) if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimensino variable "' + vdimyn + '" !!' quit(-1) objdimx = objsf.variables[vdimxn] objdimy = objsf.variables[vdimyn] odimxu = objdimx.getncattr('units') odimyu = objdimy.getncattr('units') if len(objdimx.shape) <= 2: # odimxv = objdimx[valshad.shape] # odimyv = objdimy[valshad.shape] odimxv = objdimx[:] odimyv = objdimy[:] elif len(objdimx.shape) == 3: # dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])] # odimxv = objdimx[tuple(dimcut)] # odimyv = objdimy[tuple(dimcut)] odimxv = objdimx[0,:] odimyv = objdimy[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \ ' not ready!!' quit(-1) if countlabelfmt == 'None': countlfmt = None else: countlfmt = countlabelfmt shading_nx = np.zeros((2), dtype=np.float) shading_nx[0] = np.float(shadminmax.split(',')[0]) shading_nx[1] = np.float(shadminmax.split(',')[1]) clevmin = np.float(contlevels.split(',')[0]) clevmax = np.float(contlevels.split(',')[1]) Nclevels = int(contlevels.split(',')[2]) levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels) if len(levels_cont) <= 1: print warnmsg print ' ' + fname + ': wrong contour levels:', levels_cont,' !!' del(levels_cont) levels_cont = np.zeros((Nclevels), dtype=np.float) levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1) print ' generating default ones: ',levels_cont if mapvalue == 'None': mapvalue = None drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu, \ odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont, \ varunits, figtitle, figkind, revals, mapvalue) return def draw_2D_shad_cont_time(ncfile, values, varn): """ plotting two fields, one with shading and the other with contour lines being one of the dimensions of time characteristics draw_2D_shad_cont(ncfile, values, varn) ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file) values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv] [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables [dimvals/c]: list of [dimname]|[value] telling at which dimension of the variable a given value is required (no dimension name, all the length) [dimx/yvn]: ',' list with the name of the variables with the values of the dimensions [colorbar]: name of the color bar [ckind]: kind of contours 'cmap': as it gets from colorbar 'fixc,[colname]': fixed color [colname], all stright lines 'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line [clabfmt]: format of the labels in the contour (None, also possible) [smin/axv]: minimum and maximum value for the shading [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: modification to the dimensions: 'transposed': transpose matrices 'flip',[x/y]: flip only the dimension [x] or [y] [timevals]: [timen]|[units]|[kind]|[tfmt]|[label] time labels characteristics [timen]; name of the time variable [units]; units string according to CF conventions ([tunits] since [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces) [kind]; kind of output 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [tfmt]; desired format [label]; label at the graph ('!' for spaces) [mapv]: map characteristics: [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lamvbert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full valules= 'rh,ta:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal|average|of|rh|&|ta:pdf:flip@y:time!hours!since!1949/12/01|exct,5d|%d|date!([DD]):None' varn= [varsn],[varcn] name of the variable to plot with shading variable with contour """ fname = 'draw_2D_shad_cont_time' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad_cont_time.__doc__ quit() expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]', \ '[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]', \ '[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', \ '[mapv]'] drw.check_arguments(fname,len(expectargs),values,':',expectargs) vnamesfig = values.split(':')[0].split(',') dimvals= values.split(':')[1].replace('|',':') dimvalc= values.split(':')[2].replace('|',':') vdimxn = values.split(':')[3] vdimyn = values.split(':')[4] colbarn = values.split(':')[5] countkind = values.split(':')[6] countlabelfmt = values.split(':')[7] shadminmax = values.split(':')[8] contlevels = values.split(':')[9] figtitle = values.split(':')[10].replace('|',' ') figkind = values.split(':')[11] revals = values.split(':')[12] timevals = values.split(':')[13] mapvalue = values.split(':')[14] if2filenames = ncfile.find(',') if if2filenames != -1: ncfiles = ncfile.split(',')[0] ncfilec = ncfile.split(',')[1] else: ncfiles = ncfile ncfilec = ncfile if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) if not os.path.isfile(ncfilec): print errormsg print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') objcf = NetCDFFile(ncfilec, 'r') varns = varn.split(',')[0] varnc = varn.split(',')[1] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) if not objcf.variables.has_key(varnc): print errormsg print ' ' + fname + ': contour file "' + ncfilec + \ '" does not have variable "' + varnc + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] objvarc = objcf.variables[varnc] if len(objvars.shape) != len(objvarc.shape): print errormsg print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \ objvars.shape, 'different than contour variable "' + varnc + '": ', \ objvarc.shape,' !!!' quit(-1) for idim in range(len(objvars.shape)): if objvars.shape[idim] != objvarc.shape[idim]: print errormsg print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \ objvars.shape, 'different than contour variable "' + varnc + '": ', \ objvarc.shape,' !!!' quit(-1) valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|')) # Dimensions names ## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(',')) ## dimnamesv = [] ## for idd in range(len(objvars.dimensions)): ## cutdim = False ## for idc in range(len(dimvals.split(','))): ## dimcutn = dimvals.split(',')[idc].split(':')[0] ## print objvars.dimensions[idd], dimcutn ## if objvars.dimensions[idd] == dimcutn: ## cutdim = True ## break ## if not cutdim: dimnamesv.append(objvars.dimensions[idd]) dimnamesv = [vdimyn, vdimxn] varunits = [] varunits.append(objvars.getncattr('units')) varunits.append(objvarc.getncattr('units')) if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimyn + '" !!' quit(-1) timename = timevals.split('|')[0] timeunit = timevals.split('|')[1].replace('!',' ') timekind = timevals.split('|')[2] timefmt = timevals.split('|')[3] timelabel = timevals.split('|')[4].replace('!',' ') if vdimxn == timename: timevals = objsf.variables[vdimxn][:] dimt = 'x' ovalaxis = objsf.variables[vdimyn] ovalu = ovalaxis.getncattr('units') elif vdimyn == timename: timevals = objsf.variables[vdimyn][:] dimt = 'y' ovalaxis = objsf.variables[vdimxn] ovalu = ovalaxis.getncattr('units') else: print errormsg print ' ' + fname + ": time variable '" + timename + "' not found!!" quit(-1) timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt) if len(ovalaxis.shape) <= 2: ovalaxisv = ovalaxis[:] elif len(ovalaxis.shape) == 3: ovalaxisv = ovalaxis[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', ovalaxis.shape, \ ' not ready!!' quit(-1) if countlabelfmt == 'None': countlfmt = None else: countlfmt = countlabelfmt shading_nx = np.zeros((2), dtype=np.float) shading_nx[0] = np.float(shadminmax.split(',')[0]) shading_nx[1] = np.float(shadminmax.split(',')[1]) clevmin = np.float(contlevels.split(',')[0]) clevmax = np.float(contlevels.split(',')[1]) Nclevels = int(contlevels.split(',')[2]) levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels) if len(levels_cont) <= 1: print warnmsg print ' ' + fname + ': wrong contour levels:', levels_cont,' !!' del(levels_cont) levels_cont = np.zeros((Nclevels), dtype=np.float) levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1) print ' generating default ones: ',levels_cont if mapvalue == 'None': mapvalue = None drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv, \ timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn, \ countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, \ revals, mapvalue) return def draw_2D_shad_line(ncfile, values, varn): """ plotting a fields with shading and another with line draw_2D_shad_line(ncfile, values, varn) ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]: [kindfig]:[reverse]:[mapv]:[close] [vnamefs]: Name in the figure of the variable to be shaded [vnamefl]: Name in the figure of the variable to be lined [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the variable a given value is required (-1, all the length) [dimx/yvn]: name of the variables with the values of the final dimensions (x,y) [colorbar]: name of the color bar [colline]: name of the color for the line [smin/axv]: minimum and maximum value for the shading [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: Transformation of the values * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y [mapv]: map characteristics: [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lamvbert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None' varn= [varsn],[varnl] name of the variable to plot with shading and with line """ fname = 'draw_2D_shad_line' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad_line.__doc__ quit() farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \ '[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \ '[mapv]', '[close]'] drw.check_arguments(fname,11,values,':',farguments) vnamesfig = values.split(':')[0].split(',')[0] dimvals= values.split(':')[1].replace('|',':') vdimxn = values.split(':')[2] vdimyn = values.split(':')[3] colbarn = values.split(':')[4].split(',')[0] shadminmax = values.split(':')[5] figtitle = values.split(':')[6].replace('|',' ') figkind = values.split(':')[7] revals = values.split(':')[8] mapvalue = values.split(':')[9] # varn = values.split(':')[10] ncfiles = ncfile.split(',')[0] if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') varns = varn.split(',')[0] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) # Dimensions names ## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(',')) ## dimnamesv = [] ## for idd in range(len(objvars.dimensions)): ## cutdim = False ## for idc in range(len(dimvals.split(','))): ## dimcutn = dimvals.split(',')[idc].split(':')[0] ## print objvars.dimensions[idd], dimcutn ## if objvars.dimensions[idd] == dimcutn: ## cutdim = True ## break ## if not cutdim: dimnamesv.append(objvars.dimensions[idd]) dimnamesv = [vdimyn, vdimxn] varunits = objvars.getncattr('units') if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimensino variable "' + vdimyn + '" !!' quit(-1) objdimx = objsf.variables[vdimxn] objdimy = objsf.variables[vdimyn] odimxu = objdimx.getncattr('units') odimyu = objdimy.getncattr('units') if len(objdimx.shape) <= 2: # odimxv = objdimx[valshad.shape] # odimyv = objdimy[valshad.shape] odimxv = objdimx[:] odimyv = objdimy[:] elif len(objdimx.shape) == 3: # dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])] # odimxv = objdimx[tuple(dimcut)] # odimyv = objdimy[tuple(dimcut)] odimxv = objdimx[0,:] odimyv = objdimy[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \ ' not ready!!' quit(-1) shading_nx = np.zeros((2), dtype=np.float) shading_nx[0] = np.float(shadminmax.split(',')[0]) shading_nx[1] = np.float(shadminmax.split(',')[1]) if mapvalue == 'None': mapvalue = None # line plot ## ncfilel = ncfile.split(',')[1] vnamelfig = values.split(':')[0].split(',')[1] varnl = varn.split(',')[1] colline = values.split(':')[4].split(',')[1] objlf = NetCDFFile(ncfilel,'r') objlvar = objlf.variables[varnl] linevals = objlvar[:] varlunits = objlvar.units drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \ odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits, \ figtitle, figkind, revals, mapvalue, True) objsf.close() objlf.close() return def draw_2D_shad_line_time(ncfile, values, varn): """ plotting a fields with shading and a line with time values draw_2D_shad_line(ncfile, values, varn) ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]: [kindfig]:[reverse]:[timevals]:[close] [vnamefs]: Name in the figure of the variable to be shaded [vnamefl]: Name in the figure of the variable to be lined [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the variable a given value is required (-1, all the length) [dimx/yvn]: name of the variables with the values of the final dimensions (x,y) [colorbar]: name of the color bar [smin/axv]: minimum and maximum value for the shading [figt]: title of the figure ('|' for spaces) [kindfig]: kind of figure [reverse]: Transformation of the values * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics [timen]; name of the time variable [units]; units string according to CF conventions ([tunits] since [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces) [kind]; kind of output 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [tfmt]; desired format [label]; label at the graph ('!' for spaces) [close]: should figure be closed (finished) values='dtcon,prc:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|' 'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True varn= [varsn].[varln] name of the variable to plot with shading and to plot with line """ fname = 'draw_2D_shad_line_time' if values == 'h': print fname + '_____________________________________________________________' print draw_2D_shad__line_time.__doc__ quit() farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \ '[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \ '[timevals]', '[close]'] drw.check_arguments(fname,11,values,':',farguments) vnamesfig = values.split(':')[0].split(',')[0] dimvals= values.split(':')[1].replace('|',':') vdimxn = values.split(':')[2] vdimyn = values.split(':')[3] colbarn = values.split(':')[4] shadminmax = values.split(':')[5] figtitle = values.split(':')[6].replace('|',' ') figkind = values.split(':')[7] revals = values.split(':')[8] timevals = values.split(':')[9] close = values.split(':')[10] ncfiles = ncfile.split(',')[0] if not os.path.isfile(ncfiles): print errormsg print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!' quit(-1) objsf = NetCDFFile(ncfiles, 'r') varns = varn.split(',')[0] if not objsf.variables.has_key(varns): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have variable "' + varns + '" !!' quit(-1) # Variables' values objvars = objsf.variables[varns] valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|')) dimnamesv = [vdimyn, vdimxn] varunits = objvars.getncattr('units') if not objsf.variables.has_key(vdimxn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimension variable "' + vdimxn + '" !!' quit(-1) if not objsf.variables.has_key(vdimyn): print errormsg print ' ' + fname + ': shading file "' + ncfiles + \ '" does not have dimensino variable "' + vdimyn + '" !!' quit(-1) objdimx = objsf.variables[vdimxn] objdimy = objsf.variables[vdimyn] odimxu = objdimx.getncattr('units') odimyu = objdimy.getncattr('units') if len(objdimx.shape) <= 2: odimxv = objdimx[:] odimyv = objdimy[:] elif len(objdimx.shape) == 3: odimxv = objdimx[0,:] odimyv = objdimy[0,:] else: print errormsg print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \ ' not ready!!' quit(-1) timename = timevals.split('|')[0] timeunit = timevals.split('|')[1].replace('!',' ') timekind = timevals.split('|')[2] timefmt = timevals.split('|')[3] timelabel = timevals.split('|')[4].replace('!',' ') if vdimxn == timename: odimxv = objsf.variables[vdimxn][:] odimxu = timelabel timeaxis = 'x' odimyv = objsf.variables[vdimyn] odimyu = odimyv.getncattr('units') timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt) elif vdimyn == timename: odimyv = objsf.variables[vdimyn][:] odimyu = timelabel timeaxis = 'y' odimxv = objsf.variables[vdimxn] odimxu = odimxv.getncattr('units') timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt) else: print errormsg print ' ' + fname + ": time variable '" + timename + "' not found!!" quit(-1) shading_nx = np.zeros((2), dtype=np.float) shading_nx[0] = np.float(shadminmax.split(',')[0]) shading_nx[1] = np.float(shadminmax.split(',')[1]) closeval = drw.Str_Bool(close) drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \ dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \ timepos, timelabels, False) # Line values ## ncfilel = ncfile.split(',')[1] vnamelfig = values.split(':')[0].split(',')[1] varnl = varn.split(',')[1] objlf = NetCDFFile(ncfilel,'r') objlvar = objlf.variables[varnl] linevals = objlvar[:] if reva0 == 'tranpose': plt.plot (linevals, odimxv, '-', color='k') else: plt.plot (odimxv, linevals, '-', color='k') objsf.close() objsl.close() return def draw_topo_geogrid(ncfile, values): """ plotting geo_em.d[nn].nc topography from WPS files draw_topo_geogrid(ncfile, values) ncfile= geo_em.d[nn].nc file to use values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues] [min/max]Topo: minimum and maximum values of topography to draw lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None title: title of the graph graphic_kind: kind of figure (jpg, pdf, png) mapvalues: map characteristics [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lambert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full """ fname = 'draw_topo_geogrid' if values == 'h': print fname + '_____________________________________________________________' print draw_topo_geogrid.__doc__ quit() expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \ '[mapvalues]'] drw.check_arguments(fname,5,values,':',expectargs) mintopo = values.split(':')[0].split(',')[0] maxtopo = values.split(':')[0].split(',')[1] lonlatLS = values.split(':')[1] lonlatLv = lonlatLS.split(',')[0] if lonlatLv == 'None': lonlatL = None else: lonlatL = np.zeros((4), dtype=np.float) lonlatL[0] = np.float(lonlatLS.split(',')[0]) lonlatL[1] = np.float(lonlatLS.split(',')[1]) lonlatL[2] = np.float(lonlatLS.split(',')[2]) lonlatL[3] = np.float(lonlatLS.split(',')[3]) grtit = values.split(':')[2] kindfig = values.split(':')[3] mapvalues = values.split(':')[4] if not os.path.isfile(ncfile): print errormsg print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!' quit(-1) objdomf = NetCDFFile(ncfile, 'r') objhgt = objdomf.variables['HGT_M'] objlon = objdomf.variables['XLONG_M'] objlat = objdomf.variables['XLAT_M'] topography = objhgt[0,:,:] drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL, \ grtit, kindfig, mapvalues, True) objdomf.close() return def draw_topo_geogrid_boxes(ncfiles, values): """ plotting different geo_em.d[nn].nc topography from WPS files draw_topo_geogrid_boxes(ncfile, values) ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference) values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels] [min/max]Topo: minimum and maximum values of topography to draw lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None title: title of the graph graphic_kind: kind of figure (jpg, pdf, png) mapvalues: map characteristics [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric * 'lcc', lambert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full labels= labels to write in the graph """ # import matplotlib as mpl # mpl.use('Agg') import matplotlib.pyplot as plt fname = 'draw_topo_geogrid_boxes' if values == 'h': print fname + '_____________________________________________________________' print draw_topo_geogrid_boxes.__doc__ quit() mintopo = values.split(':')[0].split(',')[0] maxtopo = values.split(':')[0].split(',')[1] lonlatLS = values.split(':')[1] lonlatLv = lonlatLS.split(',')[0] if lonlatLv == 'None': lonlatL = None else: lonlatL = np.zeros((4), dtype=np.float) lonlatL[0] = np.float(lonlatLS.split(',')[0]) lonlatL[1] = np.float(lonlatLS.split(',')[1]) lonlatL[2] = np.float(lonlatLS.split(',')[2]) lonlatL[3] = np.float(lonlatLS.split(',')[3]) grtit = values.split(':')[2] kindfig = values.split(':')[3] mapvalues = values.split(':')[4] labels = values.split(':')[5] ncfile = ncfiles.split(',')[0] if not os.path.isfile(ncfile): print errormsg print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!' quit(-1) objdomf = NetCDFFile(ncfile, 'r') objhgt = objdomf.variables['HGT_M'] objlon0 = objdomf.variables['XLONG_M'] objlat0 = objdomf.variables['XLAT_M'] topography = objhgt[0,:,:] Nfiles = len(ncfiles.split(',')) boxlabels = labels.split(',') Xboxlines = [] Yboxlines = [] for ifile in range(Nfiles): ncfile = ncfiles.split(',')[ifile] # print ifile, ncfile if not os.path.isfile(ncfile): print errormsg print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!' quit(-1) objdomfi = NetCDFFile(ncfile, 'r') objlon = objdomfi.variables['XLONG_M'] objlat = objdomfi.variables['XLAT_M'] dx = objlon.shape[2] dy = objlon.shape[1] Xboxlines.append(objlon[0,0,:]) Yboxlines.append(objlat[0,0,:]) Xboxlines.append(objlon[0,dy-1,:]) Yboxlines.append(objlat[0,dy-1,:]) Xboxlines.append(objlon[0,:,0]) Yboxlines.append(objlat[0,:,0]) Xboxlines.append(objlon[0,:,dx-1]) Yboxlines.append(objlat[0,:,dx-1]) objdomfi.close() drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels, \ objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True) objdomf.close() return def movievalslice(origslice, dimmovien, framenum): """ Function to provide variable slice according to a geneation of a movie movievals(origslice, dimmovien, framenum) [origslice]= slice original as [dimname1]|[val1],[...,[dimnameN]|[valN]] ([val] = -1, full length) [dimmovien]= name of the dimension to produce the movie [framenum]= value of the frame to substitue in [origslice] as [dimmovien]|[framenum] >>> movievalslice('East_West|-1,North_South|-1,Time|2','Time',0) East_West|-1,North_South|-1,Time|0 """ fname = 'movievalslice' if origslice == 'h': print fname + '_____________________________________________________________' print movievalslice.__doc__ quit() dims = origslice.split(',') movieslice = '' idim = 0 for dimn in dims: dn = dimn.split('|')[0] if dn == dimmovien: movieslice = movieslice + dn + '|' + str(framenum) else: movieslice = movieslice + dimn if idim < len(dims)-1: movieslice = movieslice + ',' idim = idim + 1 return movieslice class Capturing(list): """ Class to capture function output as a list from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call """ # from cStringIO import StringIO def __enter__(self): self._stdout = sys.stdout sys.stdout = self._stringio = StringIO() return self def __exit__(self, *args): self.extend(self._stringio.getvalue().splitlines()) sys.stdout = self._stdout def create_movie(netcdfile, values, variable): """ Function to create a movie assuming ImageMagick installed! values= [graph]#[movie_dimension]#[graph_values] [graph]: which graphic [movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval] [dimnmovie]; name of the dimension from which make the movie [dimvmovie]; name of the variable with the values of the dimension [moviedelay]; delay between frames [interval]; [beg]@[end]@[freq] or -1 (all) [graph_values]: values to generate the graphic netcdfile= netCDF file variable= variable to use (when applicable) """ fname = 'create_movie' if values == 'h': print fname + '_____________________________________________________________' print create_movie.__doc__ quit() graph = values.split('#')[0] movie_dim = values.split('#')[1] graph_vals = values.split('#')[2] ncobj = NetCDFFile(netcdfile, 'r') # Movie dimension ## dimnmovie = movie_dim.split('@')[0] dimvmovie = movie_dim.split('@')[1] moviedelay = movie_dim.split('@')[2] moviebeg = int(movie_dim.split('@')[3]) if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie): print errormsg print ' ' + fname + ": file '" + netcdfile + "' has not dimension '" + \ dimnmovie + "' !!!" quit(-1) objdmovie = ncobj.dimensions[dimnmovie] dmovie = len(objdmovie) if moviebeg != -1: moviend = int(movie_dim.split('@')[4]) moviefreq = int(movie_dim.split('@')[5]) else: moviebeg = 0 moviend = dmovie moviefreq = 1 if dimvmovie == 'WRFTimes': objvdmovie = ncobj.variables['Times'] vdmovieunits = '' valsdmovie = [] for it in range(objvdmovie.shape[0]): valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:], \ 'WRFdatetime', 'Y/m/d H-M-S')) elif dimvmovie == 'CFtime': objvdmovie = ncobj.variables['time'] vdmovieunits = '' print objvdmovie.units valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard', \ objvdmovie[:]) valsdmovie = [] for it in range(objvdmovie.shape[0]): valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:], \ 'matYmdHMS', 'Y/m/d H-M-S')) else: if not drw.searchInlist(ncobj.variables.keys(),dimvmovie): print errormsg print ' ' + fname + ": file '" + netcdfile + "' has not variable '" + \ dimvmovie + "' !!!" quit(-1) vdmovieunits = objvdmovie.getncattr('units') objvdmovie = ncobj.variables[dimvmovie] if len(objvdmovie.shape) == 1: vasldmovie = objvdmovie[:] else: print errormsg print ' ' + fname + ': shape', objvdmovie.shape, 'of variable with ' + \ 'dimension movie values not ready!!!' quit(-1) ncobj.close() os.system('rm frame_*.png > /dev/null') # graphic ## if graph == 'draw_2D_shad': graphvals = graph_vals.split(':') for iframe in range(moviebeg,moviend,moviefreq): iframeS = str(iframe).zfill(4) drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \ 5, 'frames') titgraph = dimnmovie + '|=|' + str(valsdmovie[iframe]) + '|' + \ vdmovieunits graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe) graphvals[6] = titgraph graphvals[7] = 'png' graphv = drw.numVector_String(graphvals, ":") with Capturing() as output: draw_2D_shad(netcdfile, graphv, variable) os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png') else: print errormsg print ' ' + fname + ": graphic '" + graph + "' not defined !!!" quit(-1) os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif') print "Succesfuly creation of movie file 'create_movie.gif' !!!" return def draw_lines(ncfilens, values, varname): """ Function to draw different lines at the same time from different files draw_lines(ncfilens, values, varname): ncfilens= [filen] ',' separated list of netCDF files values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk] [dimvname]: name of the variable with he values of the common dimension [valuesaxis]: which axis will be used for the values ('x', or 'y') [dimtit]: title for the common dimension [leglabels]: ',' separated list of names for the legend [vartit]: name of the variable in the graph [title]: title of the plot ('|' for spaces) [locleg]: location of the legend (-1, autmoatic) 1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right', 5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center', 9: 'upper center', 10: 'center' [graphk]: kind of the graphic varname= variable to plot values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf' """ fname = 'draw_lines' if values == 'h': print fname + '_____________________________________________________________' print draw_lines.__doc__ quit() expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]' drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs) ncfiles = ncfilens.split(',') dimvname = values.split(':')[0] valuesaxis = values.split(':')[1] dimtit = values.split(':')[2] leglabels = values.split(':')[3].replace('_','\_') vartit = values.split(':')[4] title = values.split(':')[5].replace('|',' ') graphk = values.split(':')[6] Nfiles = len(ncfiles) # Getting trajectotries ## varvalues = [] dimvalues = [] print ' ' + fname ifn = 0 for ifile in ncfiles: filen = ifile.split('@')[0] print ' filen:',filen if not os.path.isfile(filen): print errormsg print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!" quit(-1) objfile = NetCDFFile(filen, 'r') if not objfile.variables.has_key(dimvname): print errormsg print ' ' + fname + ": netCDF file '" + filen + \ "' does not have variable '" + dimvname + "' !!" quit(-1) if not objfile.variables.has_key(varname): print errormsg print ' ' + fname + ": netCDF file '" + filen + \ "' does not have variable '" + varname + "' !!" quit(-1) vvobj = objfile.variables[varname] if len(vvobj.shape) != 1: print errormsg print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \ varname + "' !!" quit(-1) vdobj = objfile.variables[dimvname] if len(vdobj.shape) != 1: print errormsg print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \ dimvname + "' !!" quit(-1) varvalues.append(vvobj[:]) dimvalues.append(vdobj[:]) if ifn == 0: varunits = vvobj.units objfile.close() ifn = ifn + 1 drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \ vartit, varunits, title, graphk) return def draw_lines_time(ncfilens, values, varname): """ Function to draw different lines at the same time from different files with times draw_lines_time(ncfilens, values, varname): ncfilens= [filen] ',' separated list of netCDF files values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[timevals];[locleg];[graphk] [dimvname]: name of the variable with he values of the common dimension [valuesaxis]: which axis will be used for the values ('x', or 'y') [dimtit]: title for the common dimension [leglabels]: ',' separated list of names for the legend [vartit]: name of the variable in the graph [title]: title of the plot ('|' for spaces) [timevals]: [timen]|[units]|[kind]|[tfmt] time labels characteristics [timen]; name of the time variable [units]; units string according to CF conventions ([tunits] since [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces) [kind]; kind of output 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [tfmt]; desired format [locleg]: location of the legend (-1, autmoatic) 1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right', 5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center', 9: 'upper center', 10: 'center' [graphk]: kind of the graphic varname= variable to plot values= 'time;y;time ([DD]${[HH]}$);32x32;$wss^{*}$;wss Taylor's turbulence term;time|hours!since!1949-12-01_00:00:00;exct,12,h|%d$^{%H}$;2;pdf' """ fname = 'draw_lines_time' if values == 'h': print fname + '_____________________________________________________________' print draw_lines.__doc__ quit() expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];' expectargs = expectargs + '[timevals];[locleg];[graphk]' drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs) ncfiles = ncfilens.split(',') dimvname = values.split(';')[0] valuesaxis = values.split(';')[1] dimtit = values.split(';')[2] leglabels = values.split(';')[3].replace('_','\_') vartit = values.split(';')[4] title = values.split(';')[5].replace('|',' ') timevals = values.split(';')[6] locleg = int(values.split(';')[7]) graphk = values.split(';')[8] Nfiles = len(ncfiles) # Getting values ## varvalues = [] dimvalues = [] timvalues = [] timvals0 = timvalues print ' ' + fname ifn = 0 mintval = 1.e20 maxtval = -1.e20 for ifile in ncfiles: filen = ifile.split('@')[0] print ' filen:',filen if not os.path.isfile(filen): print errormsg print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!" quit(-1) objfile = NetCDFFile(filen, 'r') if not objfile.variables.has_key(dimvname): print errormsg print ' ' + fname + ": netCDF file '" + filen + \ "' does not have variable '" + dimvname + "' !!" quit(-1) if not objfile.variables.has_key(varname): print errormsg print ' ' + fname + ": netCDF file '" + filen + \ "' does not have variable '" + varname + "' !!" quit(-1) vvobj = objfile.variables[varname] if len(vvobj.shape) != 1: print errormsg print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \ varname + "' !!" quit(-1) vdobj = objfile.variables[dimvname] if len(vdobj.shape) != 1: print errormsg print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \ dimvname + "' !!" quit(-1) varvalues.append(vvobj[:]) dimvalues.append(vdobj[:]) mindvals = np.min(vdobj[:]) maxdvals = np.max(vdobj[:]) if mindvals < mintval: mintval = mindvals if maxdvals > maxtval: maxtval = maxdvals if ifn == 0: varunits = vvobj.units objfile.close() ifn = ifn + 1 # Times timename = timevals.split('|')[0] timeunit = timevals.split('|')[1].replace('!',' ') timekind = timevals.split('|')[2] timefmt = timevals.split('|')[3] dtvals = (maxtval - mintval)/5. tvals = np.arange(mintval, maxtval, dtvals/2.) timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt) drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \ vartit, varunits, timepos, timelabels, title, locleg, graphk) return def draw_Neighbourghood_evol(filen, values, variable): """ Function to draw the temporal evolution of a neighbourghood around a point draw_Neighbourghood_evol(filen, values, variable) filen= netCDF file name values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]: [timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile] [dimsval]: [dimn1]|[val1]|[dimv1],...,[dimnN]|[valN]|[dimvN] dimension names, values to get (-1, for all; no name/value pair given full length) and variable with values of the dimension NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2 [neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined [Nneig]: Number of grid points of the full side of the box (odd value) [Ncol]: Number of columns ('auto': square final plot) [gvarname]: name of the variable to appear in the graph [timetits]: [titX],[titY] titles of the axes ('|' for spaces) [tkinds]: [tkindX]|[tkindY] kinds of time to appear in the graph 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [timefmts]: [tfmtX],[tfmtY] format of the time labels [gtitle]: title of the graphic ('|' for spaces) [shadxtrms]: Extremes for the shading [cbar]: colorbar to use [gkind]: kind of graphical output [ofile]: True/False whether the netcdf with data should be created or not variable= name of the variable values = 'q:Time|-1|Times,bottom_top|6|ZNU,south_north|3|XLAT,west_east|26|XLONG:south_north,west_east:5:auto:time|($[DD]^{[HH]}$),time|($[DD]^{[HH]}$):exct,2,h|exct,1,d:$%d^{%H}$,$%d^{%H}$:5|pts|neighbourghood|temporal|evolution:0.0,0.004:BuPu:pdf:True' """ fname = 'draw_Neighbourghood_evol' if values == 'h': print fname + '_____________________________________________________________' print draw_Neighbourghood_evol.__doc__ quit() expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]', \ '[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]', \ '[gkind]', '[ofile]'] drw.check_arguments(fname,len(expectargs),values,':',expectargs) gvarname = values.split(':')[0] dimsval = values.split(':')[1].split(',') neigdims = values.split(':')[2].split(',') Nneig = int(values.split(':')[3]) Ncol0 = values.split(':')[4] timetits = values.split(':')[5].split(',') timekinds = values.split(':')[6].split('|') timefmts = values.split(':')[7].split(',') gtitle = values.split(':')[8].replace('|',' ') shadxtrms = values.split(':')[9].split(',') cbar = values.split(':')[10] gkind = values.split(':')[11] ofile = values.split(':')[12] if Ncol0 != 'auto': Ncol = int(Ncol0) else: Ncol = Ncol0 timetits[0] = timetits[0].replace('|',' ') timetits[1] = timetits[1].replace('|',' ') if np.mod(Nneig,2) == 0: print errormsg print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!' quit(-1) Nneig2 = int(Nneig/2) # Values to slice the variable dimvslice = {} dimvvalues = {} for dimvs in dimsval: dimn = dimvs.split('|')[0] dimv = int(dimvs.split('|')[1]) dimnv = dimvs.split('|')[2] dimvvalues[dimn] = dimnv dimvslice[dimn] = dimv ncobj = NetCDFFile(filen, 'r') varobj = ncobj.variables[variable] slicevar = [] newdimn = [] newdimsvar = {} for dimn in varobj.dimensions: if not drw.searchInlist(dimvslice.keys(), dimn): dimsize = len(ncobj.dimensions[dimn]) slicevar.append(slice(0, dimsize+1)) newdimn.append(dimn) newdimsvar[dimn] = dimseize for dimslicen in dimvslice.keys(): if dimn == dimslicen: if dimvslice[dimn] != -1: if drw.searchInlist(neigdims, dimn): slicevar.append(slice(dimvslice[dimn]-Nneig2, \ dimvslice[dimn]+Nneig2+1)) newdimn.append(dimn) newdimsvar[dimn] = Nneig break else: slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1)) break else: dimsize = len(ncobj.dimensions[dimn]) slicevar.append(slice(0, dimsize+1)) newdimn.append(dimn) newdimsvar[dimn] = dimsize break varv = varobj[tuple(slicevar)] if len(newdimn) != 3: print errormsg print ' ' + fname + ': sliced variable with shape=', varv.shape, \ ' must have three dimensions',len(varv.shape),'given !!' quit(-1) newdims = [] for nwdims in newdimn: newdims.append(newdimsvar[nwdims]) # The dimension which is not in the neighbourhood dimensions must be time! for dim1 in newdimn: if not drw.searchInlist(neigdims, dim1): dimt = newdimsvar[dim1] dimtime = dim1 if Ncol == 'auto': dimtsqx = int(np.sqrt(dimt)) + 1 dimtsqy = int(np.sqrt(dimt)) + 1 else: dimtsqx = int(Ncol) dimtsqy = dimt/dimtsqx + 1 neighbourghood = np.ones((dimtsqy*Nneig,dimtsqx*Nneig), dtype=np.float)*fillValue for it in range(dimt): ity = int(it/dimtsqx) itx = it-ity*dimtsqx itty = (dimtsqy - ity - 1)*Nneig + Nneig2 ittx = itx*Nneig + Nneig2 neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \ varv[it,::-1,:] variablevals = drw.variables_values(variable) if drw.searchInlist(varobj.ncattrs(), 'units'): vunits = varobj.units else: vunits = variablevals[5] # Time values at the X/Y axes if ncobj.variables[dimvvalues[dimtime]].dtype == '|S1': print ' ' + fname + ': WRF time variable!' refdate = '19491201000000' tunitsval = 'hours' dimtvalues = np.zeros((dimt), dtype=np.float) tvals = ncobj.variables[dimvvalues[dimtime]] yrref=refdate[0:4] monref=refdate[4:6] dayref=refdate[6:8] horref=refdate[8:10] minref=refdate[10:12] secref=refdate[12:14] refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \ minref + ':' + secref tunits = tunitsval + ' since ' + refdateS for it in range(dimt): wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS') dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval) else: dimtvalues = ncobj.variables[dimvvalues[dimtime]][:] tunits = ncobj.variables[newdimsvar[dimtime]].units dimxv = dimtvalues[0:dimtsqx] dimyv = dimtvalues[0:dimt:dimtsqx] print 'Lluis lens: dimxv, dimyv:',len(dimxv), len(dimyv) dimn = ['time','time'] if ofile == 'True': ofilen = 'Neighbourghood_evol.nc' newnc = NetCDFFile(ofilen, 'w') # Dimensions newdim = newnc.createDimension('time',None) newdim = newnc.createDimension('y',dimtsqy*Nneig) newdim = newnc.createDimension('x',dimtsqx*Nneig) # Dimension values newvar = newnc.createVariable('time','f8',('time')) newvar[:] = np.arange(dimt) newattr = drw.basicvardef(newvar, 'time','time',tunits) # Neighbourhghood variable newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \ fill_value=fillValue) newvar[:] = neighbourghood newnc.sync() newnc.close() print fname + ": Successfull generation of file '" + ofilen + "' !!" # Time ticks timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0]) timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1]) timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]] timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]] for i in range(2): if shadxtrms[i][0:1] != 'S': shadxtrms[i] = np.float(shadxtrms[i]) drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \ timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True) def draw_timeSeries(filen, values, variables): """ Function to draw a time-series draw_timeSeries(filen, values, variable): filen= name of the file values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind] [gvarname]: name of the variable to appear in the graph [timetit]: title of the time axis (assumed x-axis, '|' for spaces) [tkind]: kind of time to appear in the graph (assumed x-axis) 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [timefmt]: format of the time labels [title]: title of the graphic ('|' for spaces) [locleg]: location of the legend (-1, autmoatic) 1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right', 5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center', 9: 'upper center', 10: 'center' [gkind]: kind of graphical output variables= [varname],[timename] names of variable and variable with times draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf', 'LDQCON,time') """ fname = 'draw_timeSeries' if values == 'h': print fname + '_____________________________________________________________' print draw_timeSeries.__doc__ quit() expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]', \ '[locleg]', '[gkind]'] drw.check_arguments(fname,len(expectargs),values,':',expectargs) gvarname = values.split(':')[0] timetit = values.split(':')[1].replace('|',' ') tkind = values.split(':')[2] timefmt = values.split(':')[3] title = values.split(':')[4].replace('|',' ') locleg = int(values.split(':')[5]) gkind = values.split(':')[6] ncobj = NetCDFFile(filen, 'r') variable = variables.split(',')[0] timevar = variables.split(',')[1] if not ncobj.variables.has_key(variable): print errormsg print ' ' + fname + ": file '" + filen + "' does not have variable '" + \ variable + "' !!" quit(-1) if not ncobj.variables.has_key(timevar): print errormsg print ' ' + fname + ": file '" + filen + "' does not have variable time '" \ + timevar + "' !!" quit(-1) varobj = ncobj.variables[variable] timeobj = ncobj.variables[timevar] dimt = len(timeobj[:]) varvals = np.zeros((2,dimt), dtype=np.float) gunits = varobj.getncattr('units') tunits = timeobj.getncattr('units') varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits) varvals[1,:] = timeobj[:] tseriesvals = [] tseriesvals.append(varvals) drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \ 'TimeSeries', gvarname, timetit, tkind, timefmt, title, \ gvarname.replace('_','\_'), locleg, gkind) return #draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf', 'LDQCON,time') def draw_trajectories(trjfilens, values, observations): """ Function to draw different trajectories at the same time draw_trajectories(trjfilens, values, observations): trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories, time intervals and reference maps (first one will be used to plot) [filen]: name of the file to use (lines with '#', not readed) as: [t-step] [x] [y] [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval [map]: [file]#[lonname]#[latname] [file]; with the [lon],[lat] matrices [lonname],[latname]; names of the longitudes and latitudes variables values=[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind] [leglabels]: ',' separated list of names for the legend [lonlatlims]: limits of the map [lonmin, latmin, lonmax, latmax] or None [title]: title of the plot [graphk]: kind of the graphic [mapkind]: drawing coastaline ([proj],[res]) or None [proj]: projection * 'cyl', cilindric * 'lcc', lambert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full obsevations= [obsfile],[obsname],[Tint],[null] [obsfile]: name fo the File with the observations as [t-step] [lat] [lon] [obsname]: name of the observations in the graph [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval [null]: null value for the observed trajectory """ fname = 'draw_trajectories' if values == 'h': print fname + '_____________________________________________________________' print draw_trajectories.__doc__ quit() trjfiles = trjfilens.split(',') leglabels = values.split('|')[0] lonlatlims = values.split('|')[1] title = values.split('|')[2] graphk = values.split('|')[3] mapkind = values.split('|')[4] Nfiles = len(trjfiles) # Getting trajectotries ## lontrjvalues = [] lattrjvalues = [] print ' ' + fname ifn = 0 for ifile in trjfiles: filen = ifile.split('@')[0] Tint = ifile.split('@')[1] print ' trajectory:',filen if Tint != '-1': Tbeg = Tint Tend = ifile.split('@')[2] mapv = ifile.split('@')[3] else: mapv = ifile.split('@')[2] if not os.path.isfile(filen): print errormsg print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!" quit(-1) # Charging longitude and latitude values ## lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \ mapv.split('#')[2]) if ifn == 0: mapref = mapv ifn = ifn + 1 objfile = open(filen, 'r') trjtimev = [] trjxv = [] trjyv = [] for line in objfile: if line[0:1] != '#': trjtimev.append(int(line.split(' ')[0])) trjxv.append(int(line.split(' ')[1])) trjyv.append(int(line.split(' ')[2])) objfile.close() if Tint != '-1': lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]]) lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]]) else: lontrjvalues.append(lonvals[trjyv[:], trjxv[:]]) lattrjvalues.append(latvals[trjyv[:], trjxv[:]]) # lonlatlimits ## if lonlatlims == 'None': lonlatlimsv = None else: lonlatlimsv = np.zeros((4), dtype=np.float) lonlatlimsv[0] = np.float(lonlatlims.split(',')[0]) lonlatlimsv[1] = np.float(lonlatlims.split(',')[1]) lonlatlimsv[2] = np.float(lonlatlims.split(',')[2]) lonlatlimsv[3] = np.float(lonlatlims.split(',')[3]) # lon/lat objects ## objnc = NetCDFFile(mapref.split('#')[0]) lonobj = objnc.variables[mapref.split('#')[1]] latobj = objnc.variables[mapref.split('#')[2]] # map ## if mapkind == 'None': mapkindv = None else: mapkindv = mapkind if observations is None: obsname = None else: obsfile = observations.split(',')[0] obsname = observations.split(',')[1] Tint = observations.split(',')[2] null = np.float(observations.split(',')[3]) print ' observational trajectory:',obsfile if not os.path.isfile(obsfile): print errormsg print ' ' + fname + ": observations trajectory file '" + obsfile + \ "' does not exist !!" quit(-1) objfile = open(obsfile, 'r') obstrjtimev = [] obstrjxv = [] obstrjyv = [] for line in objfile: if line[0:1] != '#': lon = np.float(line.split(' ')[2]) lat = np.float(line.split(' ')[1]) if not lon == null and not lat == null: obstrjtimev.append(int(line.split(' ')[0])) obstrjxv.append(lon) obstrjyv.append(lat) else: obstrjtimev.append(int(line.split(' ')[0])) obstrjxv.append(None) obstrjyv.append(None) objfile.close() if Tint != '-1': Tint = int(observations.split(',')[2].split('@')[0]) Tend = int(observations.split(',')[2].split('@')[1]) lontrjvalues.append(obstrjxv[Tint:Tend+1]) lattrjvalues.append(obstrjyv[Tint:Tend+1]) else: lontrjvalues.append(obstrjxv[:]) lattrjvalues.append(obstrjyv[:]) drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \ lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname) objnc.close() return def draw_vals_trajectories(ncfile, values, variable): """ Function to draw values from the outputs from 'compute_tevolboxtraj' draw_vals_trajectories(ncfile, values, variable) ncfile= [ncfile] ',' list of files to use values= [statistics]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind] [statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean', 'mean2', 'stdev' [Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times [labels]: ',' separated list of labels for the legend [locleg]: location of the legend (-1, autmoatic) 1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right', 5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center', 9: 'upper center', 10: 'center' [gvarname]: name of the variable to appear in the graph [timetit]: title of the time axis (assumed x-axis, '|' for spaces) [tkind]: kind of time to appear in the graph (assumed x-axis) 'Nval': according to a given number of values as 'Nval',[Nval] 'exct': according to an exact time unit as 'exct',[tunit]; tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month, 'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 'l': milisecond [timefmt]: format of the time labels [title]: title of the graphic ('|' for spaces) [gkind]: kind of graphical output variable= variable to use """ fname = 'draw_vals_trajectories' if values == 'h': print fname + '_____________________________________________________________' print draw_vals_trajectories.__doc__ quit() sims = ncfile.split(',') if len(values.split(':')) != 9: print errormsg print ' ' + fname + ': wrong number of values!', len(values.split(':')), \ 'given 9 needed!!' print ' ',values.split(':') quit(-1) statistics = values.split(':')[0] Tint = values.split(':')[1] labels = values.split(':')[2] gvarname = values.split(':')[3] timetit = values.split(':')[4].replace('|',' ') tkind = values.split(':')[5] timefmt = values.split(':')[6] title = values.split(':')[7].replace('|',' ') gkind = values.split(':')[8] leglabels = labels.split('@')[0].split(',') locleg = int(labels.split('@')[1]) Nsims = len(sims) if Tint != '-1': tini = np.float(Tint.split('@')[0]) tend = np.float(Tint.split('@')[1]) else: tini = -1. tend = -1. vartimetrjv = [] print ' ' + fname for trjfile in sims: print ' ' + trjfile + ' ...' if not os.path.isfile(trjfile): print errormsg print ' ' + fname + ": trajectory file: '" + trjfile + \ "' does not exist !!" quit(-1) trjobj = NetCDFFile(trjfile, 'r') otim = trjobj.variables['time'] if not trjobj.variables.has_key(statistics + 'box_' + variable): print errormsg print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\ statistics + 'box_' + variable + "' !!" quit(-1) ovar = trjobj.variables[statistics + 'box_' + variable] dimt = otim.shape[0] if trjfile == sims[0]: gunits = ovar.getncattr('units') lname = ovar.getncattr('long_name') tunits = otim.getncattr('units') if tini != -1: tiniid = -1 tendid = -1 for itv in range(dimt): if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv if tiniid == -1 or tendid == -1: print errormsg print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \ tendid, ',', tiniid, ' !!' print ' data interval [',otim[0], otim[dimt-1],']' quit(-1) dimt = tendid - tiniid + 1 else: dimt = otim.shape[0] valsv = np.zeros((2,dimt), dtype=np.float) if tini == -1: valsv[1,:] = otim[:] valsv[0,:] = ovar[:] else: valsv[1,:] = otim[tiniid:tendid+1] valsv[0,:] = ovar[tiniid:tendid+1] vartimetrjv.append(valsv) trjobj.close() drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \ 'val_trajectories_' + statistics, gvarname, timetit, tkind, timefmt, title, \ leglabels, locleg, gkind) def variable_values(values): """ Function to give back values for a given variable values= [varname] name of the variable """ fname = 'variable_values' values = drw.variables_values(values) print fname,'values:',values print fname,'variable_name:',values[0] print fname,'standard_name:',values[1] print fname,'min,max:',str(values[2]) + ',' + str(values[3]) print fname,'long_name:',values[4] print fname,'units:',values[5] print fname,'shad_colors:',values[6] print fname,'all_values:',drw.numVector_String(values,',') return ####### ###### ##### #### ### ## # ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'" ### Options ##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]" string_operation="""operation to make: draw_topo_geogrid, draws topography from a WPS geo_em.d[nn].nc: -S [minTopo],[maxTopo]:[SW_lon],[SW_lat],[NE_lon],[NE_lat]:[title]:[graphic_kind]:[projection],[res_coastline] draw_2D_shad_cont, draws two 2D fields, first with shading second with contour lines: -v [varns],[varnc] -S [vnamefs],[vnamefc],[dimxvn],[dimyvn],[colorbar],[ckind],[clabfmt],[sminv]:[smaxv],[sminc]:[smaxv]:[Nlev],[figt],[kindfig],[reverse] [ckind]: 'cmap': as it gets from colorbar 'fixc,[colname]': fixed color [colname], all stright lines 'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line """ #print string_operation parser = OptionParser() parser.add_option("-f", "--netCDF_file", dest="ncfile", help="file to use", metavar="FILE") parser.add_option("-o", "--operation", type='choice', dest="operation", choices=namegraphics, help="operation to make: " + ngraphics, metavar="OPER") parser.add_option("-S", "--valueS", dest="values", help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES") parser.add_option("-v", "--variable", dest="varname", help="[WHEN APPLICABLE] variable to check", metavar="VAR") (opts, args) = parser.parse_args() ####### ####### ## MAIN ####### # Not checking file operation Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \ 'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \ 'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \ 'variable_values'] ####### ###### ##### #### ### ## # errormsg='ERROR -- error -- ERROR -- error' varn=opts.varname oper=opts.operation if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \ not drw.searchInlist(Notcheckingfile, oper): print errormsg print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!' quit(-1) if oper == 'create_movie': create_movie(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad': draw_2D_shad(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad_time': draw_2D_shad_time(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad_cont': draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad_cont_time': draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad_line': draw_2D_shad_line(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_2D_shad_line_time': draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_Neighbourghood_evol': draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_lines': draw_lines(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_lines_time': draw_lines_time(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_timeSeries': draw_timeSeries(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_topo_geogrid': draw_topo_geogrid(opts.ncfile, opts.values) elif oper == 'draw_topo_geogrid_boxes': draw_topo_geogrid_boxes(opts.ncfile, opts.values) elif oper == 'draw_trajectories': draw_trajectories(opts.ncfile, opts.values, opts.varname) elif oper == 'draw_vals_trajectories': draw_vals_trajectories(opts.ncfile, opts.values, opts.varname) elif oper == 'list_graphics': # From: http://www.diveintopython.net/power_of_introspection/all_together.html import drawing as myself object = myself for opern in namegraphics: if opern != 'list_graphics': print opern + '_______ ______ _____ ____ ___ __ _' print getattr(object, opern).__doc__ elif oper == 'variable_values': variable_values(opts.values) else: print errormsg print ' ' + main + ": the graphic '" + oper + "' is not ready !!" print errormsg quit()