# -*- coding: iso-8859-15 -*- # Python tools to plot netCDF files. # From L. Fita work in different places: LMD (France) # More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot # # pyNCplot and its component drawing_tools.py comes with ABSOLUTELY NO WARRANTY. # This work is licendes under a Creative Commons # Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0) # #import pylab as plt # From http://stackoverflow.com/questions/13336823/matplotlib-python-error import numpy as np import matplotlib as mpl mpl.use('Agg') from matplotlib.pylab import * import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap import os from netCDF4 import Dataset as NetCDFFile import nc_var_tools as ncvar import generic_tools as gen errormsg = 'ERROR -- error -- ERROR -- error' warnmsg = 'WARNING -- waring -- WARNING -- warning' fillValue = 1.e20 fillValueF = 1.e20 colorsauto = ['#FF0000', '#00FF00', '#0000FF', '#FF00FF', '#00FFFF', '#FFAA00', \ '#AA0000', '#00AA00', '#0000AA', '#AA00AA', '#00AAAA', '#AA3200'] pointkindsauto = ['.', 'x', '+', '*', '|', '_', '1', '2', '3', '4', 'o', '<', '>', \ 'v', '^', 's', 'D', 'p' ,'h' ,'H'] linekindsauto = ['-', '--', '-.', ':'] linewidthsauto = [1.] pointsizesauto = [7.] ####### Funtions # check_colorBar: Check if the given colorbar exists in matplotlib # colorbar_vals: Function to provide the colorbar values for a figure # format_axes: Function to provide the format for the ticks of the axes in a figure # graphic_range: Function to provide the ranges of a figure following different options # legend_values: Function to determine the values for the legend # searchInlist: # datetimeStr_datetime: # dateStr_date: # numVector_String: # timeref_datetime: # slice_variable: Function to return a slice of a given variable according to values to its dimension # interpolate_locs: # datetimeStr_conversion: # percendone: # netCDFdatetime_realdatetime: # file_nlines: # variables_values: # units_lunits: # ASCII_LaTeX: # DegGradSec_deg: # intT2dt: # lonlat2D: Function to return lon, lat 2D matrices from any lon,lat matrix # lonlat_values: # date_CFtime: # pot_values: # CFtimes_plot: # color_lines: # output_kind: # check_arguments: # Str_Bool: # plot_points: # plot_2Dfield: # plot_2Dfield_easy: # plot_topo_geogrid: # plot_topo_geogrid_boxes: # plot_2D_shadow: # plot_2D_shadow_time: Plotting a 2D field with one of the axes being time # plot_Neighbourghood_evol:Plotting neighbourghood evolution# plot_Trajectories # plot_2D_shadow_contour: # plot_2D_shadow_contour_time: # dxdy_lonlat: Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values # plot_2D_shadow_line: # plot_2lines: Function to plot two lines in different axes (x/x2 or y/y2) # plot_2lines_time: Function to plot two time-lines in different axes (x/x2 or y/y2) # plot_lines: Function to plot a collection of lines # plot_ZQradii: Function to plot following radial averages only at exact grid poins # transform: Function to transform the values and the axes # From nc_var_tools.py def reduce_spaces(string): """ Function to give words of a line of text removing any extra space """ values = string.replace('\n','').split(' ') vals = [] for val in values: if len(val) > 0: vals.append(val) return vals def searchInlist(listname, nameFind): """ Function to search a value within a list listname = list nameFind = value to find >>> searInlist(['1', '2', '3', '5'], '5') True """ for x in listname: if x == nameFind: return True break return False def datetimeStr_datetime(StringDT): """ Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object >>> datetimeStr_datetime('1976-02-17_00:00:00') 1976-02-17 00:00:00 """ import datetime as dt fname = 'datetimeStr_datetime' dateD = np.zeros((3), dtype=int) timeT = np.zeros((3), dtype=int) # quit() dateD[0] = int(StringDT[0:4]) dateD[1] = int(StringDT[5:7]) dateD[2] = int(StringDT[8:10]) trefT = StringDT.find(':') if not trefT == -1: # print ' ' + fname + ': refdate with time!' timeT[0] = int(StringDT[11:13]) timeT[1] = int(StringDT[14:16]) timeT[2] = int(StringDT[17:19]) if int(dateD[0]) == 0: print warnmsg print ' ' + fname + ': 0 reference year!! changing to 1' dateD[0] = 1 newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2]) return newdatetime def dateStr_date(StringDate): """ Function to transform a string date ([YYYY]-[MM]-[DD] format) to a date object >>> dateStr_date('1976-02-17') 1976-02-17 """ import datetime as dt dateD = StringDate.split('-') if int(dateD[0]) == 0: print warnmsg print ' dateStr_date: 0 reference year!! changing to 1' dateD[0] = 1 newdate = dt.date(int(dateD[0]), int(dateD[1]), int(dateD[2])) return newdate def numVector_String(vec,char): """ Function to transform a vector of numbers to a single string [char] separated numVector_String(vec,char) vec= vector with the numerical values char= single character to split the values >>> print numVector_String(np.arange(10),' ') 0 1 2 3 4 5 6 7 8 9 """ fname = 'numVector_String' if vec == 'h': print fname + '_____________________________________________________________' print numVector_String.__doc__ quit() Nvals = len(vec) string='' for i in range(Nvals): if i == 0: string = str(vec[i]) else: string = string + char + str(vec[i]) return string def timeref_datetime(refd, timeval, tu): """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to datetime object refd: time of reference (as datetime object) timeval: time value (as [tu] from [tref]) tu: time units >>> timeref = date(1949,12,1,0,0,0) >>> timeref_datetime(timeref, 229784.36, hours) 1976-02-17 08:21:36 """ import datetime as dt import numpy as np ## Not in timedelta # if tu == 'years': # realdate = refdate + dt.timedelta(years=float(timeval)) # elif tu == 'months': # realdate = refdate + dt.timedelta(months=float(timeval)) if tu == 'weeks': realdate = refd + dt.timedelta(weeks=float(timeval)) elif tu == 'days': realdate = refd + dt.timedelta(days=float(timeval)) elif tu == 'hours': realdate = refd + dt.timedelta(hours=float(timeval)) elif tu == 'minutes': realdate = refd + dt.timedelta(minutes=float(timeval)) elif tu == 'seconds': realdate = refd + dt.timedelta(seconds=float(timeval)) elif tu == 'milliseconds': realdate = refd + dt.timedelta(milliseconds=float(timeval)) else: print errormsg print ' timeref_datetime: time units "' + tu + '" not ready!!!!' quit(-1) return realdate def slice_variable(varobj, dimslice): """ Function to return a slice of a given variable according to values to its dimensions slice_variable(varobj, dims) varobj= object wit the variable dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension [value]: * [integer]: which value of the dimension * -1: all along the dimension * [beg]:[end] slice from [beg] to [end] * -9: last value of the dimension """ fname = 'slice_variable' if varobj == 'h': print fname + '_____________________________________________________________' print slice_variable.__doc__ quit() vardims = varobj.dimensions Ndimvar = len(vardims) Ndimcut = len(dimslice.split('|')) if Ndimcut == 0: Ndimcut = 1 dimcut = list(dimslice) dimsl = dimslice.split('|') varvalsdim = [] dimnslice = [] monodim = [] for idd in range(Ndimvar): found = False for idc in range(Ndimcut): dimcutn = dimsl[idc].split(':')[0] dimcutv = dimsl[idc].split(':')[1] if vardims[idd] == dimcutn: posfrac = dimcutv.find('@') if posfrac != -1: inifrac = int(dimcutv.split('@')[0]) endfrac = int(dimcutv.split('@')[1]) varvalsdim.append(slice(inifrac,endfrac)) dimnslice.append(vardims[idd]) else: if int(dimcutv) == -1: varvalsdim.append(slice(0,varobj.shape[idd])) dimnslice.append(vardims[idd]) elif int(dimcutv) == -9: varvalsdim.append(varobj.shape[idd]-1) monodim.append(vardims[idd]) else: if int(dimcutv) > varobj.shape[idd]-1: print errormsg print ' ' + fname + ': given value:', int(dimcutv), \ "is larger than length of dimension '" + dimcutn + \ "':", varobj.shape[idd],'(C-like!) !!' quit(-1) varvalsdim.append(int(dimcutv)) monodim.append(vardims[idd]) found = True break if not found and not searchInlist(dimnslice,vardims[idd]) and \ not searchInlist(monodim,vardims[idd]): varvalsdim.append(slice(0,varobj.shape[idd])) dimnslice.append(vardims[idd]) varvalues = varobj[tuple(varvalsdim)] varvalues = np.squeeze(varobj[tuple(varvalsdim)]) return varvalues, dimnslice def get_str_nc(varo, Lchar): """ Function to get string values in a netCDF variable as a chains of 1char values varo= netCDF variable object Lchar = length of the string in the netCDF file """ fname = 'get_str_nc' Nvalues = varo.shape[0] strings = [] for ival in range(Nvalues): stringv = varo[ival,:] string = str('') for ich in range(Lchar): charval = str(stringv[ich:ich+1]).replace('[','').replace(']','').replace("'","") if charval == '--' or len(charval) == 0: break else: if ich == 0: string = charval else: string = string + charval strings.append(string.strip()) return strings def interpolate_locs(locs,coords,kinterp): """ Function to provide interpolate locations on a given axis interpolate_locs(locs,axis,kinterp) locs= locations to interpolate coords= axis values with the reference of coordinates kinterp: kind of interpolation 'lin': linear >>> coordinates = np.arange((10), dtype=np.float) >>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0]) >>> interpolate_locs(values,coordinates,'lin') [ -1.2 2.4 5.6 7.8 13. ] >>> coordinates[0] = 0.5 >>> coordinates[2] = 2.5 >>> interpolate_locs(values,coordinates,'lin') [ -3.4 1.93333333 5.6 7.8 13. ] """ fname = 'interpolate_locs' if locs == 'h': print fname + '_____________________________________________________________' print interpolate_locs.__doc__ quit() Nlocs = locs.shape[0] Ncoords = coords.shape[0] dcoords = coords[Ncoords-1] - coords[0] intlocs = np.zeros((Nlocs), dtype=np.float) minc = np.min(coords) maxc = np.max(coords) for iloc in range(Nlocs): for icor in range(Ncoords-1): if locs[iloc] < minc and dcoords > 0.: a = 0. b = 1. / (coords[1] - coords[0]) c = coords[0] elif locs[iloc] > maxc and dcoords > 0.: a = (Ncoords-1)*1. b = 1. / (coords[Ncoords-1] - coords[Ncoords-2]) c = coords[Ncoords-2] elif locs[iloc] < minc and dcoords < 0.: a = (Ncoords-1)*1. b = 1. / (coords[Ncoords-1] - coords[Ncoords-2]) c = coords[Ncoords-2] elif locs[iloc] > maxc and dcoords < 0.: a = 0. b = 1. / (coords[1] - coords[0]) c = coords[0] elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.: a = icor*1. b = 1. / (coords[icor+1] - coords[icor]) c = coords[icor] print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.: a = icor*1. b = 1. / (coords[icor+1] - coords[icor]) c = coords[icor] if kinterp == 'lin': intlocs[iloc] = a + (locs[iloc] - c)*b else: print errormsg print ' ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!" quit(-1) return intlocs def datetimeStr_conversion(StringDT,typeSi,typeSo): """ Function to transform a string date to an another date object StringDT= string with the date and time typeSi= type of datetime string input typeSo= type of datetime string output [typeSi/o] 'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate] 'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]] 'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format 'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format 'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format 'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format 'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H], [H], ':', [M], [M], ':', [S], [S] >>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS') [1976 2 17 8 32 5] >>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S') 1980/03/05 18-00-00 """ import datetime as dt fname = 'datetimeStr_conversion' if StringDT[0:1] == 'h': print fname + '_____________________________________________________________' print datetimeStr_conversion.__doc__ quit() if typeSi == 'cfTime': timeval = np.float(StringDT.split(',')[0]) tunits = StringDT.split(',')[1].split(' ')[0] Srefdate = StringDT.split(',')[1].split(' ')[2] # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] ## yrref=Srefdate[0:4] monref=Srefdate[5:7] dayref=Srefdate[8:10] trefT = Srefdate.find(':') if not trefT == -1: # print ' ' + fname + ': refdate with time!' horref=Srefdate[11:13] minref=Srefdate[14:16] secref=Srefdate[17:19] refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \ '_' + horref + ':' + minref + ':' + secref) else: refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \ '_00:00:00') if tunits == 'weeks': newdate = refdate + dt.timedelta(weeks=float(timeval)) elif tunits == 'days': newdate = refdate + dt.timedelta(days=float(timeval)) elif tunits == 'hours': newdate = refdate + dt.timedelta(hours=float(timeval)) elif tunits == 'minutes': newdate = refdate + dt.timedelta(minutes=float(timeval)) elif tunits == 'seconds': newdate = refdate + dt.timedelta(seconds=float(timeval)) elif tunits == 'milliseconds': newdate = refdate + dt.timedelta(milliseconds=float(timeval)) else: print errormsg print ' timeref_datetime: time units "' + tunits + '" not ready!!!!' quit(-1) yr = newdate.year mo = newdate.month da = newdate.day ho = newdate.hour mi = newdate.minute se = newdate.second elif typeSi == 'matYmdHMS': yr = StringDT[0] mo = StringDT[1] da = StringDT[2] ho = StringDT[3] mi = StringDT[4] se = StringDT[5] elif typeSi == 'YmdHMS': yr = int(StringDT[0:4]) mo = int(StringDT[4:6]) da = int(StringDT[6:8]) ho = int(StringDT[8:10]) mi = int(StringDT[10:12]) se = int(StringDT[12:14]) elif typeSi == 'Y-m-d_H:M:S': dateDT = StringDT.split('_') dateD = dateDT[0].split('-') timeT = dateDT[1].split(':') yr = int(dateD[0]) mo = int(dateD[1]) da = int(dateD[2]) ho = int(timeT[0]) mi = int(timeT[1]) se = int(timeT[2]) elif typeSi == 'Y-m-d H:M:S': dateDT = StringDT.split(' ') dateD = dateDT[0].split('-') timeT = dateDT[1].split(':') yr = int(dateD[0]) mo = int(dateD[1]) da = int(dateD[2]) ho = int(timeT[0]) mi = int(timeT[1]) se = int(timeT[2]) elif typeSi == 'Y/m/d H-M-S': dateDT = StringDT.split(' ') dateD = dateDT[0].split('/') timeT = dateDT[1].split('-') yr = int(dateD[0]) mo = int(dateD[1]) da = int(dateD[2]) ho = int(timeT[0]) mi = int(timeT[1]) se = int(timeT[2]) elif typeSi == 'WRFdatetime': yr = int(StringDT[0])*1000 + int(StringDT[1])*100 + int(StringDT[2])*10 + \ int(StringDT[3]) mo = int(StringDT[5])*10 + int(StringDT[6]) da = int(StringDT[8])*10 + int(StringDT[9]) ho = int(StringDT[11])*10 + int(StringDT[12]) mi = int(StringDT[14])*10 + int(StringDT[15]) se = int(StringDT[17])*10 + int(StringDT[18]) else: print errormsg print ' ' + fname + ': type of String input date "' + typeSi + \ '" not ready !!!!' quit(-1) if typeSo == 'matYmdHMS': dateYmdHMS = np.zeros((6), dtype=int) dateYmdHMS[0] = yr dateYmdHMS[1] = mo dateYmdHMS[2] = da dateYmdHMS[3] = ho dateYmdHMS[4] = mi dateYmdHMS[5] = se elif typeSo == 'YmdHMS': dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) + \ str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2) elif typeSo == 'Y-m-d_H:M:S': dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \ str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \ str(se).zfill(2) elif typeSo == 'Y-m-d H:M:S': dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \ str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \ str(se).zfill(2) elif typeSo == 'Y/m/d H-M-S': dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' + \ str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \ str(se).zfill(2) elif typeSo == 'WRFdatetime': dateYmdHMS = [] yM = yr/1000 yC = (yr-yM*1000)/100 yD = (yr-yM*1000-yC*100)/10 yU = yr-yM*1000-yC*100-yD*10 mD = mo/10 mU = mo-mD*10 dD = da/10 dU = da-dD*10 hD = ho/10 hU = ho-hD*10 miD = mi/10 miU = mi-miD*10 sD = se/10 sU = se-sD*10 dateYmdHMS.append(str(yM)) dateYmdHMS.append(str(yC)) dateYmdHMS.append(str(yD)) dateYmdHMS.append(str(yU)) dateYmdHMS.append('-') dateYmdHMS.append(str(mD)) dateYmdHMS.append(str(mU)) dateYmdHMS.append('-') dateYmdHMS.append(str(dD)) dateYmdHMS.append(str(dU)) dateYmdHMS.append('_') dateYmdHMS.append(str(hD)) dateYmdHMS.append(str(hU)) dateYmdHMS.append(':') dateYmdHMS.append(str(miD)) dateYmdHMS.append(str(miU)) dateYmdHMS.append(':') dateYmdHMS.append(str(sD)) dateYmdHMS.append(str(sU)) else: print errormsg print ' ' + fname + ': type of output date "' + typeSo + '" not ready !!!!' quit(-1) return dateYmdHMS def percendone(nvals,tot,percen,msg): """ Function to provide the percentage of an action across the matrix nvals=number of values tot=total number of values percen=percentage frequency for which the message is wanted msg= message """ from sys import stdout num = int(tot * percen/100) if (nvals%num == 0): print '\r ' + msg + '{0:8.3g}'.format(nvals*100./tot) + ' %', stdout.flush() return '' def netCDFdatetime_realdatetime(units, tcalendar, times): """ Function to transfrom from netCDF CF-compilant times to real time """ import datetime as dt txtunits = units.split(' ') tunits = txtunits[0] Srefdate = txtunits[len(txtunits) - 1] # Calendar type ## is360 = False if tcalendar is not None: print ' netCDFdatetime_realdatetime: There is a calendar attribute' if tcalendar == '365_day' or tcalendar == 'noleap': print ' netCDFdatetime_realdatetime: No leap years!' isleapcal = False elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian': isleapcal = True elif tcalendar == '360_day': is360 = True isleapcal = False else: print errormsg print ' netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!' quit(-1) # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] ## timeval = Srefdate.find(':') if not timeval == -1: print ' netCDFdatetime_realdatetime: refdate with time!' refdate = datetimeStr_datetime(Srefdate) else: refdate = dateStr_date(Srefdate + '_00:00:00') dimt = len(times) # datetype = type(dt.datetime(1972,02,01)) # realdates = np.array(dimt, datetype) # print realdates ## Not in timedelta # if tunits == 'years': # for it in range(dimt): # realdate = refdate + dt.timedelta(years=float(times[it])) # realdates[it] = int(realdate.year) # elif tunits == 'months': # for it in range(dimt): # realdate = refdate + dt.timedelta(months=float(times[it])) # realdates[it] = int(realdate.year) # realdates = [] realdates = np.zeros((dimt, 6), dtype=int) if tunits == 'weeks': for it in range(dimt): realdate = refdate + dt.timedelta(weeks=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'days': for it in range(dimt): realdate = refdate + dt.timedelta(days=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'hours': for it in range(dimt): realdate = refdate + dt.timedelta(hours=float(times[it])) # if not isleapcal: # Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year)) # realdate = realdate - dt.timedelta(days=Nleapdays) # if is360: # Nyears360 = int(realdate.year) - int(refdate.year) + 1 # realdate = realdate -dt.timedelta(days=Nyears360*5) # realdates[it] = realdate # realdates = refdate + dt.timedelta(hours=float(times)) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'minutes': for it in range(dimt): realdate = refdate + dt.timedelta(minutes=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'seconds': for it in range(dimt): realdate = refdate + dt.timedelta(seconds=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'milliseconds': for it in range(dimt): realdate = refdate + dt.timedelta(milliseconds=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] elif tunits == 'microseconds': for it in range(dimt): realdate = refdate + dt.timedelta(microseconds=float(times[it])) realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second] else: print errormsg print ' netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!' quit(-1) return realdates def file_nlines(filen): """ Function to provide the number of lines of a file filen= name of the file >>> file_nlines('trajectory.dat') 49 """ fname = 'file_nlines' if not os.path.isfile(filen): print errormsg print ' ' + fname + ' file: "' + filen + '" does not exist !!' quit(-1) fo = open(filen,'r') nlines=0 for line in fo: nlines = nlines + 1 fo.close() return nlines def realdatetime1_CFcompilant(time, Srefdate, tunits): """ Function to transform a matrix with a real time value ([year, month, day, hour, minute, second]) to a netCDF one time= matrix with time Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format) tunits= units of time respect to Srefdate >>> realdatetime1_CFcompilant([1976, 2, 17, 8, 20, 0], '19491201000000', 'hours') 229784.33333333 """ import datetime as dt yrref=int(Srefdate[0:4]) monref=int(Srefdate[4:6]) dayref=int(Srefdate[6:8]) horref=int(Srefdate[8:10]) minref=int(Srefdate[10:12]) secref=int(Srefdate[12:14]) refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref) if tunits == 'weeks': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5])-refdate cfdates = (cfdate.days + cfdate.seconds/(3600.*24.))/7. elif tunits == 'days': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate cfdates = cfdate.days + cfdate.seconds/(3600.*24.) elif tunits == 'hours': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate cfdates = cfdate.days*24. + cfdate.seconds/3600. elif tunits == 'minutes': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate cfdates = cfdate.days*24.*60. + cfdate.seconds/60. elif tunits == 'seconds': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate cfdates = cfdate.days*24.*3600. + cfdate.seconds elif tunits == 'milliseconds': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate cfdates = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000. elif tunits == 'microseconds': cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],times[5]) - refdate cfdates = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000. else: print errormsg print ' ' + fname + ': time units "' + tunits + '" is not ready!!!' quit(-1) return cfdates def basicvardef(varobj, vstname, vlname, vunits): """ Function to give the basic attributes to a variable varobj= netCDF variable object vstname= standard name of the variable vlname= long name of the variable vunits= units of the variable """ attr = varobj.setncattr('standard_name', vstname) attr = varobj.setncattr('long_name', vlname) attr = varobj.setncattr('units', vunits) return def variables_values(varName): """ Function to provide values to plot the different variables values from ASCII file 'variables_values.dat' variables_values(varName) [varName]= name of the variable return: [var name], [std name], [minimum], [maximum], [long name]('|' for spaces), [units], [color palette] (following: http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html) [varn]: original name of the variable NOTE: It might be better doing it with an external ASII file. But then we got an extra dependency... >>> variables_values('WRFght') ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow'] """ import subprocess as sub fname='variables_values' if varName == 'h': print fname + '_____________________________________________________________' print variables_values.__doc__ quit() # This does not work.... # folderins = sub.Popen(["pwd"], stdout=sub.PIPE) # folder = list(folderins.communicate())[0].replace('\n','') # From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python folder = os.path.dirname(os.path.realpath(__file__)) infile = folder + '/variables_values.dat' if not os.path.isfile(infile): print errormsg print ' ' + fname + ": File '" + infile + "' does not exist !!" quit(-1) # Variable name might come with a statistical surname... stats=['min','max','mean','stdv', 'sum'] # Variables with a statistical section on their name... NOstatsvars = ['zmaxth', 'zmax_th', 'lmax_th', 'lmaxth'] ifst = False if not searchInlist(NOstatsvars, varName.lower()): for st in stats: if varName.find(st) > -1: print ' '+ fname + ": varibale '" + varName + "' with a " + \ "statistical surname: '",st,"' !!" Lst = len(st) LvarName = len(varName) varn = varName[0:LvarName - Lst] ifst = True break if not ifst: varn = varName ncf = open(infile, 'r') for line in ncf: if line[0:1] != '#': values = line.replace('\n','').split(',') if len(values) != 8: print errormsg print "problem in varibale: '", values[0], \ "' it should have 8 values and it has:", len(values) quit(-1) if varn[0:6] == 'varDIM': # Variable from a dimension (all with 'varDIM' prefix) Lvarn = len(varn) varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1., \ "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1', \ 'rainbow'] else: varvals = [values[1].replace(' ',''), values[2].replace(' ',''), \ np.float(values[3]), np.float(values[4]),values[5].replace(' ',''),\ values[6].replace(' ',''), values[7].replace(' ','')] if values[0] == varn: ncf.close() return varvals break print errormsg print ' ' + fname + ": variable '" + varn + "' not defined !!!" ncf.close() quit(-1) return def variables_values_old(varName): """ Function to provide values to plot the different variables variables_values(varName) [varName]= name of the variable return: [var name], [std name], [minimum], [maximum], [long name]('|' for spaces), [units], [color palette] (following: http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html) [varn]: original name of the variable NOTE: It might be better doing it with an external ASII file. But then we got an extra dependency... >>> variables_values('WRFght') ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow'] """ fname='variables_values' if varName == 'h': print fname + '_____________________________________________________________' print variables_values.__doc__ quit() # Variable name might come with a statistical surname... stats=['min','max','mean','stdv', 'sum'] ifst = False for st in stats: if varName.find(st) > -1: print ' '+ fname + ": varibale '" + varName + "' with a statistical "+\ " surname: '",st,"' !!" Lst = len(st) LvarName = len(varName) varn = varName[0:LvarName - Lst] ifst = True break if not ifst: varn = varName if varn[0:6] == 'varDIM': # Variable from a dimension (all with 'varDIM' prefix) Lvarn = len(varn) varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1., \ "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1', 'rainbox'] elif varn == 'a_tht' or varn == 'LA_THT': varvals = ['ath', 'total_thermal_plume_cover', 0., 1., \ 'total|column|thermal|plume|cover', '1', 'YlGnBu'] elif varn == 'acprc' or varn == 'RAINC': varvals = ['acprc', 'accumulated_cmulus_precipitation', 0., 3.e4, \ 'accumulated|cmulus|precipitation', 'mm', 'Blues'] elif varn == 'acprnc' or varn == 'RAINNC': varvals = ['acprnc', 'accumulated_non-cmulus_precipitation', 0., 3.e4, \ 'accumulated|non-cmulus|precipitation', 'mm', 'Blues'] elif varn == 'bils' or varn == 'LBILS': varvals = ['bils', 'surface_total_heat_flux', -100., 100., \ 'surface|total|heat|flux', 'Wm-2', 'seismic'] elif varn == 'landcat' or varn == 'category': varvals = ['landcat', 'land_categories', 0., 22., 'land|categories', '1', \ 'rainbow'] elif varn == 'c' or varn == 'QCLOUD' or varn == 'oliq' or varn == 'OLIQ': varvals = ['c', 'condensed_water_mixing_ratio', 0., 3.e-4, \ 'condensed|water|mixing|ratio', 'kgkg-1', 'BuPu'] elif varn == 'ci' or varn == 'iwcon' or varn == 'LIWCON': varvals = ['ci', 'cloud_iced_water_mixing_ratio', 0., 0.0003, \ 'cloud|iced|water|mixing|ratio', 'kgkg-1', 'Purples'] elif varn == 'cl' or varn == 'lwcon' or varn == 'LLWCON': varvals = ['cl', 'cloud_liquidwater_mixing_ratio', 0., 0.0003, \ 'cloud|liquid|water|mixing|ratio', 'kgkg-1', 'Blues'] elif varn == 'cld' or varn == 'CLDFRA' or varn == 'rneb' or varn == 'lrneb' or \ varn == 'LRNEB': varvals = ['cld', 'cloud_area_fraction', 0., 1., 'cloud|fraction', '1', \ 'gist_gray'] elif varn == 'cldc' or varn == 'rnebcon' or varn == 'lrnebcon' or \ varn == 'LRNEBCON': varvals = ['cldc', 'convective_cloud_area_fraction', 0., 1., \ 'convective|cloud|fraction', '1', 'gist_gray'] elif varn == 'cldl' or varn == 'rnebls' or varn == 'lrnebls' or varn == 'LRNEBLS': varvals = ['cldl', 'large_scale_cloud_area_fraction', 0., 1., \ 'large|scale|cloud|fraction', '1', 'gist_gray'] elif varn == 'clt' or varn == 'CLT' or varn == 'cldt' or \ varn == 'Total cloudiness': varvals = ['clt', 'cloud_area_fraction', 0., 1., 'total|cloud|cover', '1', \ 'gist_gray'] elif varn == 'cll' or varn == 'cldl' or varn == 'LCLDL' or \ varn == 'Low-level cloudiness': varvals = ['cll', 'low_level_cloud_area_fraction', 0., 1., \ 'low|level|(p|>|680|hPa)|cloud|fraction', '1', 'gist_gray'] elif varn == 'clm' or varn == 'cldm' or varn == 'LCLDM' or \ varn == 'Mid-level cloudiness': varvals = ['clm', 'mid_level_cloud_area_fraction', 0., 1., \ 'medium|level|(440|<|p|<|680|hPa)|cloud|fraction', '1', 'gist_gray'] elif varn == 'clh' or varn == 'cldh' or varn == 'LCLDH' or \ varn == 'High-level cloudiness': varvals = ['clh', 'high_level_cloud_area_fraction', 0., 1., \ 'high|level|(p|<|440|hPa)|cloud|fraction', '1', 'gist_gray'] elif varn == 'clmf' or varn == 'fbase' or varn == 'LFBASE': varvals = ['clmf', 'cloud_base_max_flux', -0.3, 0.3, 'cloud|base|max|flux', \ 'kgm-2s-1', 'seismic'] elif varn == 'clp' or varn == 'pbase' or varn == 'LPBASE': varvals = ['clp', 'cloud_base_pressure', -0.3, 0.3, 'cloud|base|pressure', \ 'Pa', 'Reds'] elif varn == 'cpt' or varn == 'ptconv' or varn == 'LPTCONV': varvals = ['cpt', 'convective_point', 0., 1., 'convective|point', '1', \ 'seismic'] elif varn == 'dqajs' or varn == 'LDQAJS': varvals = ['dqajs', 'dry_adjustment_water_vapor_tendency', -0.0003, 0.0003, \ 'dry|adjustment|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqcon' or varn == 'LDQCON': varvals = ['dqcon', 'convective_water_vapor_tendency', -3e-8, 3.e-8, \ 'convective|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqdyn' or varn == 'LDQDYN': varvals = ['dqdyn', 'dynamics_water_vapor_tendency', -3.e-7, 3.e-7, \ 'dynamics|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqeva' or varn == 'LDQEVA': varvals = ['dqeva', 'evaporation_water_vapor_tendency', -3.e-6, 3.e-6, \ 'evaporation|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqlscst' or varn == 'LDQLSCST': varvals = ['dqlscst', 'stratocumulus_water_vapor_tendency', -3.e-7, 3.e-7, \ 'stratocumulus|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqlscth' or varn == 'LDQLSCTH': varvals = ['dqlscth', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7, \ 'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqlsc' or varn == 'LDQLSC': varvals = ['dqlsc', 'condensation_water_vapor_tendency', -3.e-6, 3.e-6, \ 'condensation|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqphy' or varn == 'LDQPHY': varvals = ['dqphy', 'physics_water_vapor_tendency', -3.e-7, 3.e-7, \ 'physics|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqthe' or varn == 'LDQTHE': varvals = ['dqthe', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7, \ 'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqvdf' or varn == 'LDQVDF': varvals = ['dqvdf', 'vertical_difussion_water_vapor_tendency', -3.e-8, 3.e-8,\ 'vertical|difussion|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dqwak' or varn == 'LDQWAK': varvals = ['dqwak', 'wake_water_vapor_tendency', -3.e-7, 3.e-7, \ 'wake|water|vapor|tendency', 'kg/kg/s', 'seismic'] elif varn == 'dta' or varn == 'tnt' or varn == 'LTNT': varvals = ['dta', 'tendency_air_temperature', -3.e-3, 3.e-3, \ 'tendency|of|air|temperature', 'K/s', 'seismic'] elif varn == 'dtac' or varn == 'tntc' or varn == 'LTNTC': varvals = ['dtac', 'moist_convection_tendency_air_temperature', -3.e-3, \ 3.e-3, 'moist|convection|tendency|of|air|temperature', 'K/s', 'seismic'] elif varn == 'dtar' or varn == 'tntr' or varn == 'LTNTR': varvals = ['dtar', 'radiative_heating_tendency_air_temperature', -3.e-3, \ 3.e-3, 'radiative|heating|tendency|of|air|temperature', 'K/s', 'seismic'] elif varn == 'dtascpbl' or varn == 'tntscpbl' or varn == 'LTNTSCPBL': varvals = ['dtascpbl', \ 'stratiform_cloud_precipitation_BL_mixing_tendency_air_temperature', \ -3.e-6, 3.e-6, \ 'stratiform|cloud|precipitation|Boundary|Layer|mixing|tendency|air|' + 'temperature', 'K/s', 'seismic'] elif varn == 'dtajs' or varn == 'LDTAJS': varvals = ['dtajs', 'dry_adjustment_thermal_tendency', -3.e-5, 3.e-5, \ 'dry|adjustment|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtcon' or varn == 'LDTCON': varvals = ['dtcon', 'convective_thermal_tendency', -3.e-5, 3.e-5, \ 'convective|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtdyn' or varn == 'LDTDYN': varvals = ['dtdyn', 'dynamics_thermal_tendency', -3.e-4, 3.e-4, \ 'dynamics|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dteva' or varn == 'LDTEVA': varvals = ['dteva', 'evaporation_thermal_tendency', -3.e-3, 3.e-3, \ 'evaporation|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtlscst' or varn == 'LDTLSCST': varvals = ['dtlscst', 'stratocumulus_thermal_tendency', -3.e-4, 3.e-4, \ 'stratocumulus|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtlscth' or varn == 'LDTLSCTH': varvals = ['dtlscth', 'thermals_thermal_tendency', -3.e-4, 3.e-4, \ 'thermal|plumes|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtlsc' or varn == 'LDTLSC': varvals = ['dtlsc', 'condensation_thermal_tendency', -3.e-3, 3.e-3, \ 'condensation|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtlwr' or varn == 'LDTLWR': varvals = ['dtlwr', 'long_wave_thermal_tendency', -3.e-3, 3.e-3, \ 'long|wave|radiation|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtphy' or varn == 'LDTPHY': varvals = ['dtphy', 'physics_thermal_tendency', -3.e-4, 3.e-4, \ 'physics|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtsw0' or varn == 'LDTSW0': varvals = ['dtsw0', 'cloudy_sky_short_wave_thermal_tendency', -3.e-4, 3.e-4, \ 'cloudy|sky|short|wave|radiation|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtthe' or varn == 'LDTTHE': varvals = ['dtthe', 'thermals_thermal_tendency', -3.e-4, 3.e-4, \ 'thermal|plumes|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtvdf' or varn == 'LDTVDF': varvals = ['dtvdf', 'vertical_difussion_thermal_tendency', -3.e-5, 3.e-5, \ 'vertical|difussion|thermal|tendency', 'K/s', 'seismic'] elif varn == 'dtwak' or varn == 'LDTWAK': varvals = ['dtwak', 'wake_thermal_tendency', -3.e-4, 3.e-4, \ 'wake|thermal|tendency', 'K/s', 'seismic'] elif varn == 'ducon' or varn == 'LDUCON': varvals = ['ducon', 'convective_eastward_wind_tendency', -3.e-3, 3.e-3, \ 'convective|eastward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'dudyn' or varn == 'LDUDYN': varvals = ['dudyn', 'dynamics_eastward_wind_tendency', -3.e-3, 3.e-3, \ 'dynamics|eastward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'duvdf' or varn == 'LDUVDF': varvals = ['duvdf', 'vertical_difussion_eastward_wind_tendency', -3.e-3, \ 3.e-3, 'vertical|difussion|eastward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'dvcon' or varn == 'LDVCON': varvals = ['dvcon', 'convective_difussion_northward_wind_tendency', -3.e-3, \ 3.e-3, 'convective|northward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'dvdyn' or varn == 'LDVDYN': varvals = ['dvdyn', 'dynamics_northward_wind_tendency', -3.e-3, \ 3.e-3, 'dynamics|difussion|northward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'dvvdf' or varn == 'LDVVDF': varvals = ['dvvdf', 'vertical_difussion_northward_wind_tendency', -3.e-3, \ 3.e-3, 'vertical|difussion|northward|wind|tendency', 'ms-2', 'seismic'] elif varn == 'etau' or varn == 'ZNU': varvals = ['etau', 'etau', 0., 1, 'eta values on half (mass) levels', '-', \ 'reds'] elif varn == 'evspsbl' or varn == 'LEVAP' or varn == 'evap' or varn == 'SFCEVPde': varvals = ['evspsbl', 'water_evaporation_flux', 0., 1.5e-4, \ 'water|evaporation|flux', 'kgm-2s-1', 'Blues'] elif varn == 'evspsbl' or varn == 'SFCEVPde': varvals = ['evspsblac', 'water_evaporation_flux_ac', 0., 1.5e-4, \ 'accumulated|water|evaporation|flux', 'kgm-2', 'Blues'] elif varn == 'g' or varn == 'QGRAUPEL': varvals = ['g', 'grauepl_mixing_ratio', 0., 0.0003, 'graupel|mixing|ratio', \ 'kgkg-1', 'Purples'] elif varn == 'h2o' or varn == 'LH2O': varvals = ['h2o', 'water_mass_fraction', 0., 3.e-2, \ 'mass|fraction|of|water', '1', 'Blues'] elif varn == 'h' or varn == 'QHAIL': varvals = ['h', 'hail_mixing_ratio', 0., 0.0003, 'hail|mixing|ratio', \ 'kgkg-1', 'Purples'] elif varn == 'hfls' or varn == 'LH' or varn == 'LFLAT' or varn == 'flat': varvals = ['hfls', 'surface_upward_latent_heat_flux', -400., 400., \ 'upward|latnt|heat|flux|at|the|surface', 'Wm-2', 'seismic'] elif varn == 'hfss' or varn == 'LSENS' or varn == 'sens' or varn == 'HFX': varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150., \ 'upward|sensible|heat|flux|at|the|surface', 'Wm-2', 'seismic'] elif varn == 'hfso' or varn == 'GRDFLX': varvals = ['hfso', 'downward_heat_flux_in_soil', -150., 150., \ 'Downward|soil|heat|flux', 'Wm-2', 'seismic'] elif varn == 'hus' or varn == 'WRFrh' or varn == 'LMDZrh' or varn == 'rhum' or \ varn == 'LRHUM': varvals = ['hus', 'specific_humidity', 0., 1., 'specific|humidty', '1', \ 'BuPu'] elif varn == 'huss' or varn == 'WRFrhs' or varn == 'LMDZrhs' or varn == 'rh2m' or\ varn == 'LRH2M': varvals = ['huss', 'specific_humidity', 0., 1., 'specific|humidty|at|2m', \ '1', 'BuPu'] elif varn == 'i' or varn == 'QICE': varvals = ['i', 'iced_water_mixing_ratio', 0., 0.0003, \ 'iced|water|mixing|ratio', 'kgkg-1', 'Purples'] elif varn == 'lat' or varn == 'XLAT' or varn == 'XLAT_M' or varn == 'latitude': varvals = ['lat', 'latitude', -90., 90., 'latitude', 'degrees North', \ 'seismic'] elif varn == 'lcl' or varn == 's_lcl' or varn == 'ls_lcl' or varn == 'LS_LCL': varvals = ['lcl', 'condensation_level', 0., 2500., 'level|of|condensation', \ 'm', 'Greens'] elif varn == 'lambdath' or varn == 'lambda_th' or varn == 'LLAMBDA_TH': varvals = ['lambdath', 'thermal_plume_vertical_velocity', -30., 30., \ 'thermal|plume|vertical|velocity', 'm/s', 'seismic'] elif varn == 'lmaxth' or varn == 'LLMAXTH': varvals = ['lmaxth', 'upper_level_thermals', 0., 100., 'upper|level|thermals'\ , '1', 'Greens'] elif varn == 'lon' or varn == 'XLONG' or varn == 'XLONG_M': varvals = ['lon', 'longitude', -180., 180., 'longitude', 'degrees East', \ 'seismic'] elif varn == 'longitude': varvals = ['lon', 'longitude', 0., 360., 'longitude', 'degrees East', \ 'seismic'] elif varn == 'orog' or varn == 'HGT' or varn == 'HGT_M': varvals = ['orog', 'orography', 0., 3000., 'surface|altitude', 'm','terrain'] elif varn == 'pfc' or varn == 'plfc' or varn == 'LPLFC': varvals = ['pfc', 'pressure_free_convection', 100., 1100., \ 'pressure|free|convection', 'hPa', 'BuPu'] elif varn == 'plcl' or varn == 'LPLCL': varvals = ['plcl', 'pressure_lifting_condensation_level', 700., 1100., \ 'pressure|lifting|condensation|level', 'hPa', 'BuPu'] elif varn == 'pr' or varn == 'RAINTOT' or varn == 'precip' or \ varn == 'LPRECIP' or varn == 'Precip Totale liq+sol': varvals = ['pr', 'precipitation_flux', 0., 1.e-4, 'precipitation|flux', \ 'kgm-2s-1', 'BuPu'] elif varn == 'prprof' or varn == 'vprecip' or varn == 'LVPRECIP': varvals = ['prprof', 'precipitation_profile', 0., 1.e-3, \ 'precipitation|profile', 'kg/m2/s', 'BuPu'] elif varn == 'prprofci' or varn == 'pr_con_i' or varn == 'LPR_CON_I': varvals = ['prprofci', 'precipitation_profile_convective_i', 0., 1.e-3, \ 'precipitation|profile|convective|i', 'kg/m2/s', 'BuPu'] elif varn == 'prprofcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L': varvals = ['prprofcl', 'precipitation_profile_convective_l', 0., 1.e-3, \ 'precipitation|profile|convective|l', 'kg/m2/s', 'BuPu'] elif varn == 'prprofli' or varn == 'pr_lsc_i' or varn == 'LPR_LSC_I': varvals = ['prprofli', 'precipitation_profile_large_scale_i', 0., 1.e-3, \ 'precipitation|profile|large|scale|i', 'kg/m2/s', 'BuPu'] elif varn == 'prprofll' or varn == 'pr_lsc_l' or varn == 'LPR_LSC_L': varvals = ['prprofll', 'precipitation_profile_large_scale_l', 0., 1.e-3, \ 'precipitation|profile|large|scale|l', 'kg/m2/s', 'BuPu'] elif varn == 'pracc' or varn == 'ACRAINTOT': varvals = ['pracc', 'precipitation_amount', 0., 100., \ 'accumulated|precipitation', 'kgm-2', 'BuPu'] elif varn == 'prc' or varn == 'LPLUC' or varn == 'pluc' or varn == 'WRFprc' or \ varn == 'RAINCde': varvals = ['prc', 'convective_precipitation_flux', 0., 2.e-4, \ 'convective|precipitation|flux', 'kgm-2s-1', 'Blues'] elif varn == 'prci' or varn == 'pr_con_i' or varn == 'LPR_CON_I': varvals = ['prci', 'convective_ice_precipitation_flux', 0., 0.003, \ 'convective|ice|precipitation|flux', 'kgm-2s-1', 'Purples'] elif varn == 'prcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L': varvals = ['prcl', 'convective_liquid_precipitation_flux', 0., 0.003, \ 'convective|liquid|precipitation|flux', 'kgm-2s-1', 'Blues'] elif varn == 'pres' or varn == 'presnivs' or varn == 'pressure' or \ varn == 'lpres' or varn == 'LPRES': varvals = ['pres', 'air_pressure', 0., 103000., 'air|pressure', 'Pa', \ 'Blues'] elif varn == 'prls' or varn == 'WRFprls' or varn == 'LPLUL' or varn == 'plul' or \ varn == 'RAINNCde': varvals = ['prls', 'large_scale_precipitation_flux', 0., 2.e-4, \ 'large|scale|precipitation|flux', 'kgm-2s-1', 'Blues'] elif varn == 'prsn' or varn == 'SNOW' or varn == 'snow' or varn == 'LSNOW': varvals = ['prsn', 'snowfall', 0., 1.e-4, 'snowfall|flux', 'kgm-2s-1', 'BuPu'] elif varn == 'prw' or varn == 'WRFprh': varvals = ['prw', 'atmosphere_water_vapor_content', 0., 10., \ 'water|vapor"path', 'kgm-2', 'Blues'] elif varn == 'ps' or varn == 'psfc' or varn =='PSFC' or varn == 'psol' or \ varn == 'Surface Pressure': varvals=['ps', 'surface_air_pressure', 85000., 105400., 'surface|pressure', \ 'hPa', 'cool'] elif varn == 'psl' or varn == 'mslp' or varn =='WRFmslp': varvals=['psl', 'air_pressure_at_sea_level', 85000., 104000., \ 'mean|sea|level|pressure', 'Pa', 'Greens'] elif varn == 'qth' or varn == 'q_th' or varn == 'LQ_TH': varvals = ['qth', 'thermal_plume_total_water_content', 0., 25., \ 'total|water|cotent|in|thermal|plume', 'mm', 'YlOrRd'] elif varn == 'r' or varn == 'QVAPOR' or varn == 'ovap' or varn == 'LOVAP': varvals = ['r', 'water_mixing_ratio', 0., 0.03, 'water|mixing|ratio', \ 'kgkg-1', 'BuPu'] elif varn == 'r2' or varn == 'Q2': varvals = ['r2', 'water_mixing_ratio_at_2m', 0., 0.03, 'water|mixing|' + \ 'ratio|at|2|m','kgkg-1', 'BuPu'] elif varn == 'rsds' or varn == 'SWdnSFC' or varn == 'SWdn at surface' or \ varn == 'SWDOWN': varvals=['rsds', 'surface_downwelling_shortwave_flux_in_air', 0., 1200., \ 'downward|SW|surface|radiation', 'Wm-2' ,'Reds'] elif varn == 'rsdsacc': varvals=['rsdsacc', 'accumulated_surface_downwelling_shortwave_flux_in_air', \ 0., 1200., 'accumulated|downward|SW|surface|radiation', 'Wm-2' ,'Reds'] elif varn == 'rvor' or varn == 'WRFrvor': varvals = ['rvor', 'air_relative_vorticity', -2.5E-3, 2.5E-3, \ 'air|relative|vorticity', 's-1', 'seismic'] elif varn == 'rvors' or varn == 'WRFrvors': varvals = ['rvors', 'surface_air_relative_vorticity', -2.5E-3, 2.5E-3, \ 'surface|air|relative|vorticity', 's-1', 'seismic'] elif varn == 's' or varn == 'QSNOW': varvals = ['s', 'snow_mixing_ratio', 0., 0.0003, 'snow|mixing|ratio', \ 'kgkg-1', 'Purples'] elif varn == 'stherm' or varn == 'LS_THERM': varvals = ['stherm', 'thermals_excess', 0., 0.8, 'thermals|excess', 'K', \ 'Reds'] elif varn == 'ta' or varn == 'WRFt' or varn == 'temp' or varn == 'LTEMP' or \ varn == 'Air temperature': varvals = ['ta', 'air_temperature', 195., 320., 'air|temperature', 'K', \ 'YlOrRd'] elif varn == 'tah' or varn == 'theta' or varn == 'LTHETA': varvals = ['tah', 'potential_air_temperature', 195., 320., \ 'potential|air|temperature', 'K', 'YlOrRd'] elif varn == 'tas' or varn == 'T2' or varn == 't2m' or varn == 'T2M' or \ varn == 'Temperature 2m': varvals = ['tas', 'air_temperature', 240., 310., 'air|temperature|at|2m', ' \ K', 'YlOrRd'] elif varn == 'tds' or varn == 'TH2': varvals = ['tds', 'air_dew_point_temperature', 240., 310., \ 'air|dew|point|temperature|at|2m', 'K', 'YlGnBu'] elif varn == 'tke' or varn == 'TKE' or varn == 'tke' or varn == 'LTKE': varvals = ['tke', 'turbulent_kinetic_energy', 0., 0.003, \ 'turbulent|kinetic|energy', 'm2/s2', 'Reds'] elif varn == 'time'or varn == 'time_counter': varvals = ['time', 'time', 0., 1000., 'time', \ 'hours|since|1949/12/01|00:00:00', 'Reds'] elif varn == 'tmla' or varn == 's_pblt' or varn == 'LS_PBLT': varvals = ['tmla', 'atmosphere_top_boundary_layer_temperature', 250., 330., \ 'atmosphere|top|boundary|layer|temperature', 'K', 'Reds'] elif varn == 'ua' or varn == 'vitu' or varn == 'U' or varn == 'Zonal wind' or \ varn == 'LVITU': varvals = ['ua', 'eastward_wind', -30., 30., 'eastward|wind', 'ms-1', \ 'seismic'] elif varn == 'uas' or varn == 'u10m' or varn == 'U10' or varn =='Vent zonal 10m': varvals = ['uas', 'eastward_wind', -30., 30., 'eastward|2m|wind', \ 'ms-1', 'seismic'] elif varn == 'va' or varn == 'vitv' or varn == 'V' or varn == 'Meridional wind' \ or varn == 'LVITV': varvals = ['va', 'northward_wind', -30., 30., 'northward|wind', 'ms-1', \ 'seismic'] elif varn == 'vas' or varn == 'v10m' or varn == 'V10' or \ varn =='Vent meridien 10m': varvals = ['vas', 'northward_wind', -30., 30., 'northward|2m|wind', 'ms-1', \ 'seismic'] elif varn == 'wakedeltaq' or varn == 'wake_deltaq' or varn == 'lwake_deltaq' or \ varn == 'LWAKE_DELTAQ': varvals = ['wakedeltaq', 'wake_delta_vapor', -0.003, 0.003, \ 'wake|delta|mixing|ratio', '-', 'seismic'] elif varn == 'wakedeltat' or varn == 'wake_deltat' or varn == 'lwake_deltat' or \ varn == 'LWAKE_DELTAT': varvals = ['wakedeltat', 'wake_delta_temp', -0.003, 0.003, \ 'wake|delta|temperature', '-', 'seismic'] elif varn == 'wakeh' or varn == 'wake_h' or varn == 'LWAKE_H': varvals = ['wakeh', 'wake_height', 0., 1000., 'height|of|the|wakes', 'm', \ 'YlOrRd'] elif varn == 'wakeomg' or varn == 'wake_omg' or varn == 'lwake_omg' or \ varn == 'LWAKE_OMG': varvals = ['wakeomg', 'wake_omega', 0., 3., 'wake|omega', \ '-', 'BuGn'] elif varn == 'wakes' or varn == 'wake_s' or varn == 'LWAKE_S': varvals = ['wakes', 'wake_area_fraction', 0., 0.5, 'wake|spatial|fraction', \ '1', 'BuGn'] elif varn == 'wa' or varn == 'W' or varn == 'Vertical wind': varvals = ['wa', 'upward_wind', -10., 10., 'upward|wind', 'ms-1', \ 'seismic'] elif varn == 'wap' or varn == 'vitw' or varn == 'LVITW': varvals = ['wap', 'upward_wind', -3.e-10, 3.e-10, 'upward|wind', 'mPa-1', \ 'seismic'] elif varn == 'wss' or varn == 'SPDUV': varvals = ['wss', 'air_velocity', 0., 30., 'surface|horizontal|wind|speed', \ 'ms-1', 'Reds'] # Water budget # Water budget de-accumulated elif varn == 'ccond' or varn == 'CCOND' or varn == 'ACCCONDde': varvals = ['ccond', 'cw_cond', 0., 30., \ 'cloud|water|condensation', 'mm', 'Reds'] elif varn == 'wbr' or varn == 'ACQVAPORde': varvals = ['wbr', 'wbr', 0., 30., 'Water|Budget|water|wapor', 'mm', 'Blues'] elif varn == 'diabh' or varn == 'DIABH' or varn == 'ACDIABHde': varvals = ['diabh', 'diabh', 0., 30., 'diabatic|heating', 'K', 'Reds'] elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde': varvals = ['wbpw', 'water_budget_pw', 0., 30., 'Water|Budget|water|content',\ 'mms-1', 'Reds'] elif varn == 'wbf' or varn == 'WBACF' or varn == 'WBACFde': varvals = ['wbf', 'water_budget_hfcqv', 0., 30., \ 'Water|Budget|horizontal|convergence|of|water|vapour|(+,|' + \ 'conv.;|-,|div.)', 'mms-1', 'Reds'] elif varn == 'wbfc' or varn == 'WBFC' or varn == 'WBACFCde': varvals = ['wbfc', 'water_budget_fc', 0., 30., \ 'Water|Budget|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\ 'div.)', 'mms-1', 'Reds'] elif varn == 'wbfp' or varn == 'WBFP' or varn == 'WBACFPde': varvals = ['wbfp', 'water_budget_cfp', 0., 30., \ 'Water|Budget|horizontal|convergence|of|precipitation|(+,|' + \ 'conv.;|-,|div.)', 'mms-1', 'Reds'] elif varn == 'wbz' or varn == 'WBZ' or varn == 'WBACZde': varvals = ['wbz', 'water_budget_z', 0., 30., \ 'Water|Budget|vertical|convergence|of|water|vapour|(+,|conv.' +\ ';|-,|div.)', 'mms-1', 'Reds'] elif varn == 'wbc' or varn == 'WBC' or varn == 'WBACCde': varvals = ['wbc', 'water_budget_c', 0., 30., \ 'Water|Budget|Cloud|water|species','mms-1', 'Reds'] elif varn == 'wbqvd' or varn == 'WBQVD' or varn == 'WBACQVDde': varvals = ['wbqvd', 'water_budget_qvd', 0., 30., \ 'Water|Budget|water|vapour|divergence', 'mms-1', 'Reds'] elif varn == 'wbqvblten' or varn == 'WBQVBLTEN' or varn == 'WBACQVBLTENde': varvals = ['wbqvblten', 'water_budget_qv_blten', 0., 30., \ 'Water|Budget|QV|tendency|due|to|pbl|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbqvcuten' or varn == 'WBQVCUTEN' or varn == 'WBACQVCUTENde': varvals = ['wbqvcuten', 'water_budget_qv_cuten', 0., 30., \ 'Water|Budget|QV|tendency|due|to|cu|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbqvshten' or varn == 'WBQVSHTEN' or varn == 'WBACQVSHTENde': varvals = ['wbqvshten', 'water_budget_qv_shten', 0., 30., \ 'Water|Budget|QV|tendency|due|to|shallow|cu|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbpr' or varn == 'WBP' or varn == 'WBACPde': varvals = ['wbpr', 'water_budget_pr', 0., 30., \ 'Water|Budget|recipitation', 'mms-1', 'Reds'] elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde': varvals = ['wbpw', 'water_budget_pw', 0., 30., \ 'Water|Budget|water|content', 'mms-1', 'Reds'] elif varn == 'wbcondt' or varn == 'WBCONDT' or varn == 'WBACCONDTde': varvals = ['wbcondt', 'water_budget_condt', 0., 30., \ 'Water|Budget|condensation|and|deposition', 'mms-1', 'Reds'] elif varn == 'wbqcm' or varn == 'WBQCM' or varn == 'WBACQCMde': varvals = ['wbqcm', 'water_budget_qcm', 0., 30., \ 'Water|Budget|hydrometeor|change|and|convergence', 'mms-1', 'Reds'] elif varn == 'wbsi' or varn == 'WBSI' or varn == 'WBACSIde': varvals = ['wbsi', 'water_budget_si', 0., 30., \ 'Water|Budget|hydrometeor|sink', 'mms-1', 'Reds'] elif varn == 'wbso' or varn == 'WBSO' or varn == 'WBACSOde': varvals = ['wbso', 'water_budget_so', 0., 30., \ 'Water|Budget|hydrometeor|source', 'mms-1', 'Reds'] # Water Budget accumulated elif varn == 'ccondac' or varn == 'ACCCOND': varvals = ['ccondac', 'cw_cond_ac', 0., 30., \ 'accumulated|cloud|water|condensation', 'mm', 'Reds'] elif varn == 'rac' or varn == 'ACQVAPOR': varvals = ['rac', 'ac_r', 0., 30., 'accumualted|water|wapor', 'mm', 'Blues'] elif varn == 'diabhac' or varn == 'ACDIABH': varvals = ['diabhac', 'diabh_ac', 0., 30., 'accumualted|diabatic|heating', \ 'K', 'Reds'] elif varn == 'wbpwac' or varn == 'WBACPW': varvals = ['wbpwac', 'water_budget_pw_ac', 0., 30., \ 'Water|Budget|accumulated|water|content', 'mm', 'Reds'] elif varn == 'wbfac' or varn == 'WBACF': varvals = ['wbfac', 'water_budget_hfcqv_ac', 0., 30., \ 'Water|Budget|accumulated|horizontal|convergence|of|water|vapour|(+,|' + \ 'conv.;|-,|div.)', 'mm', 'Reds'] elif varn == 'wbfcac' or varn == 'WBACFC': varvals = ['wbfcac', 'water_budget_fc_ac', 0., 30., \ 'Water|Budget|accumulated|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\ 'div.)', 'mm', 'Reds'] elif varn == 'wbfpac' or varn == 'WBACFP': varvals = ['wbfpac', 'water_budget_cfp_ac', 0., 30., \ 'Water|Budget|accumulated|horizontal|convergence|of|precipitation|(+,|' + \ 'conv.;|-,|div.)', 'mm', 'Reds'] elif varn == 'wbzac' or varn == 'WBACZ': varvals = ['wbzac', 'water_budget_z_ac', 0., 30., \ 'Water|Budget|accumulated|vertical|convergence|of|water|vapour|(+,|conv.' +\ ';|-,|div.)', 'mm', 'Reds'] elif varn == 'wbcac' or varn == 'WBACC': varvals = ['wbcac', 'water_budget_c_ac', 0., 30., \ 'Water|Budget|accumulated|Cloud|water|species','mm', 'Reds'] elif varn == 'wbqvdac' or varn == 'WBACQVD': varvals = ['wbqvdac', 'water_budget_qvd_ac', 0., 30., \ 'Water|Budget|accumulated|water|vapour|divergence', 'mm', 'Reds'] elif varn == 'wbqvbltenac' or varn == 'WBACQVBLTEN': varvals = ['wbqvbltenac', 'water_budget_qv_blten_ac', 0., 30., \ 'Water|Budget|accumulated|QV|tendency|due|to|pbl|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbqvcutenac' or varn == 'WBACQVCUTEN': varvals = ['wbqvcutenac', 'water_budget_qv_cuten_ac', 0., 30., \ 'Water|Budget|accumulated|QV|tendency|due|to|cu|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbqvshtenac' or varn == 'WBACQVSHTEN': varvals = ['wbqvshtenac', 'water_budget_qv_shten_ac', 0., 30., \ 'Water|Budget|accumulated|QV|tendency|due|to|shallow|cu|parameterization', \ 'kg kg-1 s-1', 'Reds'] elif varn == 'wbprac' or varn == 'WBACP': varvals = ['wbprac', 'water_budget_pr_ac', 0., 30., \ 'Water|Budget|accumulated|precipitation', 'mm', 'Reds'] elif varn == 'wbpwac' or varn == 'WBACPW': varvals = ['wbpwac', 'water_budget_pw_ac', 0., 30., \ 'Water|Budget|accumulated|water|content', 'mm', 'Reds'] elif varn == 'wbcondtac' or varn == 'WBACCONDT': varvals = ['wbcondtac', 'water_budget_condt_ac', 0., 30., \ 'Water|Budget|accumulated|condensation|and|deposition', 'mm', 'Reds'] elif varn == 'wbqcmac' or varn == 'WBACQCM': varvals = ['wbqcmac', 'water_budget_qcm_ac', 0., 30., \ 'Water|Budget|accumulated|hydrometeor|change|and|convergence', 'mm', 'Reds'] elif varn == 'wbsiac' or varn == 'WBACSI': varvals = ['wbsiac', 'water_budget_si_ac', 0., 30., \ 'Water|Budget|accumulated|hydrometeor|sink', 'mm', 'Reds'] elif varn == 'wbsoac' or varn == 'WBACSO': varvals = ['wbsoac', 'water_budget_so_ac', 0., 30., \ 'Water|Budget|accumulated|hydrometeor|source', 'mm', 'Reds'] elif varn == 'xtime' or varn == 'XTIME': varvals = ['xtime', 'time', 0., 1.e5, 'time', \ 'minutes|since|simulation|start', 'Reds'] elif varn == 'x' or varn == 'X': varvals = ['x', 'x', 0., 100., 'x', '-', 'Reds'] elif varn == 'y' or varn == 'Y': varvals = ['y', 'y', 0., 100., 'y', '-', 'Blues'] elif varn == 'z' or varn == 'Z': varvals = ['z', 'z', 0., 100., 'z', '-', 'Greens'] elif varn == 'zg' or varn == 'WRFght' or varn == 'Geopotential height' or \ varn == 'geop' or varn == 'LGEOP': varvals = ['zg', 'geopotential_height', 0., 80000., 'geopotential|height', \ 'm2s-2', 'rainbow'] elif varn == 'zmaxth' or varn == 'zmax_th' or varn == 'LZMAX_TH': varvals = ['zmaxth', 'thermal_plume_height', 0., 4000., \ 'maximum|thermals|plume|height', 'm', 'YlOrRd'] elif varn == 'zmla' or varn == 's_pblh' or varn == 'LS_PBLH': varvals = ['zmla', 'atmosphere_boundary_layer_thickness', 0., 2500., \ 'atmosphere|boundary|layer|thickness', 'm', 'Blues'] else: print errormsg print ' ' + fname + ": variable '" + varn + "' not defined !!!" quit(-1) return varvals def lonlat2D(lon,lat): """ Function to return lon, lat 2D matrices from any lon,lat matrix lon= matrix with longitude values lat= matrix with latitude values """ fname = 'lonlat2D' if len(lon.shape) != len(lat.shape): print errormsg print ' ' + fname + ': longitude values with shape:', lon.shape, \ 'is different than latitude values with shape:', lat.shape, '(dif. size) !!' quit(-1) if len(lon.shape) == 3: lonvv = lon[0,:,:] latvv = lat[0,:,:] elif len(lon.shape) == 2: lonvv = lon[:] latvv = lat[:] elif len(lon.shape) == 1: lonlatv = np.meshgrid(lon[:],lat[:]) lonvv = lonlatv[0] latvv = lonlatv[1] return lonvv, latvv def transform(vals, trans, dxv='0', dyv='0', dxt='0', dyt='0', dxl='0', \ dyl='0', dxtit='0', dytit='0', axxkind='fixpixel', axykind='fixpixel'): """ Function to transform the values and the axes vals= values to transform trans= '|' separated list of operations of transformation 'transpose': Transpose matrix of values (x-->y, y-->x) 'flip@[x/y]': Flip the given axis d[x/y]v= original values for the [x/y]-axis d[x/y]t= original ticks for the [x/y]-axis d[x/y]l= original tick-labels for the [x/y]-axis d[x/y]tit= original titles for the [x/y]-axis ax[x/y]kind= kind of axis reference 'fixpixel': pixel preserves size all along axis 'data': pixel follows data values (as in a georeferenced map) >>> mat = np.arange(15).reshape(3,5) >>> xvals = np.zeros((3,5), dtype=np.float) >>> yvals = np.zeros((3,5), dtype=np.float) >>> for i in range(5): >>> xvals[:,i] = -1.25 + i*0.5 >>> for j in range(3): >>> yvals[j,:] = 40.25 + j*0.5 >>> xticks = np.array([-1.25, -0.25, 0.25]) >>> yticks = np.array([40.25, 41.25]) >>> xticklab = ['1.25 W', '0.25 W', '0,25 E'] >>> yticklab = ['40.25 N', '41.25 N'] >>> xtit = 'Longitude' >>> ytit = 'Latitude' >>> nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT = transform(mat, 'transpose', \ xvals, yvals, xticks, yticks, xticklab, yticklab, xtit, ytit) [[ 0 5 10] [ 1 6 11] [ 2 7 12] [ 3 8 13] [ 4 9 14]] [[ 40.25 40.75 41.25] [ 40.25 40.75 41.25] [ 40.25 40.75 41.25] [ 40.25 40.75 41.25] [ 40.25 40.75 41.25]] [[-1.25 -1.25 -1.25] [-0.75 -0.75 -0.75] [-0.25 -0.25 -0.25] [ 0.25 0.25 0.25] [ 0.75 0.75 0.75]] [ 40.25 41.25] [-1.25 -0.25 0.25] ['41.25 N', '40.25 N'] ['0,25 E', '0.25 W', '1.25 W'] Latitude Longitude """ fname = 'transform' transforms = trans.split('|') Ntransforms = len(transforms) newvals = vals.copy() # values wdvx = False if dxv != '0': newdxv = dxv.copy() wdvx = True else: newdxv = None wdvy = False if dyv != '0': newdyv = dyv.copy() wdvy = True else: newdyv = None # ticks wdtx = False if dxt != '0': newdxt = dxt.copy() wdtx = True else: newdxt = None wdty = False if dyt != '0': newdyt = dyt.copy() wdty = True else: newdyt = None # label in ticks wdlx = False if dxl != '0': newdxl = list(dxl) wdlx = True else: newdxl = None wdly = False if dyl != '0': newdyl = list(dyl) wdly = True else: newdyl = None # axes' title wdTx = False if dxt != '0': newdxtit = str(dxtit) wdTx = True else: newdxtit = None wdTy = False if dyt != '0': newdytit = str(dytit) wdTy = True else: newdytit = None dx = newvals.shape[0] dy = newvals.shape[1] for transform in transforms: if transform == 'transpose': newvals = np.transpose(newvals) if wdvx and wdvy: copy = newdxv.copy() newdxv = np.transpose(newdyv) newdyv = np.transpose(copy) # else: # print errormsg # print ' ' + fname + ": it can not be '" + transform + "' without " +\ # ' both dimxvals:', dxv, ' and dimyvals:', dyv, ' !!' # quit(-1) if wdtx and wdtx: copy = newdxt.copy() newdxt = np.transpose(newdyt) newdyt = np.transpose(copy) # else: # print errormsg # print ' ' + fname + ": it can not be '" + transform + "' without " +\ # ' both dimxtickvals:', dxt, ' and dimytickvals:', dyt, ' !!' # quit(-1) if wdlx and wdly: copy = list(newdxl) newdxl = list(newdyl[::-1]) newdyl = list(copy[::-1]) # else: # print errormsg # print ' ' + fname + ": it can not be '" + transform + "' without " +\ # ' both dimxlabels:', dxl, ' and dimylabels:', dyl, ' !!' # quit(-1) if wdTx and wdTy: copy = str(newdxtit) newdxtit = str(newdytit) newdytit = str(copy) # else: # print errormsg # print ' ' + fname + ": it can not be '" + transform + "' without " +\ # ' both dimx title:', dxtit, ' and dimytitle:', dytit, ' !!' # quit(-1) elif transform[0:4] == 'flip': if transform.find('@') == -1: print errormsg print ' ' + fname + ": no ax given for '" + transform + "' !!" quit(-1) flip = transform.split('@')[1] if flip == 'x': newvals = newvals[...,::-1] if axxkind == 'fixpixel': if wdtx: # No pixel deformation newdxt = 1.-newdxt elif axxkind == 'data': if wdvx: if len(newdxv.shape) > 1: newdxv = newdxv[...,::-1] else: newdxv = newdxv[::-1] if wdtx: newdxt = newdxt[::-1] # if wdlx: # newdxl = newdxl[::-1] elif flip == 'y': newvals = newvals[...,::-1,:] # No pixel deformation if axykind == 'fixpixel': if wdty: # No pixel deformation newdyt = 1.-newdyt elif axykind == 'data': if wdvy: if len(newdyv.shape) > 1: newdyv = newdyv[...,::-1] else: newdyv = newdyv[::-1] if wdty: newdyt = newdyt[::-1] # if wdly: # newdyl = newdyl[::-1] elif flip == 'z': newvals = newvals[...,::-1,:,:] else: print errormsg print ' '+ fname + ": transformation '" + transform + "' not ready!!" quit(-1) return newvals, newdxv, newdyv, newdxt, newdyt, newdxl, newdyl, newdxtit, newdytit #mat = np.arange(15).reshape(3,5) #xvals = np.zeros((3,5), dtype=np.float) #yvals = np.zeros((3,5), dtype=np.float) #for i in range(5): # xvals[:,i] = -1.25 + i*0.5 #for j in range(3): # yvals[j,:] = 40.25 + j*0.5 #xticks = np.array([-1.25, -0.25, 0.25]) #yticks = np.array([40.25, 41.25]) #xticklab = ['1.25 W', '0.25 W', '0,25 E'] #yticklab = ['40.25 N', '41.25 N'] #xtit = 'Longitude' #ytit = 'Latitude' #nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT = transform(mat, 'flip@y', \ # xvals, yvals, xticks, yticks, xticklab, yticklab, xtit, ytit) # #print nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT #quit() def legend_values(legstring, char): """ Function to determine the values for the legend legstring= char separated list of [locleg][char][fontsize]: [locleg]: location of the legend (0, 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' [fontsize]: font size for the legend ('auto' for 12) char= separation character """ fname = 'legend_values' locleg = int(legstring.split('|')[0]) fontsize0 = legstring.split('|')[1] if fontsize0 == 'auto': fontsize = 12 else: fontsize = int(fontsize0) return locleg, fontsize ####### ####### ####### ####### ####### ####### ####### ####### ####### ####### def check_colorBar(cbarn): """ Check if the given colorbar exists in matplotlib """ fname = 'check_colorBar' # FROM: http://matplotlib.org/1.2.1/examples/pylab_examples/show_colormaps.html # Get a list of the colormaps in matplotlib. Ignore the ones that end with # '_r' because these are simply reversed versions of ones that don't end # with '_r' maps = sorted(m for m in plt.cm.datad if not m.endswith("_r")) nmaps = len(maps) + 1 # Possible color bars # colorbars = ['binary', 'Blues', 'BuGn', 'BuPu', 'gist_yarg', 'GnBu', 'Greens', \ # 'Greys', 'Oranges', 'OrRd', 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu', \ # 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'afmhot', 'autumn', 'bone', \ # 'cool', 'copper', 'gist_gray', 'gist_heat', 'gray', 'hot', 'pink', 'spring', \ # 'summer', 'winter', 'BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr', 'RdBu', \ # 'RdGy', 'RdYlBu', 'RdYlGn', 'seismic', 'Accent', 'Dark2', 'hsv', 'Paired', \ # 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3', 'spectral', 'gist_earth', \ # 'gist_ncar', 'gist_rainbow', 'gist_stern', 'jet', 'brg', 'CMRmap', 'cubehelix',\ # 'gnuplot', 'gnuplot2', 'ocean', 'rainbow', 'terrain', 'flag', 'prism'] if not searchInlist(maps,cbarn): print warnmsg print ' ' + fname + ' color bar: "' + cbarn + '" does not exist !!' print " the standard 'rainbow' will be use instead !!" print ' colorbars availables with this version:', maps return def colorbar_vals(colorbarvs,splitchar): """ Function to provide the colorbar values for a figure colorbarvs= [splitchar] separated list of different colorbar values [cbn]: name of the color bar to use [fmtcb]: format of the tick labels in the colorbar ('C'-like) ('auto' for '%6g') [cbor]: orientation of the color bar ('vertical'/'horizontal') ('auto' for vertical) """ fname = 'colorbar_vals' colbarn = colorbarvs.split(splitchar)[0] fmtcolbar0 = colorbarvs.split(splitchar)[1] colbaror0 = colorbarvs.split(splitchar)[2] if fmtcolbar0 == 'auto': fmtcolbar = '%6g' else: fmtcolbar = fmtcolbar0 if colbaror0 == 'auto': colbaror = 'vertical' else: colbaror = colbaror0 # check colorbar existence check_colorBar(colbarn) return colbarn, fmtcolbar, colbaror def units_lunits(u): """ Fucntion to provide LaTeX equivalences from a given units u= units to transform >>> units_lunits('kgkg-1') '$kgkg^{-1}$' """ fname = 'units_lunits' if u == 'h': print fname + '_____________________________________________________________' print units_lunits.__doc__ quit() # Units which does not change same = ['1', 'category', 'day', 'deg', 'degree', 'degrees', 'degrees East', \ 'degrees Nord', 'degrees North', 'g', 'gpm', 'hour', 'hPa', 'K', 'Km', 'kg', \ 'km', 'm', 'minute', 'mm', 'month', 'Pa', 's', 'second', 'um', 'year', '-'] if searchInlist(same,u): lu = '$' + u + '$' elif len(u.split(' ')) > 1 and u.split(' ')[1] == 'since': uparts = u.split(' ') ip=0 for up in uparts: if ip == 0: lu = '$' + up else: lu = lu + '\ ' + up ip=ip+1 lu = lu + '$' else: if u == '': lu='-' elif u == 'C': lu='$^{\circ}C$' elif u == 'days': lu='$day$' elif u == 'degrees_East': lu='$degrees\ East$' elif u == 'degrees_east': lu='$degrees\ East$' elif u == 'degree_east': lu='$degrees\ East$' elif u == 'degrees longitude': lu='$degrees\ East$' elif u == 'degrees latitude': lu='$degrees\ North$' elif u == 'degrees_North': lu='$degrees\ North$' elif u == 'degrees_north': lu='$degrees\ North$' elif u == 'degree_north': lu='$degrees\ North$' elif u == 'deg C': lu='$^{\circ}C$' elif u == 'degC': lu='$^{\circ}C$' elif u == 'deg K': lu='$K$' elif u == 'degK': lu='$K$' elif u == 'g/g': lu='$gg^{-1}$' elif u == 'hours': lu='$hour$' elif u == 'J/kg': lu='$Jkg^{-1}$' elif u == 'Jkg-1': lu='$Jkg^{-1}$' elif u == 'K/h': lu='$Kh^{-1}$' elif u == 'Kh-1': lu='$Kh^{-1}$' elif u == 'K/m': lu='$Km^{-1}$' elif u == 'Km-1': lu='$Km^{-1}$' elif u == 'Km-3': lu='$Km^{-3}$' elif u == 'K/s': lu='$Ks^{-1}$' elif u == 'Ks-1': lu='$Ks^{-1}$' elif u == 'K s-1': lu='$Ks^{-1}$' elif u == 'K/day': lu='$Kday^{-1}$' elif u == 'Kday-1': lu='$Kday^{-1}$' elif u == 'K day-1': lu='$Kday^{-1}$' elif u == 'kg/kg': lu='$kgkg^{-1}$' elif u == 'kgkg-1': lu='$kgkg^{-1}$' elif u == 'kg kg-1': lu='$kgkg^{-1}$' elif u == '(kg/kg)/s': lu='$kgkg^{-1}s^{-1}$' elif u == 'kgkg-1s-1': lu='$kgkg^{-1}s^{-1}$' elif u == 'kg kg-1 s-1': lu='$kgkg^{-1}s^{-1}$' elif u == 'kg/m2': lu='$kgm^{-2}$' elif u == 'kgm-2': lu='$kgm^{-2}$' elif u == 'kg m-2': lu='$kgm^{-2}$' elif u == 'Kg m-2': lu='$kgm^{-2}$' elif u == 'kg/m2/s': lu='$kgm^{-2}s^{-1}$' elif u == 'kg/m^2/s': lu='$kgm^{-2}s^{-1}$' elif u == 'kg/(m2*s)': lu='$kgm^{-2}s^{-1}$' elif u == 'kg/(s*m2)': lu='$kgm^{-2}s^{-1}$' elif u == 'kgm-2s-1': lu='$kgm^{-2}s^{-1}$' elif u == 'kg m-2 s-1': lu='$kgm^{-2}s^{-1}$' elif u == '1/m': lu='$m^{-1}$' elif u == 'm-1': lu='$m^{-1}$' elif u == 'm^2': lu='$m^{2}$' elif u == 'm2/s': lu='$m2s^{-1}$' elif u == 'm2s-1': lu='$m2s^{-1}$' elif u == 'm2/s2': lu='$m2s^{-2}$' elif u == 'm3/s': lu='$m^{3}s^{-1}$' elif u == 'm^3/s': lu='$m^{3}s^{-1}$' elif u == 'm/s': lu='$ms^{-1}$' elif u == 'mmh-3': lu='$mmh^{-3}$' elif u == 'ms-1': lu='$ms^{-1}$' elif u == 'm s-1': lu='$ms^{-1}$' elif u == 'm/s2': lu='$ms^{-2}$' elif u == 'ms-2': lu='$ms^{-2}$' elif u == 'minutes': lu='$minute$' elif u == 'meters MSL': lu='$m$' elif u == 'none': lu='-' elif u == 'Pa/s': lu='$Pas^{-1}$' elif u == 'Pas-1': lu='$Pas^{-1}$' elif u == 'W m-2': lu='$Wm^{-2}$' elif u == 'Wm-2': lu='$Wm^{-2}$' elif u == 'W/m2': lu='$Wm^{-2}$' elif u == '1/s': lu='$s^{-1}$' elif u == 's-1': lu='$s^{-1}$' elif u == 'seconds': lu='$second$' elif u == '%': lu='\%' elif u == '?': lu='-' else: print errormsg print ' ' + fname + ': units "' + u + '" not ready!!!!' quit(-1) return lu def ASCII_LaTeX(ln): """ Function to transform from an ASCII line to LaTeX codification >>> ASCII_LaTeX('Laboratoire de Météorologie Dynamique però Hovmöller') Laboratoire de M\'et\'eorologie Dynamique per\`o Hovm\"oller """ fname='ASCII_LaTeX' if ln == 'h': print fname + '_____________________________________________________________' print ASCII_LaTeX.__doc__ quit() newln = ln.replace('\\', '\\textbackslash') newln = newln.replace('á', "\\'a") newln = newln.replace('é', "\\'e") newln = newln.replace('í', "\\'i") newln = newln.replace('ó', "\\'o") newln = newln.replace('ú', "\\'u") newln = newln.replace('à', "\\`a") newln = newln.replace('è', "\\`e") newln = newln.replace('ì', "\\`i") newln = newln.replace('ò', "\\`o") newln = newln.replace('ù', "\\`u") newln = newln.replace('â', "\\^a") newln = newln.replace('ê', "\\^e") newln = newln.replace('î', "\\^i") newln = newln.replace('ô', "\\^o") newln = newln.replace('û', "\\^u") newln = newln.replace('ä', '\\"a') newln = newln.replace('ë', '\\"e') newln = newln.replace('ï', '\\"i') newln = newln.replace('ö', '\\"o') newln = newln.replace('ü', '\\"u') newln = newln.replace('ç', '\c{c}') newln = newln.replace('ñ', '\~{n}') newln = newln.replace('Á', "\\'A") newln = newln.replace('É', "\\'E") newln = newln.replace('Í', "\\'I") newln = newln.replace('Ó', "\\'O") newln = newln.replace('Ú', "\\'U") newln = newln.replace('À', "\\`A") newln = newln.replace('È', "\\`E") newln = newln.replace('Ì', "\\`I") newln = newln.replace('Ò', "\\`O") newln = newln.replace('Ù', "\\`U") newln = newln.replace('Â', "\\^A") newln = newln.replace('Ê', "\\^E") newln = newln.replace('Î', "\\^I") newln = newln.replace('Ô', "\\^O") newln = newln.replace('Û', "\\^U") newln = newln.replace('Ä', '\\"A') newln = newln.replace('Ë', '\\"E') newln = newln.replace('Ï', '\\"I') newln = newln.replace('Ö', '\\"O') newln = newln.replace('Ü', '\\"U') newln = newln.replace('Ç', '\\c{C}') newln = newln.replace('Ñ', '\\~{N}') newln = newln.replace('¡', '!`') newln = newln.replace('¿', '¿`') newln = newln.replace('%', '\\%') newln = newln.replace('#', '\\#') newln = newln.replace('&', '\\&') newln = newln.replace('$', '\\$') newln = newln.replace('_', '\\_') newln = newln.replace('·', '\\textperiodcentered') newln = newln.replace('<', '$<$') newln = newln.replace('>', '$>$') newln = newln.replace('', '*') # newln = newln.replace('º', '$^{\\circ}$') newln = newln.replace('ª', '$^{a}$') newln = newln.replace('º', '$^{o}$') newln = newln.replace('°', '$^{\\circ}$') newln = newln.replace('\n', '\\\\\n') newln = newln.replace('\t', '\\medskip') return newln def DegGradSec_deg(grad,deg,sec): """ Function to transform from a coordinate in grad deg sec to degrees (decimal) >>> DegGradSec_deg(39.,49.,26.) 39.8238888889 """ fname = 'DegGradSec_deg' if grad == 'h': print fname + '_____________________________________________________________' print DegGradSec_deg.__doc__ quit() deg = grad + deg/60. + sec/3600. return deg def intT2dt(intT,tu): """ Function to provide an 'timedelta' object from a given interval value intT= interval value tu= interval units, [tu]= 'd': day, 'w': week, 'h': hour, 'i': minute, 's': second, 'l': milisecond >>> intT2dt(3.5,'s') 0:00:03.500000 >>> intT2dt(3.5,'w') 24 days, 12:00:00 """ import datetime as dt fname = 'intT2dt' if tu == 'w': dtv = dt.timedelta(weeks=np.float(intT)) elif tu == 'd': dtv = dt.timedelta(days=np.float(intT)) elif tu == 'h': dtv = dt.timedelta(hours=np.float(intT)) elif tu == 'i': dtv = dt.timedelta(minutes=np.float(intT)) elif tu == 's': dtv = dt.timedelta(seconds=np.float(intT)) elif tu == 'l': dtv = dt.timedelta(milliseconds=np.float(intT)) else: print errormsg print ' ' + fname + ': time units "' + tu + '" not ready!!!!' quit(-1) return dtv def lonlat_values(mapfile,lonvn,latvn): """ Function to obtain the lon/lat matrices from a given netCDF file lonlat_values(mapfile,lonvn,latvn) [mapfile]= netCDF file name [lonvn]= variable name with the longitudes [latvn]= variable name with the latitudes """ fname = 'lonlat_values' if mapfile == 'h': print fname + '_____________________________________________________________' print lonlat_values.__doc__ quit() if not os.path.isfile(mapfile): print errormsg print ' ' + fname + ": map file '" + mapfile + "' does not exist !!" quit(-1) ncobj = NetCDFFile(mapfile, 'r') lonobj = ncobj.variables[lonvn] latobj = ncobj.variables[latvn] if len(lonobj.shape) == 3: lonv = lonobj[0,:,:] latv = latobj[0,:,:] elif len(lonobj.shape) == 2: lonv = lonobj[:,:] latv = latobj[:,:] elif len(lonobj.shape) == 1: lon0 = lonobj[:] lat0 = latobj[:] lonv = np.zeros( (len(lat0),len(lon0)), dtype=np.float ) latv = np.zeros( (len(lat0),len(lon0)), dtype=np.float ) for iy in range(len(lat0)): lonv[iy,:] = lon0 for ix in range(len(lon0)): latv[:,ix] = lat0 else: print errormsg print ' ' + fname + ': lon/lat variables shape:',lonobj.shape,'not ready!!' quit(-1) return lonv, latv def date_CFtime(ind,refd,tunits): """ Function to transform from a given date object a CF-convention time ind= date object to transform refd= reference date tunits= units for time >>> date_CFtime(dt.datetime(1976,02,17,08,30,00), dt.datetime(1949,12,01,00,00,00), 'seconds') 827224200.0 """ import datetime as dt fname = 'date_CFtime' dt = ind - refd if tunits == 'weeks': value = dt.days/7. + dt.seconds/(3600.*24.*7.) elif tunits == 'days': value = dt.days + dt.seconds/(3600.*24.) elif tunits == 'hours': value = dt.days*24. + dt.seconds/(3600.) elif tunits == 'minutes': value = dt.days*24.*60. + dt.seconds/(60.) elif tunits == 'seconds': value = dt.days*24.*3600. + dt.seconds elif tunits == 'milliseconds': value = dt.days*24.*3600.*1000. + dt.seconds*1000. else: print errormsg print ' ' + fname + ': reference time units "' + trefu + '" not ready!!!!' quit(-1) return value def pot_values(values, uvals): """ Function to modify a seies of values by their potency of 10 pot_values(values, uvals) values= values to modify uvals= units of the values >>> vals = np.sin(np.arange(20)*np.pi/5.+0.01)*10.e-5 >>> pot_values(vals,'ms-1') (array([ 0.00000000e+00, 5.87785252e-01, 9.51056516e-01, 9.51056516e-01, 5.87785252e-01, 1.22464680e-16, -5.87785252e-01, -9.51056516e-01, -9.51056516e-01, -5.87785252e-01, -2.44929360e-16, 5.87785252e-01, 9.51056516e-01, 9.51056516e-01, 5.87785252e-01, 3.67394040e-16, -5.87785252e-01, -9.51056516e-01, -9.51056516e-01, -5.87785252e-01]), -4, 'x10e-4 ms-1', 'x10e-4') """ fname = 'pot_values' if np.min(values) != 0.: potmin = int( np.log10( np.abs(np.min(values)) ) ) else: potmin = 0 if np.max(values) != 0.: potmax = int( np.log10( np.abs(np.max(values)) ) ) else: potmax = 0 if potmin * potmax > 9: potval = -np.min([np.abs(potmin), np.abs(potmax)]) * np.abs(potmin) / potmin newvalues = values*10.**potval potvalue = - potval potS = 'x10e' + str(potvalue) newunits = potS + ' ' + uvals else: newvalues = values potvalue = None potS = '' newunits = uvals return newvalues, potvalue, newunits, potS def CFtimes_plot(timev,units,kind,tfmt): """ Function to provide a list of string values from a CF time values in order to use them in a plot, according to the series of characteristics. String outputs will be suited to the 'human-like' output timev= time values (real values) units= units string according to CF conventions ([tunits] since [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]]) 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 >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'Nval,5',"%Y/%m/%d %H:%M:%S") 0.0 1979/12/01 00:00:00 24.75 1979/12/02 00:45:00 49.5 1979/12/03 01:30:00 74.25 1979/12/04 02:15:00 99.0 1979/12/05 03:00:00 >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'exct,2,d',"%Y/%m/%d %H:%M:%S") 0.0 1979/12/01 00:00:00 48.0 1979/12/03 00:00:00 96.0 1979/12/05 00:00:00 144.0 1979/12/07 00:00:00 """ import datetime as dt # Seconds between 0001 and 1901 Jan - 01 secs0001_1901=59958144000. fname = 'CFtimes_plot' if timev == 'h': print fname + '_____________________________________________________________' print CFtimes_plot.__doc__ quit() secsYear = 365.*24.*3600. secsWeek = 7.*24.*3600. secsDay = 24.*3600. secsHour = 3600. secsMinute = 60. secsMilisecond = 1./1000. secsMicrosecond = 1./1000000. # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] ## trefT = units.find(':') txtunits = units.split(' ') Ntxtunits = len(txtunits) if Ntxtunits == 3: Srefdate = txtunits[Ntxtunits - 1] else: Srefdate = txtunits[Ntxtunits - 2] if not trefT == -1: # print ' ' + fname + ': refdate with time!' if Ntxtunits == 3: refdate = datetimeStr_datetime(Srefdate) else: refdate = datetimeStr_datetime(Srefdate + '_' + txtunits[Ntxtunits - 1]) else: refdate = datetimeStr_datetime(Srefdate + '_00:00:00') trefunits=units.split(' ')[0] if trefunits == 'weeks': trefu = 'w' elif trefunits == 'days': trefu = 'd' elif trefunits == 'hours': trefu = 'h' elif trefunits == 'minutes': trefu = 'm' elif trefunits == 'seconds': trefu = 's' elif trefunits == 'milliseconds': trefu = 'l' else: print errormsg print ' ' + fname + ': reference time units "' + trefu + '" not ready!!!!' quit(-1) okind=kind.split(',')[0] dtv = len(timev) if refdate.year == 1: print warnmsg print ' ' + fname + ': changing reference date: ',refdate, \ 'to 1901-01-01_00:00:00 !!!' refdate = datetimeStr_datetime('1901-01-01_00:00:00') if trefu == 'w': timev = timev - secs0001_1901/(7.*24.*3600.) if trefu == 'd': timev = timev - secs0001_1901/(24.*3600.) if trefu == 'h': timev = timev - secs0001_1901/(3600.) if trefu == 'm': timev = timev - secs0001_1901/(60.) if trefu == 's': timev = timev - secs0001_1901 if trefu == 'l': timev = timev - secs0001_1901*1000. # if refdate.year < 1901: # print warnmsg # print ' ' + fname + ': changing reference date: ',refdate, \ # 'to 1901-01-01_00:00:00 !!!' # diffrefs = dt.datetime(1901,1,1,0,0,0) - refdate # if gen.searchInlist(dir(diffrefs), 'total_seconds'): # secsref_1901 = diffrefs.total_seconds() # else: # secsref_1901 = diffrefs.days()*3600*24. + diffrefs.seconds()*1. # refdate = datetimeStr_datetime('1901-01-01_00:00:00') # print 'time 0 :', timev[0] # if trefu == 'w': timev = timev - secsref_1901/(7.*24.*3600.) # if trefu == 'd': timev = timev - secsref_1901/(24.*3600.) # if trefu == 'h': timev = timev - secsref_1901/(3600.) # if trefu == 'm': timev = timev - secsref_1901/(60.) # if trefu == 's': timev = timev - secsref_1901 # if trefu == 'l': timev = timev - secsref_1901*1000. # print 'time 0 after:', timev[0] firstT = timev[0] lastT = timev[dtv-1] # First and last times as datetime objects firstTdt = timeref_datetime(refdate, firstT, trefunits) lastTdt = timeref_datetime(refdate, lastT, trefunits) # First and last times as [year, mon, day, hour, minut, second] vectors firstTvec = np.zeros((6), dtype= np.float) lastTvec = np.zeros((6), dtype= np.float) chTvec = np.zeros((6), dtype= bool) firstTvec = np.array([firstTdt.year, firstTdt.month, firstTdt.day, firstTdt.hour,\ firstTdt.minute, firstTdt.second]) lastTvec = np.array([lastTdt.year, lastTdt.month, lastTdt.day, lastTdt.hour, \ lastTdt.minute, lastTdt.second]) chdate= lastTvec - firstTvec chTvec = np.where (chdate != 0., True, False) TOTdt = lastTdt - firstTdt TOTdtsecs = TOTdt.days*secsDay + TOTdt.seconds + TOTdt.microseconds*secsMicrosecond timeout = [] if okind == 'Nval': nvalues = int(kind.split(',')[1]) intervT = (lastT - firstT)/(nvalues-1) dtintervT = intT2dt(intervT, trefu) for it in range(nvalues): timeout.append(firstTdt + dtintervT*it) elif okind == 'exct': Nunits = int(kind.split(',')[1]) tu = kind.split(',')[2] # Generic incremental dt [seconds] according to all the possibilities ['c', 'y', 'm', # 'w', 'd', 'h', 'i', 's', 'l'], some of them approximated (because they are not # already necessary!) basedt = np.zeros((9), dtype=np.float) basedt[0] = (365.*100. + 25.)*24.*3600. basedt[1] = secsYear basedt[2] = 31.*24.*3600. basedt[3] = secsWeek basedt[4] = secsDay basedt[5] = secsHour basedt[6] = secsMinute basedt[7] = 1. basedt[8] = secsMilisecond # Increment according to the units of the CF dates if trefunits == 'weeks': basedt = basedt/(secsWeek) elif trefunits == 'days': basedt = basedt/(secsDay) elif trefunits == 'hours': basedt = basedt/(secsHour) elif trefunits == 'minutes': basedt = basedt/(secsMinute) elif trefunits == 'seconds': basedt = basedt elif trefunits == 'milliseconds': basedt = basedt*secsMilisecond if tu == 'c': ti = firstTvec[0] tf = lastTvec[0] centi = firstTvec[0] / 100 datev = firstTdt while datev < lastTdt: yr = datev.year + Nunits*100 mon = datev.month datev = dt.datetime(yr, mon, 1, 0, 0, 0) timeout.append(datev) elif tu == 'y': ti = firstTvec[0] tf = lastTvec[0] year = firstTvec[0] TOTsteps = int(TOTdtsecs/(Nunits*31*secsDay)) + 1 if np.sum(firstTvec[1:5]) > 0: firstTdt = dt.datetime(year, 1, 1, 0, 0, 0) datev = dt.datetime(year, 1, 1, 0, 0, 0) datev = firstTdt timeout.append(datev) yr = datev.year while yr + Nunits <= tf : yr = datev.year + Nunits mon = datev.month datev = dt.datetime(yr, mon, 1, 0, 0, 0) timeout.append(datev) elif tu == 'm': ti = firstTvec[1] tf = lastTvec[1] yr = firstTvec[0] mon = firstTvec[1] TOTsteps = int(TOTdtsecs/(Nunits*31*secsDay)) + 1 datev = firstTdt while datev < lastTdt: mon = datev.month + Nunits if mon > 12: yr = yr + 1 mon = 1 datev = dt.datetime(yr, mon, 1, 0, 0, 0) timeout.append(datev) elif tu == 'w': datev=firstTdt it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(days=7*Nunits*it) timeout.append(datev) it = it + 1 elif tu == 'd': # datev=firstTdt yr = firstTvec[0] mon = firstTvec[1] day = firstTvec[2] Iunits = np.mod(day,Nunits) if np.sum(firstTvec[2:5]) > 0: firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0) datev = dt.datetime(yr, mon, day+1, 0, 0, 0) elif Iunits != 0: nNunits = int(day/Nunits) firstTdt = dt.datetime(yr, mon, nNunits*Nunits, 0, 0, 0) datev = dt.datetime(yr, mon, nNunits*Nunits, 0, 0, 0) else: firstTdt = dt.datetime(yr, mon, day, 0, 0, 0) datev = dt.datetime(yr, mon, day, 0, 0, 0) it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(days=Nunits*it) timeout.append(datev) it = it + 1 elif tu == 'h': datev=firstTdt yr = firstTvec[0] mon = firstTvec[1] day = firstTvec[2] hour = firstTvec[3] Iunits = np.mod(hour,Nunits) if np.sum(firstTvec[4:5]) > 0 or Iunits != 0: tadvance = 2*Nunits if tadvance >= 24: firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0) datev = dt.datetime(yr, mon, day+1, 0, 0, 0) else: nNunits = int(hour/Nunits) firstTdt = dt.datetime(yr, mon, day, nNunits*Nunits, 0, 0) datev = dt.datetime(yr, mon, day, nNunits*Nunits, 0, 0) else: firstTdt = dt.datetime(yr, mon, day, hour, 0, 0) datev = dt.datetime(yr, mon, day, hour, 0, 0) it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(seconds=Nunits*3600*it) timeout.append(datev) it = it + 1 elif tu == 'i': datev=firstTdt yr = firstTvec[0] mon = firstTvec[1] day = firstTvec[2] hour = firstTvec[3] minu = firstTvec[4] Iunits = np.mod(minu,Nunits) if firstTvec[5] > 0 or Iunits != 0: tadvance = 2*Nunits if tadvance >= 60: firstTdt = dt.datetime(yr, mon, day, hour, 0, 0) datev = dt.datetime(yr, mon, day, hour, 0, 0) else: nNunits = int(minu/Nunits) firstTdt = dt.datetime(yr, mon, day, hour, nNunits*Nunits, 0) datev = dt.datetime(yr, mon, day, hour, nNunits*Nunits, 0) else: firstTdt = dt.datetime(yr, mon, day, hour, minu, 0) datev = dt.datetime(yr, mon, day, hour, minu, 0) it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(seconds=Nunits*60*it) timeout.append(datev) it = it + 1 elif tu == 's': datev=firstTdt yr = firstTvec[0] mon = firstTvec[1] day = firstTvec[2] hour = firstTvec[3] min = firstTvec[4] secu = firstTvec[5] Iunits = np.mod(secu,Nunits) if firstTvec[5] > 0 or Iunits != 0: tadvance = 2*Nunits if tadvance >= 60: firstTdt = dt.datetime(yr, mon, day, hour, min, 0) datev = dt.datetime(yr, mon, day, hour, min, 0) else: nNunits = int(minu/Nunits) firstTdt = dt.datetime(yr, mon, day, hour, min, nNunits*Nunits) datev = dt.datetime(yr, mon, day, hour, min, nNunits*Nunits) else: firstTdt = dt.datetime(yr, mon, day, hour, min, secu) datev = dt.datetime(yr, mon, day, hour, min, secu) it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(seconds=Nunits*it) timeout.append(datev) it = it + 1 elif tu == 'l': datev=firstTdt it=0 while datev <= lastTdt: datev = firstTdt + dt.timedelta(seconds=Nunits*it/1000.) timeout.append(datev) it = it + 1 else: print errormsg print ' ' + fname + ': exact units "' + tu + '" not ready!!!!!' quit(-1) else: print errormsg print ' ' + fname + ': output kind "' + okind + '" not ready!!!!' quit(-1) dtout = len(timeout) timeoutS = [] timeoutv = np.zeros((dtout), dtype=np.float) for it in range(dtout): timeoutS.append(timeout[it].strftime(tfmt)) timeoutv[it] = date_CFtime(timeout[it], refdate, trefunits) # print it,':',timeoutv[it], timeoutS[it] if len(timeoutv) < 1 or len(timeoutS) < 1: print errormsg print ' ' + fname + ': no time values are generated!' print ' values passed:',timev print ' units:',units print ' desired kind:',kind print ' format:',tfmt print ' function values ___ __ _' print ' reference date:',refdate print ' first date:',firstTdt print ' last date:',lastTdt print ' increment:',basedt,trefunits quit(-1) return timeoutv, timeoutS def color_lines(Nlines): """ Function to provide a color list to plot lines color_lines(Nlines) Nlines= number of lines """ fname = 'color_lines' colors = ['r', 'b', 'g', 'p', 'g'] colorv = [] colorv.append('k') for icol in range(Nlines): colorv.append(colors[icol]) return colorv def output_kind(kindf, namef, close): """ Function to generate the output of the figure kindf= kind of the output null: show in screen [jpg/pdf/png/ps]: standard output types namef= name of the figure (without extension) close= if the graph has to be close or not [True/False] """ fname = 'output_kind' if kindf == 'h': print fname + '_____________________________________________________________' print output_kind.__doc__ quit() if close: if kindf == 'null': print 'showing figure...' plt.show() elif kindf == 'gif': plt.savefig(namef + ".gif") if close: print "Successfully generation of figure '" + namef + ".jpg' !!!" elif kindf == 'jpg': plt.savefig(namef + ".jpg") if close: print "Successfully generation of figure '" + namef + ".jpg' !!!" elif kindf == 'pdf': plt.savefig(namef + ".pdf") if close: print "Successfully generation of figure '" + namef + ".pdf' !!!" elif kindf == 'png': plt.savefig(namef + ".png") if close: print "Successfully generation of figure '" + namef + ".png' !!!" elif kindf == 'ps': plt.savefig(namef + ".ps") if close: print "Successfully generation of figure '" + namef + ".ps' !!!" else: print errormsg print ' ' + fname + ' output format: "' + kindf + '" not ready !!' print errormsg quit(-1) plt.close() return def check_arguments(funcname,args,expectargs,char): """ Function to check the number of arguments if they are coincident check_arguments(funcname,Nargs,args,char) funcname= name of the function/program to check args= passed arguments expectargs= expected arguments char= character used to split the arguments """ fname = 'check_arguments' Nvals = len(args.split(char)) Nargs = len(expectargs.split(char)) if Nvals != Nargs: print errormsg print ' ' + fname + ': wrong number of arguments:',Nvals," passed to '", \ funcname, "' which requires:",Nargs,'!!' print ' # given expected _______' Nmax = np.max([Nvals, Nargs]) for ia in range(Nmax): if ia > Nvals-1: aval = ' ' else: aval = args.split(char)[ia] if ia > Nargs-1: aexp = ' ' else: aexp = expectargs.split(char)[ia] print ' ', ia, aval, aexp quit(-1) return def Str_Bool(val): """ Function to transform from a String value to a boolean one >>> Str_Bool('True') True >>> Str_Bool('0') False >>> Str_Bool('no') False """ fname = 'Str_Bool' if val == 'True' or val == 'true' or val == '1' or val == 'yes': boolv = True elif val == 'False' or val == 'false' or val == '0' or val== 'no': boolv = False else: print errormsg print ' ' + fname + ": value '" + val + "' not ready!!" quit(-1) return boolv def graphic_range(vrng,varsV): """ Function to provide the ranges of a figure following different options vrng: list of two values: maximum and minimum for the range of the figure '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) """ fname = 'graphic_range' rangegrph = np.zeros((2), dtype=np.float) # Changing limits of the colors if type(vrng[0]) == type('S') and vrng[0][0:1] == 'S': if vrng[0] == 'Srange': rangegrph[0] = np.min(varsV) rangegrph[1] = np.max(varsV) elif vrng[0][0:11] == 'Saroundmean': meanv = np.mean(varsV) permean = np.float(vrng[0].split('@')[1]) minv = np.min(varsV)*permean maxv = np.max(varsV)*permean minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)]) rangegrph[0] = meanv-minextrm rangegrph[1] = meanv+minextrm elif vrng[0][0:13] == 'Saroundminmax': permean = np.float(vrng[0].split('@')[1]) minv = np.min(varsV)*permean maxv = np.max(varsV)*permean rangegrph[0] = minv rangegrph[1] = maxv elif vrng[0][0:17] == 'Saroundpercentile': medianv = np.median(varsV) valper = np.float(vrng[0].split('@')[1]) minv = np.percentile(varsV, valper) maxv = np.percentile(varsV, 100.-valper) minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) rangegrph[0] = medianv-minextrm rangegrph[1] = medianv+minextrm elif vrng[0][0:5] == 'Smean': meanv = np.mean(varsV) permean = np.float(vrng[0].split('@')[1]) minv = np.min(varsV)*permean maxv = np.max(varsV)*permean minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)]) rangegrph[0] = -minextrm rangegrph[1] = minextrm elif vrng[0][0:7] == 'Smedian': medianv = np.median(varsV) permedian = np.float(vrng[0].split('@')[1]) minv = np.min(varsV)*permedian maxv = np.max(varsV)*permedian minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) rangegrph[0] = -minextrm rangegrph[1] = minextrm elif vrng[0][0:11] == 'Spercentile': medianv = np.median(varsV) valper = np.float(vrng[0].split('@')[1]) minv = np.percentile(varsV, valper) maxv = np.percentile(varsV, 100.-valper) minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) rangegrph[0] = -minextrm rangegrph[1] = minextrm else: print errormsg print ' ' + fname + ": range '" + vrng[0] + "' not ready!!!" quit(-1) print ' ' + fname + ': modified shadow min,max:',rangegrph else: rangegrph[0] = np.float(vrng[0]) rangegrph[1] = np.float(vrng[1]) return rangegrph def format_axes(axesS,char): """ Function to provide the format for the ticks of the axes in a figure axesS= [char] separated list of ('auto' for all 'auto' values) xs: style of x-axis ('auto' for 'pretty') 'Nfix', values computed at even 'Ndx' 'Vfix', values computed at even 'Ndx' increments 'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10) xf: format of the ticks labels at the x-axis ('auto' for '5g') Nx: number of ticks labels at the x-axis ('auto' for 5) orx: angle of orientation of ticks at the x-axis ('auto' for horizontal) ys: style of y-axis ('auto' for 'pretty') yf: format of the ticks labels at the y-axis ('auto' for '5g') Ny: number of ticks labels at the y-axis ('auto' for 5) ory: angle of orientation of ticks at the y-axis ('auto' for horizontal) >>> format_axes('auto',',') ('pretty', '5g', 5, 0.0, 'pretty', '5g', 5, 0.0) >>> format_axes('pretty,auto,7,auto,Vfix,.0f,30,45.',',') ('pretty', '5g', 7, 0.0, 'Vfix', '.0f', 30, 45.0) """ fname = 'format_axes' if axesS == 'auto': xs = 'pretty' xf = '5g' Nx = 5 orx = 0. ys = 'pretty' yf = '5g' Ny = 5 ory = 0. else: Sv = axesS.split(char) if Sv[0] == 'auto': xs = 'pretty' else: xs = Sv[0] if Sv[1] == 'auto': xf = '5g' else: xf = Sv[1] if Sv[2] == 'auto': Nx = 5 else: if Sv[0] == 'pretty' or Sv[0] == 'Nfix': Nx = int(Sv[2]) else: Nx = np.float(Sv[2]) if Sv[3] == 'auto': orx = 0. else: orx = np.float(Sv[3]) if Sv[4] == 'auto': ys = 'pretty' else: ys = Sv[4] if Sv[5] == 'auto': yf = '5g' else: yf = Sv[5] if Sv[6] == 'auto': Ny = 5 else: if Sv[4] == 'pretty' or Sv[4] == 'Nfix': Ny = int(Sv[6]) else: Ny = np.float(Sv[6]) if Sv[7] == 'auto': ory = 0. else: ory = np.float(Sv[7]) return xs, xf, Nx, orx, ys, yf, Ny, ory ####### ###### ##### #### ### ## # END of GENERIC functions def plot_TimeSeries(valtimes, vunits, tunits, hfileout, vtit, ttit, tkind, tformat, \ tit, linesn, lloc, kfig,coll,ptl,ptf): """ Function to draw time-series valtimes= list of arrays to plot [vals1[1values, 1times], [...,valsM[Mvals,Mtimes]]) vunits= units of the values tunits= units of the times hfileout= header of the output figure. Final name: [hfileout]_[vtit].[kfig] vtit= variable title to be used in the graph ttit= time title to be used in the graph tkind= kind of time values to appear in the 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 tformat= desired format of times tit= title of the graph linesn= list of values fot the legend lloc= location of the legend (0, 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' kfig= type of figure: jpg, png, pds, ps coll= ',' list of colors for the lines, None for automatic, single value all the same ptl= ',' list of type of points for the lines, None for automatic, single value all the same ptf= frequency of point plotting, 'all' for all time steps """ fname = 'plot_TimeSeries' if valtimes == 'h': print fname + '_____________________________________________________________' print plot_TimeSeries.__doc__ quit() Nlines = len(valtimes) # Canging line kinds every 7 lines (end of standard colors) Npts = len(pointkindsauto) linekinds = [] pointkinds = [] if ptl is None: if ptf is None: for ptype in range(Npts): for ip in range(7): linekinds.append(pointkindsauto[ptype] + '-') else: for ptype in range(Npts): for ip in range(7): linekinds.append('-') pointkinds.append(pointkindsauto[ptype]) else: if ptf is None: if len(ptl) > 1: for pt in ptl: linekinds.append(pt + '-') else: for il in range(Nlines): linekinds.append(ptl+'-') else: if len(ptl) > 1: for pt in ptl: linekinds.append('-') pointkinds.append(pt) else: for il in range(Nlines): linekinds.append('-') pointkinds = ptl Nvalues = [] Ntimes = [] for il in range(Nlines): array = valtimes[il] if Nlines == 1: print warnmsg print ' ' + fname + ': drawing only one line!' Nvalues.append(array.shape[1]) Ntimes.append(array.shape[0]) tmin = np.min(array[1]) tmax = np.max(array[1]) vmin = np.min(array[0]) vmax = np.max(array[0]) else: Nvalues.append(array.shape[1]) Ntimes.append(array.shape[0]) tmin = np.min(array[1,:]) tmax = np.max(array[1,:]) vmin = np.min(array[0,:]) vmax = np.max(array[0,:]) if il == 0: xmin = tmin xmax = tmax ymin = vmin ymax = vmax else: if tmin < xmin: xmin = tmin if tmax > xmax: xmax = tmax if vmin < ymin: ymin = vmin if vmax > ymax: ymax = vmax dx = np.max(Ntimes) dy = np.min(Nvalues) plt.rc('text', usetex=True) print vtit if vtit == 'ps': plt.ylim(98000.,ymax) else: plt.ylim(ymin,ymax) plt.xlim(xmin,xmax) # print 'x lim:',xmin,xmax # print 'y lim:',ymin,ymax N7lines=0 for il in range(Nlines): array = valtimes[il] if vtit == 'ps': array[0,:] = np.where(array[0,:] < 98000., None, array[0,:]) if ptf is None: if coll is None: if ptf is None: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il]) else: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il]) plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il], \ label=linesn[il]) else: if ptf is None: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \ color=coll[il]) else: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \ color=coll[il]) plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il], \ label=linesn[il], color=coll[il]) else: if coll is None: if ptf is None: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il]) else: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il]) plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il], \ label=linesn[il]) else: if ptf is None: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \ color=coll[il]) else: plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \ color=coll[il]) plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il], \ label=linesn[il], color=coll[il]) timevals = np.arange(xmin,xmax)*1. tpos, tlabels = CFtimes_plot(timevals, tunits, tkind, tformat) if len(tpos) > 10: print warnmsg print ' ' + fname + ': with "' + tkind + '" there are', len(tpos), 'xticks !' plt.xticks(tpos, tlabels) # plt.Axes.set_xticklabels(tlabels) plt.legend(loc=lloc) plt.xlabel(ttit) plt.ylabel(vtit + " (" + vunits + ")") plt.title(tit.replace('_','\_').replace('&','\&')) figname = hfileout + '_' + vtit output_kind(kfig, figname, True) return #Nt = 50 #Nlines = 3 #vtvalsv = [] #valsv = np.zeros((2,Nt), dtype=np.float) ## First #valsv[0,:] = np.arange(Nt) #valsv[1,:] = np.arange(Nt)*180. #vtvalsv.append(valsv) #del(valsv) #valsv = np.zeros((2,Nt/2), dtype=np.float) ## Second #valsv[0,:] = np.arange(Nt/2) #valsv[1,:] = np.arange(Nt/2)*180.*2. #vtvalsv.append(valsv) #del(valsv) #valsv = np.zeros((2,Nt/4), dtype=np.float) ## Third #valsv[0,:] = np.arange(Nt/4) #valsv[1,:] = np.arange(Nt/4)*180.*4. #vtvalsv.append(valsv) #del(valsv) #varu='mm' #timeu='seconds' #title='test' #linesname = ['line 1', 'line 2', 'line 3'] #plot_TimeSeries(vtvalsv, units_lunits(varu), timeu, 'test', 'vartest', 'time', title, linesname, 'png') #quit() def plot_points(xval, yval, vlon, vlat, extravals, extrapar, vtit, mapv, figk, color,\ labels, lloc, kfig, figname): """ plotting points [x/yval]: x,y values to plot vlon= 2D-matrix with the longitudes vlat= 2D-matrix with the latitudes extravals= extra values to be added into the plot (None for nothing) extrapar= [varname, min, max, cbar, varunits] of the extra variable vtit= title of the graph ('|' for spaces) mapv= 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 figK= kind of figure 'legend': only points in the map with the legend with the names 'labelled',[txtsize],[txtcol]: points with the names and size, color of text color= color for the points/labels ('auto', for "red") labels= list of labels for the points (None, no labels) lloc = localisation of the legend kfig= kind of figure (jpg, pdf, png) figname= name of the figure """ fname = 'plot_points' # Canging line kinds every 7 pts (end of standard colors) ptkinds=['.','x','o','*','+','8','>','D','h','p','s'] Npts = len(xval) if Npts > len(ptkinds)*7: print errormsg print ' ' + fname + ': too many',Npts,'points!!' print " enlarge 'ptkinds' list" quit(-1) N7pts = 0 if color == 'auto': ptcol = "red" else: ptcol = color dx=vlon.shape[1] dy=vlat.shape[0] plt.rc('text', usetex=True) if not mapv is None: # vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:]) # xvala = np.array(xval) # xvala = np.where(xvala < 0., 360. + xvala, xvala) # xval = list(xvala) map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = np.min(vlon) xlon = np.max(vlon) nlat = np.min(vlat) xlat = np.max(vlat) lon2 = vlon[dy/2,dx/2] lat2 = vlat[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": map projecion '" + map_proj + "' not ready!!" print ' available: cyl, lcc' quit(-1) # lons, lats = np.meshgrid(vlon, vlat) # lons = np.where(lons < 0., lons + 360., lons) x,y = m(vlon,vlat) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) else: x = vlon y = vlat # plt.xlim(0,dx-1) # plt.ylim(0,dy-1) # Extra values if extravals is not None: plt.pcolormesh(x, y, extravals, cmap=plt.get_cmap(extrapar[3]), \ vmin=extrapar[1], vmax=extrapar[2]) cbar = plt.colorbar() cbar.set_label(extrapar[0].replace('_','\_') +'('+ units_lunits(extrapar[4])+\ ')') if labels is not None: for iv in range(len(xval)): if np.mod(iv,7) == 0: N7pts = N7pts + 1 # print iv,xval[iv],yval[iv],labels[iv],ptkinds[N7pts] plt.plot(xval[iv],yval[iv], ptkinds[N7pts],label=labels[iv]) if figk[0:8] == 'labelled': txtsize=int(figk.split(',')[1]) txtcol=figk.split(',')[2] for iv in range(len(xval)): plt.annotate(labels[iv], xy=(xval[iv],yval[iv]), xycoords='data', \ fontsize=txtsize, color=txtcol) elif figk == 'legend': plt.legend(loc=lloc) else: plt.plot(xval, yval, '.', color=ptcol) graphtit = vtit.replace('_','\_').replace('&','\&') plt.title(graphtit.replace('|', ' ')) output_kind(kfig, figname, True) return def plot_list_points(xval, yval, varxn, varyn, vtit, figk, color, ptk, pts, labels, lloc, kfig, figname): """ plotting points [x/yval]: x,y values to plot var[x/y]n: names of the x,y variables vtit= title of the graph ('|' for spaces) figK= kind of figure 'legend': only points in the map with the legend with the names 'labelled',[txtsize],[txtcol]: points with the names and size, color of text color= color for the points/labels ('auto', for "red") ptk= kind of point pts= point size labels= list of labels for the points (None, no labels) lloc = localisation of the legend kfig= kind of figure (jpg, pdf, png) figname= name of the figure """ fname = 'plot_points' # Canging line kinds every 7 pts (end of standard colors) ptkinds=['.','x','o','*','+','8','>','D','h','p','s'] if color == 'auto': ptcol = "red" else: ptcol = color plt.rc('text', usetex=True) if labels is not None: for iv in range(len(xval)): # print iv,xval[iv],yval[iv],labels[iv],ptkinds[N7pts] plt.plot(xval[iv],yval[iv], ptk, markersize=pts, label=labels[iv]) if figk[0:8] == 'labelled': txtsize=int(figk.split(',')[1]) txtcol=figk.split(',')[2] for iv in range(len(xval)): plt.annotate(labels[iv], xy=(xval[iv],yval[iv]), xycoords='data', \ fontsize=txtsize, color=txtcol) elif figk == 'legend': plt.legend(loc=lloc) else: plt.plot(xval, yval, ptk, markersize=pts, color=ptcol) plt.xlabel(varxn) plt.ylabel(varyn) graphtit = vtit.replace('_','\_').replace('&','\&') plt.title(graphtit.replace('|', ' ')) output_kind(kfig, figname, True) return def plot_2Dfield(varv,dimn,colorbar,vn,vx,unit,olon,olat,ifile,vtit,zvalue,time,tk, \ tkt,tobj,tvals,tind,kfig,mapv,reva): """ Adding labels and other staff to the graph varv= 2D values to plot dimn= dimension names to plot colorbar= name of the color bar to use vn,vm= minmum and maximum values to plot unit= units of the variable olon,olat= longitude, latitude objects ifile= name of the input file vtit= title of the variable zvalue= value on the z axis time= value on the time axis tk= kind of time (WRF) tkt= kind of time taken tobj= tim object tvals= values of the time variable tind= time index kfig= kind of figure (jpg, pdf, png) mapv= map characteristics: [proj],[res] see full documentation: http://matplotlib.org/basemap/ [proj]: projection * 'cyl', cilindric [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full reva= reverse the axes (x-->y, y-->x) """ ## import matplotlib as mpl ## mpl.use('Agg') ## import matplotlib.pyplot as plt if reva: print ' reversing the axes of the figure (x-->y, y-->x)!!' varv = np.transpose(varv) dimn0 = [] dimn0.append(dimn[1] + '') dimn0.append(dimn[0] + '') dimn = dimn0 fname = 'plot_2Dfield' dx=varv.shape[1] dy=varv.shape[0] plt.rc('text', usetex=True) # plt.rc('font', family='serif') if not mapv is None: if len(olon[:].shape) == 3: lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,]) lat0 = olat[0,] elif len(olon[:].shape) == 2: lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:]) lat0 = olat[:] elif len(olon[:].shape) == 1: lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:]) lat00 = olat[:] lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) for iy in range(len(lat00)): lon0[iy,:] = lon00 for ix in range(len(lon00)): lat0[:,ix] = lat00 map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = lon0[0,0] xlon = lon0[dy-1,dx-1] nlat = lat0[0,0] xlat = lat0[dy-1,dx-1] lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] print ' lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) if len(olon[:].shape) == 1: lons, lats = np.meshgrid(olon[:], olat[:]) else: lons = olon[0,:] lats = olat[:,0] lons = np.where(lons < 0., lons + 360., lons) x,y = m(lons,lats) plt.pcolormesh(x,y,varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx) cbar = plt.colorbar() m.drawcoastlines() # if (nlon > 180. or xlon > 180.): # nlon0 = nlon # xlon0 = xlon # if (nlon > 180.): nlon0 = nlon - 360. # if (xlon > 180.): xlon0 = xlon - 360. # meridians = gen.pretty_int(nlon0,xlon0,5) # meridians = np.where(meridians < 0., meridians + 360., meridians) # else: # meridians = gen.pretty_int(nlon,xlon,5) meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) else: plt.xlim(0,dx-1) plt.ylim(0,dy-1) plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx) cbar = plt.colorbar() plt.xlabel(dimn[1].replace('_','\_')) plt.ylabel(dimn[0].replace('_','\_')) # set the limits of the plot to the limits of the data # plt.axis([x.min(), x.max(), y.min(), y.max()]) # plt.plot(varv) cbar.set_label(unit) figname = ifile.replace('.','_') + '_' + vtit graphtit = vtit.replace('_','\_').replace('&','\&') if zvalue != 'null': graphtit = graphtit + ' at z= ' + zvalue figname = figname + '_z' + zvalue if tkt == 'tstep': graphtit = graphtit + ' at time-step= ' + time.split(',')[1] figname = figname + '_t' + time.split(',')[1].zfill(4) elif tkt == 'CFdate': graphtit = graphtit + ' at ' + tobj.strfmt("%Y%m%d%H%M%S") figname = figname + '_t' + tobj.strfmt("%Y%m%d%H%M%S") if tk == 'WRF': # datev = str(timevals[timeind][0:9]) datev = tvals[tind][0] + tvals[tind][1] + tvals[tind][2] + \ timevals[timeind][3] + timevals[timeind][4] + timevals[timeind][5] + \ timevals[timeind][6] + timevals[timeind][7] + timevals[timeind][8] + \ timevals[timeind][9] # timev = str(timevals[timeind][11:18]) timev = timevals[timeind][11] + timevals[timeind][12] + \ timevals[timeind][13] + timevals[timeind][14] + timevals[timeind][15] + \ timevals[timeind][16] + timevals[timeind][17] + timevals[timeind][18] graphtit = vtit.replace('_','\_') + ' (' + datev + ' ' + timev + ')' plt.title(graphtit) output_kind(kfig, figname, True) return def plot_2Dfield_easy(varv,dimxv,dimyv,dimn,colorbar,vn,vx,unit,ifile,vtit,kfig,reva): """ Adding labels and other staff to the graph varv= 2D values to plot dim[x/y]v = values at the axes of x and y dimn= dimension names to plot colorbar= name of the color bar to use vn,vm= minmum and maximum values to plot unit= units of the variable ifile= name of the input file vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= reverse the axes (x-->y, y-->x) """ ## import matplotlib as mpl ## mpl.use('Agg') ## import matplotlib.pyplot as plt fname = 'plot_2Dfield' if reva: print ' reversing the axes of the figure (x-->y, y-->x)!!' varv = np.transpose(varv) dimn0 = [] dimn0.append(dimn[1] + '') dimn0.append(dimn[0] + '') dimn = dimn0 if len(dimyv.shape) == 2: x = np.transpose(dimyv) else: if len(dimxv.shape) == 2: ddx = len(dimyv) ddy = dimxv.shape[1] else: ddx = len(dimyv) ddy = len(dimxv) x = np.zeros((ddy,ddx), dtype=np.float) for j in range(ddy): x[j,:] = dimyv if len(dimxv.shape) == 2: y = np.transpose(dimxv) else: if len(dimyv.shape) == 2: ddx = dimyv.shape[0] ddy = len(dimxv) else: ddx = len(dimyv) ddy = len(dimxv) y = np.zeros((ddy,ddx), dtype=np.float) for i in range(ddx): y[:,i] = dimxv else: if len(dimxv.shape) == 2: x = dimxv else: x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float) for j in range(len(dimyv)): x[j,:] = dimxv if len(dimyv.shape) == 2: y = dimyv else: y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float) for i in range(len(dimxv)): x[:,i] = dimyv dx=varv.shape[1] dy=varv.shape[0] plt.rc('text', usetex=True) plt.xlim(0,dx-1) plt.ylim(0,dy-1) plt.pcolormesh(x, y, varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx) # plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx) cbar = plt.colorbar() plt.xlabel(dimn[1].replace('_','\_')) plt.ylabel(dimn[0].replace('_','\_')) # set the limits of the plot to the limits of the data plt.axis([x.min(), x.max(), y.min(), y.max()]) # if varv.shape[1] / varv.shape[0] > 10: # plt.axes().set_aspect(0.001) # else: # plt.axes().set_aspect(np.float(varv.shape[0])/np.float(varv.shape[1])) cbar.set_label(unit) figname = ifile.replace('.','_') + '_' + vtit graphtit = vtit.replace('_','\_').replace('&','\&') plt.title(graphtit) output_kind(kfig, figname, True) return def plot_Trajectories(lonval, latval, linesn, olon, olat, lonlatLims, gtit, kfig, \ mapv, obsname): """ plotting Trajectories [lon/latval]= lon,lat values to plot (as list of vectors) linesn: name of the lines o[lon/lat]= object with the longitudes and the latitudes of the map to plot lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] gtit= title of the graph kfig= kind of figure (jpg, pdf, png) mapv= 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 obsname= name of the observations in graph (can be None for without). Observational trajectory would be the last one """ fname = 'plot_Trajectories' if lonval == 'h': print fname + '_____________________________________________________________' print plot_Trajectories.__doc__ quit() # Canging line kinds every 7 lines (end of standard colors) linekinds=['.-','x-','o-'] Ntraj = len(lonval) if obsname is not None: Ntraj = Ntraj - 1 N7lines = 0 plt.rc('text', usetex=True) if not mapv is None: if len(olon[:].shape) == 3: # lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,]) lon0 = olon[0,] lat0 = olat[0,] elif len(olon[:].shape) == 2: # lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:]) lon0 = olon[:] lat0 = olat[:] elif len(olon[:].shape) == 1: # lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:]) lon00 = olon[:] lat00 = olat[:] lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) for iy in range(len(lat00)): lon0[iy,:] = lon00 for ix in range(len(lon00)): lat0[:,ix] = lat00 map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lon0.shape[0] nlon = lon0[0,0] xlon = lon0[dy-1,dx-1] nlat = lat0[0,0] xlat = lat0[dy-1,dx-1] lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] if lonlatLims is not None: plt.xlim(lonlatLims[0], lonlatLims[2]) plt.ylim(lonlatLims[1], lonlatLims[3]) if map_proj == 'cyl': nlon = lonlatLims[0] nlat = lonlatLims[1] xlon = lonlatLims[2] xlat = lonlatLims[3] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) if len(olon.shape) == 3: # lons, lats = np.meshgrid(olon[0,:,:], olat[0,:,:]) lons = olon[0,:,:] lats = olat[0,:,:] elif len(olon.shape) == 2: # lons, lats = np.meshgrid(olon[:,:], olat[:,:]) lons = olon[:,:] lats = olat[:,:] else: dx = olon.shape dy = olat.shape # print errormsg # print ' ' + fname + ': shapes of lon/lat objects', olon.shape, \ # 'not ready!!!' for il in range(Ntraj): plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il]) if il == 6: N7lines = N7lines + 1 m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: if len(olon.shape) == 3: dx = olon.shape[2] dy = olon.shape[1] elif len(olon.shape) == 2: dx = olon.shape[1] dy = olon.shape[0] else: dx = olon.shape dy = olat.shape # print errormsg # print ' ' + fname + ': shapes of lon/lat objects', olon.shape, \ # 'not ready!!!' if lonlatLims is not None: plt.xlim(lonlatLims[0], lonlatLims[2]) plt.ylim(lonlatLims[1], lonlatLims[3]) else: plt.xlim(np.min(olon[:]),np.max(olon[:])) plt.ylim(np.min(olat[:]),np.max(olat[:])) for il in range(Ntraj): plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il]) if il == 6: N7lines = N7lines + 1 plt.xlabel('x-axis') plt.ylabel('y-axis') figname = 'trajectories' graphtit = gtit if obsname is not None: plt.plot(lonval[Ntraj], latval[Ntraj], linestyle='-', color='k', \ linewidth=3, label= obsname) plt.title(graphtit.replace('_','\_').replace('&','\&')) plt.legend() output_kind(kfig, figname, True) return def plot_topo_geogrid(varv, olon, olat, mint, maxt, lonlatLims, gtit, kfig, mapv, \ closeif): """ plotting geo_em.d[nn].nc topography from WPS files plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv) varv= topography values o[lon/lat]= longitude and latitude objects [min/max]t: minimum and maximum values of topography to draw lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] gtit= title of the graph kfig= kind of figure (jpg, pdf, png) 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 closeif= Boolean value if the figure has to be closed """ fname = 'plot_topo_geogrid' if varv == 'h': print fname + '_____________________________________________________________' print plot_topo_geogrid.__doc__ quit() dx=varv.shape[1] dy=varv.shape[0] plt.rc('text', usetex=True) # plt.rc('font', family='serif') if not mapv is None: if len(olon[:].shape) == 3: lon0 = olon[0,] lat0 = olat[0,] elif len(olon[:].shape) == 2: lon0 = olon[:] lat0 = olat[:] elif len(olon[:].shape) == 1: lon00 = olon[:] lat00 = olat[:] lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) for iy in range(len(lat00)): lon0[iy,:] = lon00 for ix in range(len(lon00)): lat0[:,ix] = lat00 map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lon0.shape[0] if lonlatLims is not None: print ' ' + fname + ': cutting the domain to plot !!!!' print ' limits: W-E', lonlatLims[0], lonlatLims[2] print ' limits: N-S', lonlatLims[1], lonlatLims[3] nlon = lonlatLims[0] xlon = lonlatLims[2] nlat = lonlatLims[1] xlat = lonlatLims[3] if map_proj == 'lcc': lon2 = (lonlatLims[0] + lonlatLims[2])/2. lat2 = (lonlatLims[1] + lonlatLims[3])/2. else: nlon = np.min(lon0) xlon = np.max(lon0) nlat = np.min(lat0) xlat = np.max(lat0) lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] plt.xlim(nlon, xlon) plt.ylim(nlat, xlat) print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": map projection '" + map_proj + "' not ready !!" quit(-1) if len(olon[:].shape) == 1: lons, lats = np.meshgrid(olon[:], olat[:]) else: if len(olon[:].shape) == 3: lons = olon[0,:,:] lats = olat[0,:,:] else: lons = olon[:] lats = olat[:] x,y = m(lons,lats) plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt) cbar = plt.colorbar() m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: print emsg print ' ' + fname + ': A projection parameter is needed None given !!' quit(-1) figname = 'domain' graphtit = gen.latex_text(gtit) cbar.set_label('height ($m$)') plt.title(graphtit) output_kind(kfig, figname, closeif) return def plot_topo_geogrid_boxes(varv, boxesX, boxesY, boxlabels, olon, olat, mint, maxt, \ lonlatLims, gtit, kfig, mapv, gloc, gsiz, closeif): """ plotting geo_em.d[nn].nc topography from WPS files plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv) varv= topography values boxesX/Y= 4-line sets to draw the boxes boxlabels= labels for the legend of the boxes o[lon/lat]= longitude and latitude objects [min/max]t: minimum and maximum values of topography to draw lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] gtit= title of the graph kfig= kind of figure (jpg, pdf, png) 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 gloc= location of the legend (0, 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' gsiz= fontsize of the legend closeif= Boolean value if the figure has to be closed """ fname = 'plot_topo_geogrid' if varv == 'h': print fname + '_____________________________________________________________' print plot_topo_geogrid.__doc__ quit() cols = color_lines(len(boxlabels)) dx=varv.shape[1] dy=varv.shape[0] plt.rc('text', usetex=True) # plt.rc('font', family='serif') if not mapv is None: if len(olon[:].shape) == 3: lon0 = olon[0,] lat0 = olat[0,] elif len(olon[:].shape) == 2: lon0 = olon[:] lat0 = olat[:] elif len(olon[:].shape) == 1: lon00 = olon[:] lat00 = olat[:] lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) for iy in range(len(lat00)): lon0[iy,:] = lon00 for ix in range(len(lon00)): lat0[:,ix] = lat00 map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lon0.shape[0] if lonlatLims is not None: print ' ' + fname + ': cutting the domain to plot !!!!' print ' limits: W-E', lonlatLims[0], lonlatLims[2] print ' limits: N-S', lonlatLims[1], lonlatLims[3] nlon = lonlatLims[0] xlon = lonlatLims[2] nlat = lonlatLims[1] xlat = lonlatLims[3] if map_proj == 'lcc': lon2 = (lonlatLims[0] + lonlatLims[2])/2. lat2 = (lonlatLims[1] + lonlatLims[3])/2. else: nlon = np.min(lon0) xlon = np.max(lon0) nlat = np.min(lat0) xlat = np.max(lat0) lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] plt.xlim(nlon, xlon) plt.ylim(nlat, xlat) print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' does not exist!!" print ' existing ones: cyl, lcc' quit(-1) if len(olon[:].shape) == 1: lons, lats = np.meshgrid(olon[:], olat[:]) else: if len(olon[:].shape) == 3: lons = olon[0,:,:] lats = olat[0,:,:] else: lons = olon[:] lats = olat[:] x,y = m(lons,lats) plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt) cbar = plt.colorbar() Nboxes = len(boxesX)/4 for ibox in range(Nboxes): plt.plot(boxesX[ibox*4], boxesY[ibox*4], linestyle='-', linewidth=3, \ label=gen.latex_text(boxlabels[ibox].replace('!',' ')), color=cols[ibox]) plt.plot(boxesX[ibox*4+1], boxesY[ibox*4+1], linestyle='-', linewidth=3, \ color=cols[ibox]) plt.plot(boxesX[ibox*4+2], boxesY[ibox*4+2], linestyle='-', linewidth=3, \ color=cols[ibox]) plt.plot(boxesX[ibox*4+3], boxesY[ibox*4+3], linestyle='-', linewidth=3, \ color=cols[ibox]) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: print emsg print ' ' + fname + ': A projection parameter is needed None given !!' quit(-1) figname = 'domain_boxes' graphtit = gen.latex_text(gtit) cbar.set_label('height ($m$)') plt.title(graphtit) plt.legend(loc=gloc, prop={'size':gsiz}) output_kind(kfig, figname, closeif) return def plot_2D_shadow(varsv, vnames, dimxv, dimyv, dimxu, dimyu, xaxv, yaxv, dimn, \ cbarv, vs, uts, vtit, kfig, reva, mapv, ifclose): """ Adding labels and other staff to the graph varsv= 2D values to plot with shading vnames= variable names for the figure dim[x/y]v = values at the axes of x and y dim[x/y]u = units at the axes of x and y xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] dimn= dimension names to plot cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar: name of the color bar to use cbarfmt: format of the numbers in the colorbar cbaror: orientation of the colorbar vs= minmum and maximum values to plot in shadow 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) uts= units of the variable to shadow vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= ('|' for combination) * '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', lambert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full ifclose= boolean value whether figure should be close (finish) or not """ ## import matplotlib as mpl ## mpl.use('Agg') ## import matplotlib.pyplot as plt fname = 'plot_2D_shadow' # print dimyv[73,21] # dimyv[73,21] = -dimyv[73,21] # print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv) if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) # Axis ticks # Usually axis > x must be the lon, thus... dimxv0 = dimxv.copy() dimyv0 = dimyv.copy() dxn = dimxv.min() dxx = dimxv.max() dyn = dimyv.min() dyx = dimyv.max() if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2])) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) if yaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2])) elif yaxv[0] == 'Nfix': dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2])) elif yaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,yaxv[2]) dimxl0 = [] for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1])) dimyl0 = [] for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')' dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')' if mapv is not None: pixkind = 'data' else: # No following data values dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0)) dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0)) pixkind = 'fixpixel' if reva is not None: varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 if len(dimxv[:].shape) == 3: xdims = '1,2' elif len(dimxv[:].shape) == 2: xdims = '0,1' elif len(dimxv[:].shape) == 1: xdims = '0' if len(dimyv[:].shape) == 3: ydims = '1,2' elif len(dimyv[:].shape) == 2: ydims = '0,1' elif len(dimyv[:].shape) == 1: ydims = '0' # lon0 = dimxv # lat0 = dimyv lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims) if not mapv is None: map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lat0.shape[0] # nlon = lon0[0,0] # xlon = lon0[dy-1,dx-1] # nlat = lat0[0,0] # xlat = lat0[dy-1,dx-1] nlon = np.min(lon0) xlon = np.max(lon0) nlat = np.min(lat0) xlat = np.max(lat0) # Thats too much! :) # if lonlatLims is not None: # print ' ' + fname + ': cutting the domain to plot !!!!' # plt.xlim(lonlatLims[0], lonlatLims[2]) # plt.ylim(lonlatLims[1], lonlatLims[3]) # print ' limits: W-E', lonlatLims[0], lonlatLims[2] # print ' limits: N-S', lonlatLims[1], lonlatLims[3] # if map_proj == 'cyl': # nlon = lonlatLims[0] # nlat = lonlatLims[1] # xlon = lonlatLims[2] # xlat = lonlatLims[3] # elif map_proj == 'lcc': # lon2 = (lonlatLims[0] + lonlatLims[2])/2. # lat2 = (lonlatLims[1] + lonlatLims[3])/2. # nlon = lonlatLims[0] # xlon = lonlatLims[2] # nlat = lonlatLims[1] # xlat = lonlatLims[3] lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": map projection '" + map_proj + "' not defined!!!" print ' available: cyl, lcc' quit(-1) x,y = m(lon0,lat0) else: # No following data values #x = dimxv #y = dimyv x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv)) y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv)) # Changing limits of the colors vsend = graphic_range(vs,varsv) plt.rc('text', usetex=True) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], vmax=vsend[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) if not mapv is None: if cbarv[0] == 'gist_gray': m.drawcoastlines(color="red") else: m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,xaxv[2]) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,yaxv[2]) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: plt.xlabel(dimxT) plt.ylabel(dimyT) # Using output from transform # print 'x min:', x.min(),'max:',x.max() # print 'y min:', y.min(),'max:',y.max() # txpos = gen.pretty_int(x.min(),x.max(),5) # typos = gen.pretty_int(y.min(),y.max(),5) # txlabels = list(txpos) # for i in range(len(txlabels)): txlabels[i] = '{:6g}'.format(txlabels[i]) # tylabels = list(typos) # for i in range(len(tylabels)): tylabels[i] = '{:6g}'.format(tylabels[i]) plt.axis([x.min(), x.max(), y.min(), y.max()]) if mapv is None: # plt.xticks(txpos, txlabels) # plt.yticks(typos, tylabels) plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) # units labels cbar.set_label(gen.latex_text(vnames) + ' (' + units_lunits(uts) + ')') figname = '2Dfields_shadow' graphtit = gen.latex_text(vtit) plt.title(gen.latex_text(graphtit)) output_kind(kfig, figname, ifclose) return #Nvals=50 #vals1 = np.zeros((Nvals,Nvals), dtype= np.float) #for j in range(Nvals): # for i in range(Nvals): # vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.) #plot_2D_shadow(vals1, 'var1', np.arange(50)*1., np.arange(50)*1., 'ms-1', \ # 'm', ['lat','lon'], 'rainbow', [0, Nvals], 'ms-1', 'test var1', 'pdf', 'None', \ # None, True) #quit() def plot_2D_shadow_time(varsv, vnames, dimxv, dimyv, dimxu, dimyu, dimn, vaxv, cbarv,\ vs, uts, vtit, kfig, reva, taxis, tpos, tlabs, ifclose): """ Plotting a 2D field with one of the axes being time varsv= 2D values to plot with shading vnames= shading variable name for the figure dim[x/y]v= values at the axes of x and y dim[x/y]u= units at the axes of x and y dimn= dimension names to plot vaxv= list with the non-temporal axis paramteres [style, format, number and orientation] cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar: name of the color bar to use cbarfmt: format of the numbers in the colorbar cbaror: orientation of the colorbar vs= minmum and maximum values to plot in shadow 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) uts= units of the variable to shadow vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y taxis= Which is the time-axis tpos= positions of the time ticks tlabs= labels of the time ticks ifclose= boolean value whether figure should be close (finish) or not """ fname = 'plot_2D_shadow_time' if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow_time.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) dimxv0 = dimxv.copy() dimyv0 = dimyv.copy() dxn = dimxv.min() dxx = dimxv.max() dyn = dimyv.min() dyx = dimyv.max() # Definning ticks labels if taxis == 'x': dimxt0 = tpos dimxl0 = tlabs dimxT0 = dimxu if vaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,vaxv[2])) elif vaxv[0] == 'Nfix': dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*vaxv[2])) elif vaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,vaxv[2]) dimyl0 = [] for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], \ style=vaxv[1])) dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')' # No following data values dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0)) pixkind = 'fixpixel' else: if vaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,vaxv[2])) elif vaxv[0] == 'Nfix': dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*vaxv[2])) elif vaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,vaxv[2]) dimxl0 = [] for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], \ style=vaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')' dimyt0 = tpos dimyl0 = tlabs dimyT0 = dimyu # No following data values dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0)) pixkind = 'fixpixel' if reva is not None: varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 if len(dimxv.shape) == 3: dxget='1,2' elif len(dimxv.shape) == 2: dxget='0,1' elif len(dimxv.shape) == 1: dxget='0' else: print errormsg print ' ' + fname + ': shape of x-values:',dimxv.shape,'not ready!!' quit(-1) if len(dimyv.shape) == 3: dyget='1,2' elif len(dimyv.shape) == 2: dyget='0,1' elif len(dimyv.shape) == 1: dyget='0' else: print errormsg print ' ' + fname + ': shape of y-values:',dimyv.shape,'not ready!!' quit(-1) # No following data values x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv)) y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv)) # Changing limits of the colors vsend = graphic_range(vs,varsv) plt.rc('text', usetex=True) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], vmax=vsend[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) if taxis == 'x': plt.xticks((dimxt-dxn)/(dxx-dxn), dimxl) plt.yticks(dimyt, dimyl, rotation=vaxv[3]) else: plt.xticks(dimxt, dimxl, rotation=vaxv[3]) plt.yticks((dimyt-dyn)/(dyx-dyn), dimyl) plt.xlabel(dimxT) plt.ylabel(dimyT) plt.axis([x.min(), x.max(), y.min(), y.max()]) # units labels cbar.set_label(gen.latex_text(vnames) + ' (' + units_lunits(uts) + ')') figname = '2Dfields_shadow_time' graphtit = gen.latex_text(vtit) plt.title(graphtit) output_kind(kfig, figname, ifclose) return def plot_2D_shadow_contour(varsv, varcv, vnames, dimxv, dimyv, dimxu, dimyu, xaxv, \ yaxv, dimn, cbarv, ckind, clabfmt, vs, vc, uts, vtit, kfig, reva, mapv, ifclose): """ Adding labels and other staff to the graph varsv= 2D values to plot with shading varcv= 2D values to plot with contours vnames= variable names for the figure dim[x/y]v = values at the axes of x and y dim[x/y]u = units at the axes of x and y xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] dimn= dimension names to plot cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar= name of the color bar to use cbarfmt= format of the numbers in the colorbar cbaror= orientation of the colorbar ckind= contour kind '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 plot (None, no labels) vs= minmum and maximum values to plot in shadow vc= vector with the levels for the contour uts= units of the variable [u-shadow, u-contour] vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= * '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 ifclose= boolean value whether figure should be close (finish) or not """ fname = 'plot_2D_shadow_contour' if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow_contour.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) if len(varcv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable contour rank:', varcv.shape, \ 'is has to be 2D!!' quit(-1) # Axis ticks # Usually axis > x must be the lon, thus... dimxv0 = dimxv.copy() dimyv0 = dimyv.copy() dxn = dimxv.min() dxx = dimxv.max() dyn = dimyv.min() dyx = dimyv.max() if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(0.,1.,1./xaxv[2]) dimxv0 = np.arange(0.,1.,1./varsv.shape[1]) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) if yaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2])) elif yaxv[0] == 'Nfix': dimyt0 = np.arange(0.,1.,1./yaxv[2]) dimyv0 = np.arange(0.,1.,1./varsv.shape[0]) elif yaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,yaxv[2]) dimxl0 = [] if xaxv[0] != 'Nfix': for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1])) else: for i in range(len(dimxt0)): xv=dxn+i*(dxx-dxn)/(1.*xaxv[2]) dimxl0.append('{:{style}}'.format(xv, style=xaxv[1])) dimyl0 = [] if yaxv[0] != 'Nfix': for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1])) else: for i in range(len(dimyt0)): yv=dyn+i*(dyx-dyn)/(1.*yaxv[2]) dimyl0.append('{:{style}}'.format(yv, style=yaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')' dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')' if reva is not None: varcv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varcv, reva) varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 if not mapv is None: if len(dimxv[:].shape) == 3: lon0 = dimxv[0,] lat0 = dimyv[0,] elif len(dimxv[:].shape) == 2: lon0 = dimxv[:] lat0 = dimyv[:] elif len(dimxv[:].shape) == 1: lon00 = dimxv[:] lat00 = dimyv[:] lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float ) for iy in range(len(lat00)): lon0[iy,:] = lon00 for ix in range(len(lon00)): lat0[:,ix] = lat00 map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lon0.shape[0] nlon = np.min(lon0) xlon = np.max(lon0) nlat = np.min(lat0) xlat = np.max(lat0) # Thats too much! :) # if lonlatLims is not None: # print ' ' + fname + ': cutting the domain to plot !!!!' # plt.xlim(lonlatLims[0], lonlatLims[2]) # plt.ylim(lonlatLims[1], lonlatLims[3]) # print ' limits: W-E', lonlatLims[0], lonlatLims[2] # print ' limits: N-S', lonlatLims[1], lonlatLims[3] # if map_proj == 'cyl': # nlon = lonlatLims[0] # nlat = lonlatLims[1] # xlon = lonlatLims[2] # xlat = lonlatLims[3] # elif map_proj == 'lcc': # lon2 = (lonlatLims[0] + lonlatLims[2])/2. # lat2 = (lonlatLims[1] + lonlatLims[3])/2. # nlon = lonlatLims[0] # xlon = lonlatLims[2] # nlat = lonlatLims[1] # xlat = lonlatLims[3] lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) if len(dimxv.shape) == 1: lons, lats = np.meshgrid(dimxv, dimyv) else: if len(dimxv.shape) == 3: lons = dimxv[0,:,:] lats = dimyv[0,:,:] else: lons = dimxv[:] lats = dimyv[:] x,y = m(lons,lats) else: x,y = gen.lonlat2D(dimxv,dimyv) plt.rc('text', usetex=True) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vs[0], vmax=vs[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # contour ## contkind = ckind.split(',')[0] if contkind == 'cmap': cplot = plt.contour(x, y, varcv, levels=vc) elif contkind == 'fixc': plt.rcParams['contour.negative_linestyle'] = 'solid' coln = ckind.split(',')[1] cplot = plt.contour(x, y, varcv, levels=vc, colors=coln) elif contkind == 'fixsigc': coln = ckind.split(',')[1] cplot = plt.contour(x, y, varcv, levels=vc, colors=coln) else: print errormsg print ' ' + fname + ': contour kind "' + contkind + '" not defined !!!!!' quit(-1) if clabfmt is not None: plt.clabel(cplot, fmt=clabfmt) mincntS = format(vc[0], clabfmt[1:len(clabfmt)]) maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)]) else: mincntS = '{:g}'.format(vc[0]) maxcntS = '{:g}'.format(vc[len(vc)-1]) if not mapv is None: m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,xaxv[2]) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,yaxv[2]) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: plt.xlabel(dimxT) plt.ylabel(dimyT) plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) # set the limits of the plot to the limits of the data # plt.axis([x.min(), x.max(), y.min(), y.max()]) # units labels cbar.set_label(gen.latex_text(vnames[0]) + ' (' + units_lunits(uts[0]) + ')') plt.annotate(gen.latex_text(vnames[1]) +' (' + units_lunits(uts[1]) + ') [' + \ mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction', \ color=coln) figname = '2Dfields_shadow-contour' graphtit = gen.latex_text(vtit) plt.title(graphtit) output_kind(kfig, figname, ifclose) return #Nvals=50 #vals1 = np.zeros((Nvals,Nvals), dtype= np.float) #vals2 = np.zeros((Nvals,Nvals), dtype= np.float) #for j in range(Nvals): # for i in range(Nvals): # vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.) # vals2[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.) - Nvals/2 #prettylev=gen.pretty_int(-Nvals/2,Nvals/2,10) #plot_2D_shadow_contour(vals1, vals2, ['var1', 'var2'], np.arange(50)*1., \ # np.arange(50)*1., ['x-axis','y-axis'], 'rainbow', 'fixc,b', "%.2f", [0, Nvals], \ # prettylev, ['$ms^{-1}$','$kJm^{-1}s^{-1}$'], 'test var1 & var2', 'pdf', False) def plot_2D_shadow_contour_time(varsv, varcv, vnames, valv, timv, timpos, timlab, \ valu, timeu, axist, xaxv, yaxv, dimn, cbarv, ckind, clabfmt, vs, vc, uts, vtit, \ kfig, reva, ifclose): """ Adding labels and other staff to the graph varsv= 2D values to plot with shading varcv= 2D values to plot with contours vnames= variable names for the figure valv = values at the axes which is not time timv = values for the axis time timpos = positions at the axis time timlab = labes at the axis time valu = units at the axes which is not time timeu = units at the axes which is time axist = which is the axis time xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] dimn= dimension names to plot cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar= name of the color bar to use cbarfmt= format of the numbers in the colorbar cbaror= orientation of the colorbar ckind= contour kind '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 plot (None, no labels) vs= minmum and maximum values to plot in shadow vc= vector with the levels for the contour uts= units of the variable [u-shadow, u-contour] vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y ifclose= boolean value whether figure should be close (finish) or not """ fname = 'plot_2D_shadow_contour_time' if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow_contour_time.__doc__ quit() # Check variable consistency if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) if len(varcv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable contour rank:', varcv.shape, \ 'is has to be 2D!!' quit(-1) # Axis ticks # Usually time axis is the most upper dimension, thus... dimxv0 = valv.copy() dimyv0 = timv.copy() dxn = dimxv0.min() dxx = dimxv0.max() dyn = dimyv0.min() dyx = dimyv0.max() # Time ticks-labels are already given from another function if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(0.,1.+1./xaxv[2],1./xaxv[2]) dimxv0 = np.arange(0.,1.,1./valv.shape[0]) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) dimyt0 = np.array(timpos) dimxl0 = [] if xaxv[0] != 'Nfix': for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], \ style=xaxv[1])) else: for i in range(len(dimxt0)): xv=dxn+i*(dxx-dxn)/(1.*(xaxv[2])) dimxl0.append('{:{style}}'.format(xv, style=xaxv[1])) dimyl0 = list(timlab) dimxT0 = gen.latex_text(variables_values(dimn[0])[0]) +' ('+ units_lunits(valu) +\ ')' dimyT0 = timeu # Variable values usually will have time variable as 0 dimension reva = None if axist == 'x': print warnmsg print ' ' + fname + ': imposing transposition of values due to use time ' + \ 'variable at the x-axis!!' reva = 'transpose' dimxT0 = gen.latex_text(variables_values(dimn[1])[0]) + ' (' + \ units_lunits(valu) + ')' if reva is not None: varcv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varcv, reva) varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 dx=varsv.shape[1] dy=varsv.shape[0] x,y = gen.lonlat2D(dimxv,dimyv) # Changing limits of the colors vsend = graphic_range(vs,varsv) plt.rc('text', usetex=True) plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) # plt.xlim(dimxt[0],dimxt[len(dimxt)-1]) # plt.ylim(dimyt[0],dimyt[len(dimyt)-1]) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], vmax=vsend[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # contour ## contkind = ckind.split(',')[0] if contkind == 'cmap': cplot = plt.contour(x, y, varcv, levels=vc) elif contkind == 'fixc': plt.rcParams['contour.negative_linestyle'] = 'solid' coln = ckind.split(',')[1] cplot = plt.contour(x, y, varcv, levels=vc, colors=coln) elif contkind == 'fixsigc': coln = ckind.split(',')[1] cplot = plt.contour(x, y, varcv, levels=vc, colors=coln) else: print errormsg print ' ' + fname + ': contour kind "' + contkind + '" not defined !!!!!' quit(-1) if clabfmt is not None: plt.clabel(cplot, fmt=clabfmt) mincntS = format(vc[0], clabfmt[1:len(clabfmt)]) maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)]) else: mincntS = '{:g}'.format(vc[0]) maxcntS = '{:g}'.format(vc[len(vc)-1]) plt.xlabel(dimxT) plt.ylabel(dimyT) # set the limits of the plot to the limits of the data # plt.axis([dimxv.min(), dimxv.max(), dimyv.min(), dimyv.max()]) # units labels cbar.set_label(gen.latex_text(vnames[0]) + ' (' + units_lunits(uts[0]) + ')') plt.annotate(gen.latex_text(vnames[1]) +' (' + units_lunits(uts[1]) + ') [' + \ mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction', \ color=coln) figname = '2Dfields_shadow-contour' graphtit = gen.latex_text(vtit) plt.title(graphtit) output_kind(kfig, figname, ifclose) return def dxdy_lonlat(dxv,dyv,ddx,ddy): """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values dxdy_lonlat(dxv,dyv,Lv,lv) dx: values for the x dy: values for the y ddx: ',' list of which dimensions to use from values along x ddy: ',' list of which dimensions to use from values along y """ fname = 'dxdy_lonlat' if ddx.find(',') > -1: dxk = 2 ddxv = ddx.split(',') ddxy = int(ddxv[0]) ddxx = int(ddxv[1]) else: dxk = 1 ddxy = int(ddx) ddxx = int(ddx) if ddy.find(',') > -1: dyk = 2 ddyv = ddy.split(',') ddyy = int(ddyv[0]) ddyx = int(ddyv[1]) else: dyk = 1 ddyy = int(ddy) ddyx = int(ddy) ddxxv = dxv.shape[ddxx] ddxyv = dxv.shape[ddxy] ddyxv = dyv.shape[ddyx] ddyyv = dyv.shape[ddyy] slicex = [] if len(dxv.shape) > 1: for idim in range(len(dxv.shape)): if idim == ddxx or idim == ddxy: slicex.append(slice(0,dxv.shape[idim])) else: slicex.append(0) else: slicex.append(slice(0,len(dxv))) slicey = [] if len(dyv.shape) > 1: for idim in range(len(dyv.shape)): if idim == ddyx or idim == ddyy: slicey.append(slice(0,dyv.shape[idim])) else: slicey.append(0) else: slicey.append(slice(0,len(dyv))) if dxk == 2 and dyk == 2: if ddxxv != ddyxv: print errormsg print ' ' + fname + ': wrong dx dimensions! ddxx=',ddxxv,'ddyx=',ddyxv print ' choose another for x:',dxv.shape,'or y:',dyv.shape quit(-1) if ddxyv != ddyyv: print errormsg print ' ' + fname + ': wrong dy dimensions! ddxy=',ddxyv,'ddyy=',ddyv print ' choose another for x:',dxv.shape,'or y:',dyv.shape quit(-1) dx = ddxxv dy = ddxyv print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx lonv = np.zeros((dy,dx), dtype=np.float) latv = np.zeros((dy,dx), dtype=np.float) lonv = dxv[tuple(slicex)] latv = dyv[tuple(slicey)] elif dxk == 2 and dyk == 1: if not ddxxv == ddyxv and not ddxyv == ddyyv: print errormsg print ' ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv, \ 'ddyx=',ddyxv,'ddyy=',ddyyv print ' choose another for x:',dxv.shape,'or y:',dyv.shape quit(-1) dx = ddxvv dy = ddxyv print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx lonv = np.zeros((dy,dx), dtype=np.float) latv = np.zeros((dy,dx), dtype=np.float) lonv = dxv[tuple(slicex)] if ddxxv == ddyxv: for iy in range(dy): latv[iy,:] = dyv[tuple(slicey)] else: for ix in range(dx): latv[:,ix] = dyv[tuple(slicey)] elif dxk == 1 and dyk == 2: if not ddxxv == ddyxv and not ddxyv == ddyyv: print errormsg print ' ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv, \ 'ddyx=',ddyxv,'ddyy=',ddyyv print ' choose another for x:',dxv.shape,'or y:',dyv.shape quit(-1) dx = ddyxv dy = ddyyv print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx lonv = np.zeros((dy,dx), dtype=np.float) latv = np.zeros((dy,dx), dtype=np.float) latv = dyv[tuple(slicey)] if ddyxv == ddxxv: for iy in range(dy): lonv[iy,:] = dxv[tuple(slicex)] else: for ix in range(dx): lonv[:,ix] = dxv[tuple(slicex)] elif dxk == 1 and dyk == 1: dx = ddxxv dy = ddyyv # print 'dx:',dx,'dy:',dy lonv = np.zeros((dy,dx), dtype=np.float) latv = np.zeros((dy,dx), dtype=np.float) for iy in range(dy): lonv[iy,:] = dxv[tuple(slicex)] for ix in range(dx): latv[:,ix] = dyv[tuple(slicey)] return lonv,latv def dxdy_lonlatDIMS(dxv,dyv,dnx,dny,dd): """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values for a given list of values dxdy_lonlat(dxv,dyv,Lv,lv) dxv: values for the x dyv: values for the y dnx: mnames of the dimensions for values on x dny: mnames of the dimensions for values on y dd: list of [dimname]|[val] for the dimensions use [dimname]: name of the dimension [val]: value (-1 for all the range) """ fname = 'dxdy_lonlatDIMS' slicex = [] ipos=0 for dn in dnx: for idd in range(len(dd)): dname = dd[idd].split('|')[0] dvalue = dd[idd].split('|')[1] if dn == dname: if dvalue.find('@') != -1: slicex.append(slice(int(dvalue.split('@')[0]), \ int(dvalue.split('@')[1]))) else: if int(dvalue) == -1: slicex.append(slice(0,dxv.shape[ipos])) elif int(dvalue) == -9: slicex.append(dxv.shape[ipos]-1) else: slicex.append(int(dvalue)) break ipos = ipos + 1 slicey = [] ipos=0 for dn in dny: for idd in range(len(dd)): dname = dd[idd].split('|')[0] dvalue = dd[idd].split('|')[1] if dn == dname: if dvalue.find('@') != -1: slicey.append(slice(int(dvalue.split('@')[0]), \ int(dvalue.split('@')[1]))) else: if int(dvalue) == -1: slicey.append(slice(0,dyv.shape[ipos])) elif int(dvalue) == -9: slicey.append(dyv.shape[ipos]-1) else: slicey.append(int(dvalue)) break ipos = ipos + 1 lonv = dxv[tuple(slicex)] latv = dyv[tuple(slicey)] if len(lonv.shape) != len(latv.shape): print ' ' + fname + ': dimension size on x:', len(lonv.shape), 'and on y:', \ len(latv.shape),'do not coincide!!' quit(-1) return lonv,latv def plot_2D_shadow_line(varsv, varlv, vnames, vnamel, dimxv, dimyv, dimxu, dimyu, \ dimn, xaxv, yaxv, cbarv, linev, vs, uts, utl, vtit, kfig, reva, mapv, ifclose): """ Plotting a 2D field with shadows and another one with a line varsv= 2D values to plot with shading varlv= 1D values to plot with line vnames= variable names for the shadow variable in the figure vnamel= variable names for the line varibale in the figure dim[x/y]v = values at the axes of x and y dim[x/y]u = units at the axes of x and y dimn= dimension names to plot xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar: name of the color bar to use cbarfmt: format of the numbers in the colorbar cbaror: orientation of the colorbar linev= list with the values for the line [colorline, min, max, ticks: style, format, number and orientation] vs= minmum and maximum values to plot in shadow '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) uts= units of the variable to shadow utl= units of the variable to line vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= * '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', lambert conformal [res]: resolution: * 'c', crude * 'l', low * 'i', intermediate * 'h', high * 'f', full ifclose= boolean value whether figure should be close (finish) or not """ ## import matplotlib as mpl ## mpl.use('Agg') ## import matplotlib.pyplot as plt fname = 'plot_2D_shadow_line' if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow_line.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) if len(varlv.shape) != 1: print errormsg print ' ' + fname + ': wrong variable line rank:', varlv.shape, \ 'is has to be 1D!!' quit(-1) # Axis ticks # Usually axis > x must be the lon, thus... dimxv0 = dimxv.copy() dimyv0 = dimyv.copy() dxn = dimxv.min() dxx = dimxv.max() dyn = dimyv.min() dyx = dimyv.max() if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2])) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) if yaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2])) elif yaxv[0] == 'Nfix': dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2])) elif yaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,yaxv[2]) dimxl0 = [] for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], style=xaxv[1])) dimyl0 = [] for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], style=yaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')' dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')' if mapv is not None: pixkind = 'data' else: # No following data values dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0)) dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0)) pixkind = 'fixpixel' if reva is not None: varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 if len(dimxv[:].shape) == 3: xdims = '1,2' elif len(dimxv[:].shape) == 2: xdims = '0,1' elif len(dimxv[:].shape) == 1: xdims = '0' if len(dimyv[:].shape) == 3: ydims = '1,2' elif len(dimyv[:].shape) == 2: ydims = '0,1' elif len(dimyv[:].shape) == 1: ydims = '0' # lon0 = dimxv # lat0 = dimyv lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims) if not mapv is None: map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] dx = lon0.shape[1] dy = lat0.shape[0] nlon = np.min(lon0) xlon = np.max(lon0) nlat = np.min(lat0) xlat = np.max(lat0) lon2 = lon0[dy/2,dx/2] lat2 = lat0[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": map projection '" + map_proj + "' not defined!!!" print ' available: cyl, lcc' quit(-1) x,y = m(lon0,lat0) else: x = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv)) y = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv)) # Changing limits of the colors vsend = graphic_range(vs,varsv) plt.rc('text', usetex=True) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], \ vmax=vsend[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) if not mapv is None: if cbarv[0] == 'gist_gray': m.drawcoastlines(color="red") else: m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,xaxv[2]) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,xaxv[2]) m.drawparallels(parallels,labels=[False,True,True,False]) plt.xlabel('W-E') plt.ylabel('S-N') else: plt.xlabel(dimxT) plt.ylabel(dimyT) # Line ## # Line does not suffer the flipping of the y-axis it is independent to it # if reva is not None and reva[0:4] == 'flip' and reva.split('@')[1] == 'y': # b=-np.max(y)/np.max(varlv) # a=np.max(y) # else: # b=np.max(y)/np.max(varlv) # a=0. vlsend = graphic_range([linev[1],linev[2]],varlv) b=(1.-vlsend[0])/(vlsend[1]-vlsend[0]) a=0. newlinv = varlv*b+a if reva is not None and reva == 'transpose': plt.plot(newlinv, x, '-', color=linev[0], linewidth=2) else: plt.plot(x, newlinv, '-', color=linev[0], linewidth=2) plt.axis([x.min(), x.max(), y.min(), y.max()]) plt.tick_params(axis='y',right='off') if mapv is None: plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) # units labels cbar.set_label(gen.latex_text(vnames) + ' (' + units_lunits(uts) + ')') # line ticks dln = vlsend[0] dlx = vlsend[1] if linev[3] == 'pretty': dimlt = np.array(gen.pretty_int(dln,dlx,yaxv[2])) elif linev[3] == 'Nfix': dimlt = np.arange(dln,dlx,(dlx-dln)/(1.*yaxv[2])) elif linev[3] == 'Vfix': dimlt = np.arange(0,dlx,yaxv[2]) dimll = [] for i in range(len(dimlt)): dimll.append('{:{style}}'.format(dimlt[i], style=linev[4])) for it in range(len(dimlt)): yval = (dimlt[it]*b+a) plt.plot([x.max()*0.97, x.max()], [yval, yval], '-', color='k') plt.annotate(dimll[it], xy=(1.01,yval), xycoords='axes fraction',rotation=linev[6]) # plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)), \ # xycoords='axes fraction') plt.annotate(gen.latex_text(vnamel)+ ' ('+units_lunits(utl)+')', xy=(0.75,0.04), \ xycoords='figure fraction', color=linev[0]) figname = '2Dfields_shadow_line' graphtit = gen.latex_text(vtit) plt.title(graphtit) output_kind(kfig, figname, ifclose) return def plot_2D_shadow_line_time(varsv, varlv, vnames, vnamel, dimxv, dimyv, dimxu, \ dimyu, dimn, xaxv, yaxv, cbarv, linev, vs, uts, utl, vtit, kfig, reva, taxis, \ tpos, tlabs, ifclose): """ Plotting a 2D field with shadows and another one with a line evolving along time varsv= 2D values to plot with shading varlv= 1D values to plot with line vnames= variable names for the shadow variable in the figure vnamel= variable names for the line varibale in the figure dim[x/y]v = values at the axes of x and y dim[x/y]u = units at the axes of x and y dimn= dimension names to plot xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] cbarv= list with the parameters of the color bar [colorbar, cbarfmt, cbaror] colorbar: name of the color bar to use cbarfmt: format of the numbers in the colorbar cbaror: orientation of the colorbar linev= list with the values for the line [colorline, min, max, ticks: style, format, number and orientation] vs= minmum and maximum values to plot in shadow '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) uts= units of the variable to shadow utl= units of the variable to line vtit= title of the variable kfig= kind of figure (jpg, pdf, png) reva= * 'transpose': reverse the axes (x-->y, y-->x) * 'flip'@[x/y]: flip the axis x or y taxis= axis in which provide time tpos= position of the time ticks tlabs= labels of the time ticks ifclose= boolean value whether figure should be close (finish) or not """ ## import matplotlib as mpl ## mpl.use('Agg') ## import matplotlib.pyplot as plt fname = 'plot_2D_shadow_line_time' if varsv == 'h': print fname + '_____________________________________________________________' print plot_2D_shadow_line_time.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong variable shadow rank:', varsv.shape, \ 'is has to be 2D!!' quit(-1) if len(varlv.shape) != 1: print errormsg print ' ' + fname + ': wrong variable line rank:', varlv.shape, \ 'is has to be 1D!!' quit(-1) # Axis ticks dimxv0 = dimxv.copy() dimyv0 = dimyv.copy() dxn = dimxv.min() dxx = dimxv.max() dyn = dimyv.min() dyx = dimyv.max() # Definning ticks labels if taxis == 'x': dimxt0 = (tpos-np.min(tpos))/(np.max(tpos)-np.min(tpos)) dimxl0 = tlabs dimxT0 = dimxu if yaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2])) elif yaxv[0] == 'Nfix': dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2])) elif yaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,yaxv[2]) dimyl0 = [] for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], \ style=xaxv[1])) dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimyu) + ')' # No following data values dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0)) pixkind = 'fixpixel' else: if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2])) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) dimxl0 = [] for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], \ style=yaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimxu) + ')' dimyt0 = (tpos-np.min(tpos))/(np.max(tpos)-np.min(tpos)) dimyl0 = tlabs dimyT0 = dimyu # No following data values dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0)) pixkind = 'fixpixel' if reva is not None: varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 if len(dimxv[:].shape) == 3: xdims = '1,2' elif len(dimxv[:].shape) == 2: xdims = '0,1' elif len(dimxv[:].shape) == 1: xdims = '0' if len(dimyv[:].shape) == 3: ydims = '1,2' elif len(dimyv[:].shape) == 2: ydims = '0,1' elif len(dimyv[:].shape) == 1: ydims = '0' # lon0 = dimxv # lat0 = dimyv lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims) # x0 = (dimxv-np.min(dimxv))/(np.max(dimxv) - np.min(dimxv)) # y0 = (dimyv-np.min(dimyv))/(np.max(dimyv) - np.min(dimyv)) x = (lon0-np.min(lon0))/(np.max(lon0) - np.min(lon0)) y = (lat0-np.min(lat0))/(np.max(lat0) - np.min(lat0)) # Changing limits of the colors vsend = graphic_range(vs,varsv) plt.rc('text', usetex=True) plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(cbarv[0]), vmin=vsend[0], \ vmax=vsend[1]) if cbarv[2] == 'horizontal': cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) # From: http://stackoverflow.com/questions/32050030/rotation-of-colorbar-tick-labels-in-matplotlib ticklabels= cbar.ax.get_xticklabels() Nticks = len(ticklabels) ticklabs = [] for itick in range(Nticks): ticklabs.append(ticklabels[itick].get_text()) cbar.ax.set_xticklabels(ticklabs,rotation=90) else: cbar = plt.colorbar(format=cbarv[1],orientation=cbarv[2]) plt.xlabel(dimxT) plt.ylabel(dimyT) # Line ## # Line does not suffer the flipping of the y-axis it is independent to it # if reva is not None and reva[0:4] == 'flip' and reva.split('@')[1] == 'y': # b=-np.max(y)/np.max(varlv) # a=np.max(y) # else: # b=np.max(y)/np.max(varlv) # a=0. vlsend = graphic_range([linev[1],linev[2]],varlv) newlinv = (varlv-np.min(varlv))/(np.max(varlv)-np.min(varlv)) if reva is not None and reva == 'transpose': plt.plot(newlinv, y, '-', color=linev[0], linewidth=2) else: plt.plot(x, newlinv, '-', color=linev[0], linewidth=2) plt.axis([x.min(), x.max(), y.min(), y.max()]) plt.tick_params(axis='y',right='off') plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) # units labels cbar.set_label(gen.latex_text(vnames) + ' (' + units_lunits(uts) + ')') # line ticks dln = vlsend[0] dlx = vlsend[1] if linev[3] == 'pretty': dimlt = np.array(gen.pretty_int(dln,dlx,yaxv[2])) elif linev[3] == 'Nfix': dimlt = np.arange(dln,dlx,(dlx-dln)/(1.*yaxv[2])) elif linev[3] == 'Vfix': dimlt = np.arange(0,dlx,yaxv[2]) dimll = [] for i in range(len(dimlt)): dimll.append('{:{style}}'.format(dimlt[i], style=linev[4])) for it in range(1,len(dimlt)): yval = (dimlt[it]-np.min(varlv))/(np.max(varlv)-np.min(varlv)) plt.plot([x.max()*0.97, x.max()], [yval, yval], '-', color='k') plt.annotate(dimll[it], xy=(1.01,yval), xycoords='axes fraction') plt.annotate(gen.latex_text(vnamel)+ ' ('+units_lunits(utl)+')', xy=(0.75,0.04), \ xycoords='figure fraction', color=linev[0]) figname = '2Dfields_shadow_line' graphtit = gen.latex_text(vtit) plt.title(graphtit) output_kind(kfig, figname, ifclose) return def plot_Neighbourghood_evol(varsv, dxv, dyv, vnames, ttits, tpos, tlabels, colorbar, \ Nng, vs, uts, gtit, kfig, ifclose): """ Plotting neighbourghood evolution varsv= 2D values to plot with shading vnames= shading variable name for the figure d[x/y]v= values at the axes of x and y ttits= titles of both time axis tpos= positions of the time ticks tlabels= labels of the time ticks colorbar= name of the color bar to use Nng= Number of grid points of the full side of the box (odd value) vs= minmum and maximum values to plot in shadow 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) uts= units of the variable to shadow gtit= title of the graph kfig= kind of figure (jpg, pdf, png) ifclose= boolean value whether figure should be close (finish) or not """ import numpy.ma as ma fname = 'plot_Neighbourghood_evol' if varsv == 'h': print fname + '_____________________________________________________________' print plot_Neighbourghood_evol.__doc__ quit() if len(varsv.shape) != 2: print errormsg print ' ' + fname + ': wrong number of dimensions of the values: ', \ varsv.shape quit(-1) varsvmask = ma.masked_equal(varsv,fillValue) vsend = np.zeros((2), dtype=np.float) # Changing limits of the colors if type(vs[0]) != type(np.float(1.)): if vs[0] == 'Srange': vsend[0] = np.min(varsvmask) elif vs[0][0:11] == 'Saroundmean': meanv = np.mean(varsvmask) permean = np.float(vs[0].split('@')[1]) minv = np.min(varsvmask)*permean maxv = np.max(varsvmask)*permean minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)]) vsend[0] = meanv-minextrm vsend[1] = meanv+minextrm elif vs[0][0:13] == 'Saroundminmax': permean = np.float(vs[0].split('@')[1]) minv = np.min(varsvmask)*permean maxv = np.max(varsvmask)*permean vsend[0] = minv vsend[1] = maxv elif vs[0][0:17] == 'Saroundpercentile': medianv = np.median(varsvmask) valper = np.float(vs[0].split('@')[1]) minv = np.percentile(varsvmask, valper) maxv = np.percentile(varsvmask, 100.-valper) minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) vsend[0] = medianv-minextrm vsend[1] = medianv+minextrm elif vs[0][0:5] == 'Smean': meanv = np.mean(varsvmask) permean = np.float(vs[0].split('@')[1]) minv = np.min(varsvmask)*permean maxv = np.max(varsvmask)*permean minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)]) vsend[0] = -minextrm vsend[1] = minextrm elif vs[0][0:7] == 'Smedian': medianv = np.median(varsvmask) permedian = np.float(vs[0].split('@')[1]) minv = np.min(varsvmask)*permedian maxv = np.max(varsvmask)*permedian minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) vsend[0] = -minextrm vsend[1] = minextrm elif vs[0][0:11] == 'Spercentile': medianv = np.median(varsvmask) valper = np.float(vs[0].split('@')[1]) minv = np.percentile(varsvmask, valper) maxv = np.percentile(varsvmask, 100.-valper) minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)]) vsend[0] = -minextrm vsend[1] = minextrm else: print errormsg print ' ' + fname + ": range '" + vs[0] + "' not ready!!!" quit(-1) print ' ' + fname + ': modified shadow min,max:',vsend else: vsend[0] = vs[0] if type(vs[0]) != type(np.float(1.)): if vs[1] == 'range': vsend[1] = np.max(varsv) else: vsend[1] = vs[1] plt.rc('text', usetex=True) # plt.pcolormesh(dxv, dyv, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1]) plt.pcolormesh(varsvmask, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1]) cbar = plt.colorbar() newtposx = (tpos[0][:] - np.min(dxv)) * len(dxv) * Nng / (np.max(dxv) - np.min(dxv)) newtposy = (tpos[1][:] - np.min(dyv)) * len(dyv) * Nng / (np.max(dyv) - np.min(dyv)) plt.xticks(newtposx, tlabels[0]) plt.yticks(newtposy, tlabels[1]) plt.xlabel(ttits[0]) plt.ylabel(ttits[1]) plt.axes().set_aspect('equal') # From: http://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib plt.axes().xaxis.tick_top plt.axes().xaxis.set_ticks_position('top') # units labels cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')') figname = 'Neighbourghood_evol' graphtit = gtit.replace('_','\_').replace('&','\&') plt.title(graphtit, position=(0.5,1.05)) output_kind(kfig, figname, ifclose) return def plot_lines(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, gtit, gloc, gsiz, cs, \ ls, ps, ws, ss, fv, fign, kfig, ifclose): """ Function to plot a collection of lines vardv= list of set of dimension values varvv= list of set of values vaxis= which axis will be used for the values ('x', or 'y') dtit= title for the common dimension linesn= names for the legend vtit= title for the vaxis vunit= units of the vaxis gtit= main title gloc= location of the legend (0, 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' gsiz= fontsize of the legend cs= list of color names ls= list of style of lines ps= list of style of points ws= list of withs of lines ss= list of size of points fv= frequency of values fign= name of the figure kfig= kind of figure ifclose= boolean value whether figure should be close (finish) or not plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)', \ ['2.5'], 'sin', '-', 'sinus frequency dependency', 'sinus', 'pdf') """ fname = 'plot_lines' if vardv == 'h': print fname + '_____________________________________________________________' print plot_lines.__doc__ quit() # Canging line kinds every 7 lines (end of standard colors) Ntraj = len(vardv) cols, lins, pts, lws, pss = ColorsLinesPointsStyles(Ntraj, cs, ls, ps, ws, ss, fv) xmin = 100000. xmax = -100000. ymin = 100000. ymax = -100000. for il in range(Ntraj): minv = np.min(varvv[il]) maxv = np.max(varvv[il]) mind = np.min(vardv[il]) maxd = np.max(vardv[il]) if minv < xmin: xmin = minv if maxv > xmax: xmax = maxv if mind < ymin: ymin = mind if maxd > ymax: ymax = maxd plt.rc('text', usetex=True) if vaxis == 'x': for il in range(Ntraj): plt.plot(varvv[il], vardv[il], lins[il], marker=pts[il], \ linewidth=lws[il], markersize=pss[il], label= linesn[il], color=cols[il]) plt.xlabel(gen.latex_text(vtit) + ' (' + units_lunits(vunit) + ')') plt.ylabel(dtit) plt.xlim(xmin,xmax) plt.ylim(ymin,ymax) else: for il in range(Ntraj): plt.plot(vardv[il], varvv[il], lins[il], marker=pts[il], \ linewidth=lws[il], markersize=pss[il], label= linesn[il], color=cols[il]) plt.xlabel(dtit) plt.ylabel(gen.latex_text(vtit) + ' (' + units_lunits(vunit) + ')') plt.xlim(ymin,ymax) plt.ylim(xmin,xmax) graphtit = gen.latex_text(gtit) plt.title(graphtit) plt.legend(loc=gloc, prop={'size':gsiz}) output_kind(kfig, fign, ifclose) return def ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq): """ Function to provide the colors, lines & points styles Nstyles= total number of styles colors= list of colors, None for no value lines= list of lines, None for no value points= list of points, None for no value lwidths= list of line widths, None for no value psizes= list of point sizes, None for no value ptfreq= frequency for the points, None for all values >>> Nstyles = 3 >>> colors = ['blue'] >>> lines = ['-'] >>> points = ['x', '*', 'o'] >>> lwidths = ['1'] >>> psizes = ['2'] >>> ptfreq = 2 >>> colors, lines, points = ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq) ['blue', 'blue', 'blue'] ['-', '-', '-'] ['x', '*', 'o'] ['2', '2', '2'] ['1', '1', '1'] """ fname = 'ColorsLinesPointsStyles' usecolors = [] uselines = [] usepoints = [] usewlines = [] usespoints = [] # Colors if colors is None: Ncols = len(colorsauto) for ic in range(Nstyles): imc = int(ic/Ncols) iic = ic - imc * Ncols usecolors.append(colorsauto[iic]) else: Ncols = len(colors) if Ncols == 1: for ic in range(Nstyles): usecolors.append(colors[0]) else: if Ncols != Nstyles: print errormsg print ' ' + fname + ': number of provided colors:', Ncols, \ 'and required:', Nstyles,'differ !!' print ' provided:', colors quit(-1) usecolors = colors # Lines if lines is None: Nklns = len(linekindsauto) for il in range(Nstyles): iml = int(il/Nklns) iil = il - iml * Nklns uselines.append(linekindsauto[iil]) else: Nklns = len(lines) if Nklns == 1: for il in range(Nstyles): uselines.append(lines[0]) else: if Nklns != Nstyles: print errormsg print ' ' + fname + ': number of provided lines:', Nklns, \ 'and required:', Nstyles,'differ !!' print ' provided:', lines quit(-1) uselines = lines # Points if points is None: Nkpts = len(pointkindsauto) for ip in range(Nstyles): imc = int(ip/Nkpts) iip = ip - imc * Nkpts usepoints.append(pointkindsauto[iip]) else: Nkpts = len(points) if Nkpts == 1: for ip in range(Nstyles): usepoints.append(points[0]) else: if Nkpts != Nstyles: print errormsg print ' ' + fname + ': number of provided points:', Nkpts, \ 'and required:', Nstyles,'differ !!' print ' provided:', points quit(-1) usepoints = points # Line widths if lwidths is None: Nwlns = len(linewidthsauto) for il in range(Nstyles): iml = int(il/Nwlns) iil = il - iml * Nwlns usewlines.append(linewidthsauto[iil]) else: Nwlns = len(lwidths) if Nwlns == 1: for il in range(Nstyles): usewlines.append(lwidths[0]) else: if Nwlns != Nstyles: print errormsg print ' ' + fname + ': number of provided line widthss:', Nwlns, \ 'and required:', Nstyles, 'differ !!' print ' provided:', lwidths quit(-1) usewlines = lwidths # Point sizes if psizes is None: Nspts = len(pointsizesauto) for ip in range(Nstyles): ipc = int(ip/Nspts) iip = ip - ipc * Nspts usespoints.append(pointsizesauto[iip]) else: Nspts = len(psizes) if Nspts == 1: for ip in range(Nstyles): usespoints.append(psizes[0]) else: if Nspts != Nstyles: print errormsg print ' ' + fname + ': number of provided point sizes:', Nspts, \ 'and required:', Nstyles, 'differ !!' print ' provided:', psizes quit(-1) usespoints = psizes Ncols = len(usecolors) Nlins = len(uselines) Npnts = len(usepoints) lkinds = [] pkinds = [] lcolors = usecolors # Old way # if ptfreq is not None: # lkinds = uselines # pkinds = usepoints # else: # pkinds = usepoints # lkinds = uselines # for ilp in range(Nstyles): # lkinds.append(usepoints[ilp] + uselines[ilp]) return lcolors, uselines, usepoints, usewlines, usespoints #Nstyles = 3 ##colors = ['blue'] #colors = None #lines = ['-'] #points = ['x', '*', 'o'] #lwidths = ['2'] #psizes = ['1'] #ptfreq = 2 #print ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq) def LinesPointsStyles(Nstyles, lines, points, ptfreq): """ Function to provide the lines & points styles Nstyles= total number of styles lines= list of lines, None for no value points= list of points, None for no value ptfreq= frequency for the points, None for all values >>> Nstyles = 3 >>> lines = ['-'] >>> points = ['x', '*', 'o'] >>> ptfreq = 2 >>> lines, points = LinesPointsStyles(Nstyles, lines, points, ptfreq) ['-', '-', '-'] ['x', '*', 'o'] """ fname = 'LinesPointsStyles' # Canging line kinds every 7 lines (end of standard colors) uselines = [] usepoints = [] if lines is None: Nklns = len(linekindsauto) for il in range(Nstyles): iil = np.mod(il,Nklns) uselines.append(linekindsauto[iil]) else: Nklns = len(lines) if Nklns == 1: for il in range(Nstyles): uselines.append(lines[0]) else: if Nklns != Nstyles: print errormsg print ' ' + fname + ': number of provided lines:', Nklns, \ 'and required:', Nstyles,'differ !!' quit(-1) uselines = lines if points is None: Nkpts = len(pointkindsauto) for ip in range(Nstyles): iip = np.mod(ip,Nkpts) usepoints.append(pointkindsauto[iip]) else: Nkpts = len(points) if Nkpts == 1: for ip in range(Nstyles): usepoints.append(points[0]) else: if Nkpts != Nstyles: print errormsg print ' ' + fname + ': number of provided points:', Nkpts, \ 'and required:', Nstyles,'differ !!' quit(-1) usepoints = points Nlins = len(uselines) Npnts = len(usepoints) lkinds = [] pkinds = [] if ptfreq is not None: lkinds = uselines pkinds = usepoints else: pkinds = usepoints for ilp in range(Nstyles): lkinds.append(usepoints[ilp] + uselines[ilp]) return lkinds, pkinds def plot_lines_time(vardv, varvv, vaxis, dtit, linesn0, vtit, vunit, tpos, tlabs, \ gtit, gloc, gsize, kfig, vmin, vmax, lsl, coll, ptl, lwidth, psize, ptf, ifclose): """ Function to plot a collection of lines with a time axis vardv= list of set of dimension values varvv= list of set of values vaxis= which axis will be used for the time values ('x', or 'y') dtit= title for the common dimension linesn= names for the legend (None, no legend) vtit= title for the vaxis vunit= units of the vaxis tpos= positions of the time ticks tlabs= labels of the time ticks gtit= main title vmin= minimum value for the axis-value 'auto' for use the extreme vmax= maximum value for the axis-value 'auto' for use the extreme gloc= location of the legend (0, 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' gsize= fontsize of the legend kfig= kind of figure lsl= ',' list of line styles coll= ',' list of colors for the lines, None for automatic, single value all the same ptl= ',' list of type of points for the lines, None for automatic, single value all the same lwidth= ',' list of line widths psize= ',' list of point sizes ptf= frequency of point plotting, 'all' for all time steps ifclose= boolean value whether figure should be close (finish) or not plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)', \ ['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf') """ fname = 'plot_lines_time' if vardv == 'h': print fname + '_____________________________________________________________' print plot_lines_time.__doc__ quit() Ntraj = len(vardv) colvalues, linekinds, pointkinds, lwidths, psizes = ColorsLinesPointsStyles( \ Ntraj, coll, lsl, ptl, lwidth, psize, ptf) plt.rc('text', usetex=True) xtrmvv = [fillValueF,-fillValueF] xtrmdv = [fillValueF,-fillValueF] # Do we have legend? ## if linesn0 is None: linesn = [] for itrj in range(Ntraj): linesn.append(str(itrj)) else: linesn = linesn0 if vaxis == 'x': for il in range(Ntraj): plt.plot(varvv[il], vardv[il], linekinds[il], marker=pointkinds[il], \ label=linesn[il], color=colvalues[il], linewidth=lwidths[il], \ markersize=psizes[il], markevery=ptf) # Old way # if ptf is None: # plt.plot(varvv[il], vardv[il], linekinds[il], label= linesn[il], \ # color=colvalues[il], linewidth=lwidths[il], markersize=psizes[il]) # else: # plt.plot(varvv[il], vardv[il], linekinds[il], color=colvalues[il], \ # linewidth=lwidths[il]) # plt.plot(varvv[il][::ptf], vardv[il][::ptf], pointkinds[il], \ # label= linesn[il], color=colvalues[il], linewidth=lwidths[il], \ # markersize=psizes[il]) minvv = np.min(varvv[il]) maxvv = np.max(varvv[il]) mindv = np.min(vardv[il]) maxdv = np.max(vardv[il]) if minvv < xtrmvv[0]: xtrmvv[0] = minvv if maxvv > xtrmvv[1]: xtrmvv[1] = maxvv if mindv < xtrmdv[0]: xtrmdv[0] = mindv if maxdv > xtrmdv[1]: xtrmdv[1] = maxdv plt.xlabel(vtit + ' (' + vunit + ')') plt.ylabel(dtit) # plt.xlim(np.min(varTvv),np.max(varTvv)) # plt.ylim(np.min(varTdv),np.max(varTdv)) if vmin == 'auto': xmin = xtrmvv[0] else: xmin = vmin if vmax == 'auto': xmax = xtrmvv[1] else: xmax = vmax plt.xlim(xmin,xmax) plt.ylim(xtrmdv[0],xtrmdv[1]) plt.yticks(tpos, tlabs) else: for il in range(Ntraj): Nvals = len(vardv[il]) plt.plot(vardv[il], varvv[il], linekinds[il], marker=pointkinds[il], \ label=linesn[il], color=colvalues[il], linewidth=lwidths[il], \ markersize=psizes[il], markevery=ptf) # Old way # if ptf is None: # plt.plot(vardv[il], varvv[il], linekinds[il], label= linesn[il], \ # color=colvalues[il], linewidth=lwidths[il], markersize=psizes[il], markevery=ptf) # else: # plt.plot(vardv[il], varvv[il], linekinds[il], color=coll[il], \ # linewidth=lwidths[il]) # plt.plot(vardv[il][::ptf], varvv[il][::ptf], pointkinds[il], \ # label= linesn[il], color=colvalues[il], linewidth=lwidths[il], \ # markersize=psizes[il]) minvv = np.min(varvv[il]) maxvv = np.max(varvv[il]) mindv = np.min(vardv[il]) maxdv = np.max(vardv[il]) if minvv < xtrmvv[0]: xtrmvv[0] = minvv if maxvv > xtrmvv[1]: xtrmvv[1] = maxvv if mindv < xtrmdv[0]: xtrmdv[0] = mindv if maxdv > xtrmdv[1]: xtrmdv[1] = maxdv plt.xlabel(dtit) plt.ylabel(vtit + ' (' + vunit + ')') if vmin == 'auto': ymin = xtrmvv[0] else: ymin = vmin if vmax == 'auto': ymax = xtrmvv[1] else: ymax = vmax plt.xlim(xtrmdv[0],xtrmdv[1]) plt.ylim(ymin,ymax) # plt.xlim(np.min(varTdv),np.max(varTdv)) # plt.ylim(np.min(varTvv),np.max(varTvv)) plt.xticks(tpos, tlabs) figname = 'lines_time' graphtit = gtit.replace('_','\_').replace('&','\&') plt.title(graphtit) if linesn0 is not None: if Ntraj < 10: plt.legend(loc=gloc, prop={'size':gsize}) elif 10 < Ntraj < 20: plt.legend(loc=gloc, prop={'size':int(gsize*0.8)}) else: plt.legend(loc=gloc, prop={'size':int(gsize*0.6)}) print plt.xlim(),':', plt.ylim() output_kind(kfig, figname, ifclose) return def plot_barbs(xvals, yvals, uvals, vvals, vecfreq, veccolor, barblength, windn, \ wuts, xaxv, yaxv, dimn, dimu, mapv, reva, graphtit, kfig, figname, ifclose): """ Function to plot wind barbs xvals= x position of the values yvals= y position of the values uvals= values for the x-wind uvals= values for the y-wind vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 'auto', computed automatically to have 20 vectors along each axis) veccolor= color of the vectors ('auto', for 'red') two options: [colorname]: name of the color fixed for all vectors 'colormap',[colormapname]: use colormap to provide the colors tacking wind speed as reference barblength= length of the wind barbs ('auto', for 9) windn= name of the wind variable in the graph wuts= units of the wind variable in the graph xaxv= list with the x-axis paramteres [style, format, number and orientation] yaxv= list with the y-axis paramteres [style, format, number and orientation] dimn= name of the X,Y variable-dimensions dimu= units of the X,Y variable-dimensions mapv= 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 graphtit= title of the graph ('|', for spaces) kfig= kind of figure figname= name of the figure ifclose= boolean value whether figure should be close (finish) or not """ fname = 'plot_barbs' if len(xvals.shape) != 2: print errormsg print ' ' + fname + ': wrong variable xvals rank:', xvals.shape, \ 'is has to be 2D!!' quit(-1) if len(yvals.shape) != 2: print errormsg print ' ' + fname + ': wrong variable yvals rank:', yvals.shape, \ 'is has to be 2D!!' quit(-1) if len(uvals.shape) != 2: print errormsg print ' ' + fname + ': wrong variable uvals rank:', uvals.shape, \ 'is has to be 2D!!' quit(-1) if len(vvals.shape) != 2: print errormsg print ' ' + fname + ': wrong variable vvals rank:', vvals.shape, \ 'is has to be 2D!!' quit(-1) dx=xvals.shape[1] dy=xvals.shape[0] # Axis ticks # Usually axis > x must be the lon, thus... if mapv is None: dimxv0 = xvals[0,:] dimyv0 = yvals[:,0] else: dimxv0 = xvals.copy() dimyv0 = yvals.copy() dxn = dimxv0.min() dxx = dimxv0.max() dyn = dimyv0.min() dyx = dimyv0.max() if xaxv[0] == 'pretty': dimxt0 = np.array(gen.pretty_int(dxn,dxx,xaxv[2])) elif xaxv[0] == 'Nfix': dimxt0 = np.arange(dxn,dxx,(dxx-dxn)/(1.*xaxv[2])) elif xaxv[0] == 'Vfix': dimxt0 = np.arange(0,dxx,xaxv[2]) if yaxv[0] == 'pretty': dimyt0 = np.array(gen.pretty_int(dyn,dyx,yaxv[2])) elif yaxv[0] == 'Nfix': dimyt0 = np.arange(dyn,dyx,(dyx-dyn)/(1.*yaxv[2])) elif yaxv[0] == 'Vfix': dimyt0 = np.arange(0,dyx,yaxv[2]) dimxl0 = [] for i in range(len(dimxt0)): dimxl0.append('{:{style}}'.format(dimxt0[i], \ style=xaxv[1])) dimyl0 = [] for i in range(len(dimyt0)): dimyl0.append('{:{style}}'.format(dimyt0[i], \ style=yaxv[1])) dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(dimu[0]) + ')' dimyT0 = variables_values(dimn[1])[0] + ' (' + units_lunits(dimu[1]) + ')' if mapv is not None: pixkind = 'data' else: # No following data values # dimxt0 = np.arange(len(dimxt0),dtype=np.float)/(len(dimxt0)) # dimyt0 = np.arange(len(dimyt0),dtype=np.float)/(len(dimyt0)) pixkind = 'fixpixel' if reva is not None: uvals, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(uvals, reva) vvals, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT = \ transform(vvals, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0, \ dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0) else: dimxv = dimxv0 dimyv = dimyv0 dimxt = dimxt0 dimyt = dimyt0 dimxl = dimxl0 dimyl = dimyl0 dimxT = dimxT0 dimyT = dimyT0 # Frequency of vectors if vecfreq is None: xfreq = 1 yfreq = 1 elif vecfreq == 'auto': xfreq = dx/20 yfreq = dy/20 else: xfreq=int(vecfreq.split('@')[0]) yfreq=int(vecfreq.split('@')[1]) if veccolor == 'auto': bctype = 'fix' vcolor = "red" else: if veccolor.find('@') != -1: vcvals = veccolor.split('@') if vcvals[0] == 'colormap': wscol = [] bctype = 'cmap' cmap = vcvals[1] wspeed = np.sqrt(uvals[::yfreq,::xfreq]*uvals[::yfreq,::xfreq] + \ vvals[::yfreq,::xfreq]*vvals[::yfreq,::xfreq]) else: print errormsg print ' '+fname+ ": barb color type '" + vcvals[0] + "' not redy !!" else: bctype = 'fix' vcolor = veccolor if barblength == 'auto': blength = 9 else: blength = np.float(barblength) plt.rc('text', usetex=True) if not mapv is None: lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:]) lat00 = yvals[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = np.min(xvals[::yfreq,::xfreq]) xlon = np.max(xvals[::yfreq,::xfreq]) nlat = np.min(yvals[::yfreq,::xfreq]) xlat = np.max(yvals[::yfreq,::xfreq]) lon2 = xvals[dy/2,dx/2] lat2 = yvals[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') else: plt.xlabel(dimxT) plt.ylabel(dimyT) # plt.axis([x.min(), x.max(), y.min(), y.max()]) if mapv is None: plt.xticks(dimxt, dimxl, rotation=xaxv[3]) plt.yticks(dimyt, dimyl, rotation=yaxv[3]) if bctype == 'fix': plt.barbs(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], color=vcolor, pivot='tip', \ length=blength) plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ xy=(0.85,-0.10), xycoords='axes fraction', color=vcolor) elif bctype == 'cmap': plt.barbs(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], wspeed, pivot='tip', \ length=blength, cmap=plt.get_cmap(cmap)) plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ xy=(0.85,-0.10), xycoords='axes fraction', color='k') # cbar = plt.colorbar() plt.title(gen.latex_text(graphtit.replace('|',' '))) ## NOT WORKING ## # No legend so it is imposed ## windlabel=windn.replace('_','\_') +' (' + units_lunits(wuts[1]) + ')' ## vecpatch = mpatches.Patch(color=vcolor, label=windlabel) ## plt.legend(handles=[vecpatch]) ## vecline = mlines.Line2D([], [], color=vcolor, marker='.', markersize=10, label=windlabel) ## plt.legend(handles=[vecline], loc=1) output_kind(kfig, figname, ifclose) return def plot_ptZvals(vname,vunits,points,ptype,ptsize,graphlims,minmax,figtitle,cbar, \ mapv,kfig): """ Function to plot a given list of points and values vname= name of the variable in the graph vunits= units of the variable points= [lon,lat,val] matrix of values ptype= type of the point ptsize= size of the point graphlims= minLON,minLAT,maxLON,maxLAT limits of the graph, None for the full size minmax= minimum and maximum type 'auto': values taken from the extrems of the data [min],[max]: given minimum and maximum values figtitle= title of the figure cbar= color bar mapv= 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 kfig= kind of figure """ fname = 'plot_ptZvals' figname = 'pointsZval' minlon = points[:,0].min() maxlon = points[:,0].max() minlat = points[:,1].min() maxlat = points[:,1].max() minval = points[:,2].min() maxval = points[:,2].max() # print 'min/max val;',minval,maxval lonrange = (points[:,0] - minlon)/(maxlon - minlon) latrange = (points[:,1] - minlat)/(maxlat - minlat) colorrange = (points[:,2] - minval)/(maxval - minval) plt.rc('text', usetex=True) if mapv is not None: vlon = points[:,0] vlat = points[:,1] dx = len(vlon) dy = len(vlat) # vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:]) # xvala = np.array(xval) # xvala = np.where(xvala < 0., 360. + xvala, xvala) # xval = list(xvala) map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] if graphlims is not None: nlon = graphlims[0] xlon = graphlims[2] nlat = graphlims[1] xlat = graphlims[3] else: nlon = np.min(vlon) xlon = np.max(vlon) nlat = np.min(vlat) xlat = np.max(vlat) lon2 = vlon[dy/2] lat2 = vlat[dy/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": map projecion '" + map_proj + "' not ready!!" print ' available: cyl, lcc' quit(-1) # lons, lats = np.meshgrid(vlon, vlat) # lons = np.where(lons < 0., lons + 360., lons) x,y = m(vlon,vlat) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True]) parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False]) # else: # x = vlon # y = vlat # plt.xlim(0,dx-1) # plt.ylim(0,dy-1) if minmax == 'auto': plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar, \ marker=ptype) else: minv = np.float(minmax.split(',')[0]) maxv = np.float(minmax.split(',')[1]) plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar, \ marker=ptype, vmin=minv, vmax=maxv) cbar = plt.colorbar() cbar.set_label(vname.replace('_','\_') +' ('+ units_lunits(vunits) + ')') plt.title(figtitle) if graphlims is not None: plt.xlim(graphlims[0], graphlims[2]) plt.ylim(graphlims[1], graphlims[3]) output_kind(kfig, figname, True) return #pts = np.zeros((10,3), dtype=np.float) #pts[:,0] = np.arange(10,20)*1. #pts[:,1] = np.arange(30,40)*1. #pts[:,2] = np.arange(-5,5)*1. #plot_ptZvals('vals','kgm-2',pts,'.',300, 'values of values', 'seismic', 'cyl,l', 'pdf') def plot_ZQradii(Zmeans, graphtit, kfig, figname): """ Function to plot following radial averages only at exact grid poins Zmeans= radial means radii= values of the taken radii graphtit= title of the graph ('|', for spaces) kfig= kind of figure figname= name of the figure """ fname = 'plot_ZQradii' output_kind(kfig, figname, True) return def plot_vector(xvals,yvals,uvals,vvals,vecfreq,vecoln,veccolor,veclength,windn,wuts,\ mapv,graphtit,kfig,figname): """ Function to plot vectors xvals= values for the x-axis yvals= values for the y-axis uvals= values for the x-wind vvals= values for the y-wind vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 'auto', computed automatically to have 20 vectors along each axis) veccoln= name for the color of the vectors 'singlecol'@[colorn]: all the vectors same color ('auto': for 'red') 'wind'@[colorbar]: color of the vectors according to wind speed sqrt(u^2+v^2) and given [colorbar] '3rdvar'@[colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar] veccolor= color of the vectors veclength= length of the wind vectors: 'singlecol': 'auto', for 9 'wind' and '3rdvar': 'auto' length as wind speed, otherwise fix length windn= name of the wind variable in the graph wuts= units of the wind variable in the graph mapv= 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 graphtit= title of the graph ('|', for spaces) kfig= kind of figure figname= name of the figure """ fname = 'plot_vector' dx=xvals.shape[1] dy=xvals.shape[0] # Frequency of vectors if vecfreq is None: xfreq = 1 yfreq = 1 elif vecfreq == 'auto': xfreq = dx/20 yfreq = dy/20 else: xfreq=int(vecfreq.split('@')[0]) yfreq=int(vecfreq.split('@')[1]) # Vector length if veclength == 'auto': vlength = 9 else: vlength = veclength # Colors VecN = vecoln.split('@')[0] if VecN == 'singlecol': if veccolor == 'auto': vcolor = "red" else: vcolor = veccolor elif VecN == 'wind' or VecN == '3rdvar': vcolor = vecoln.split('@')[1] plt.rc('text', usetex=True) if not mapv is None: lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:]) lat00 = yvals[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = np.min(xvals[::yfreq,::xfreq]) xlon = np.max(xvals[::yfreq,::xfreq]) nlat = np.min(yvals[::yfreq,::xfreq]) xlat = np.max(yvals[::yfreq,::xfreq]) lon2 = xvals[dy/2,dx/2] lat2 = yvals[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') if VecN == 'singlecol': plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], color=vcolor, pivot='middle') plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ xy=(0.80,-0.15), xycoords='axes fraction', color=vcolor) else: if veclength != 'auto': wind = np.sqrt(uvals**2 + vvals**2) uvals = np.where(wind == 0., 0., uvals) vvals = np.where(wind == 0., 0., vvals) uvals = np.float(veclength)*uvals/wind vvals = np.float(veclength)*vvals/wind plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], veccolor[::yfreq,::xfreq], cmap=plt.get_cmap(vcolor), pivot='middle') cbar = plt.colorbar() if VecN == 'wind': cbar.set_label('$\sqrt{u^{2} + v^{2}}$ (' + units_lunits(wuts) + ')') else: vN = vecoln.split('@')[2] vU = vecoln.split('@')[3] cbar.set_label(vN + ' (' + units_lunits(vU) + ')') plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ xy=(0.80,-0.15), xycoords='axes fraction', color='black') plt.title(graphtit.replace('|',' ').replace('&','\&')) output_kind(kfig, figname, True) return def plot_basins(xvals,yvals,fvals,vecfreq,vecoln,veccolor,veclength,vcolmin,vcolmax,windn,wuts,\ mapv,drawcountry,basinid,graphtit,kfig,figname): """ Function to plot vectors xvals= values for the x-axis yvals= values for the y-axis fvals= values for the flow (1-8: N, NE, E, ..., NW; 97: sub-basin; 98: small to sea; 99: large to sea) vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 'auto', computed automatically to have 20 vectors along each axis) veccoln= name for the color of the vectors (as '3rdvar@[basinsvar]@[varn]@[units]' from plot_vector) veccolor= color of the vectors veclength= length of the wind vectors: 'singlecol': 'auto', for 9 'wind' and '3rdvar': 'auto' length as wind speed, otherwise fix length windn= name of the wind variable in the graph wuts= units of the wind variable in the graph mapv= 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 drawcountry= whether country lines should be plotted or not basinid= whether basins id should be plotted or not graphtit= title of the graph ('|', for spaces) kfig= kind of figure figname= name of the figure """ fname = 'plot_basins' dx=xvals.shape[1] dy=xvals.shape[0] # Frequency of vectors if vecfreq is None: xfreq = 1 yfreq = 1 elif vecfreq == 'auto': xfreq = dx/20 yfreq = dy/20 else: xfreq=int(vecfreq.split('@')[0]) yfreq=int(vecfreq.split('@')[1]) # Vector length if veclength == 'auto': vlength = 9 else: vlength = veclength # flow direction angle = (fvals[::yfreq,::xfreq] - 1)*np.pi/4 uvals = np.where(fvals[::yfreq,::xfreq] < 9, np.float(veclength)*np.sin(angle), 0.) vvals = np.where(fvals[::yfreq,::xfreq] < 9, np.float(veclength)*np.cos(angle), 0.) # Colors vcolor = vecoln.split('@')[0] plt.rc('text', usetex=True) ddx = np.abs(xvals[dy/2+1,dx/2] - xvals[dy/2,dx/2]) ddy = np.abs(yvals[dy/2,dx/2+1] - yvals[dy/2,dx/2]) fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3} # Setting up colors for each label # From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar my_cmap = plt.cm.get_cmap(vcolor) # vcolmin = np.min(veccolor[::yfreq,::xfreq]) # vcolmax = np.max(veccolor[::yfreq,::xfreq]) # vcolmin = 0. # vcolmax = 2500. norm = mpl.colors.Normalize(vcolmin, vcolmax) print 'min col:',vcolmin,'max col:',vcolmax xlabpos = [] ylabpos = [] labels = [] labcol = [] flow = [] flowvals = [] for j in range(0,dy,yfreq): for i in range(0,dx,xfreq): if veccolor[j,i] != '--': xlabpos.append(xvals[j,i]) ylabpos.append(yvals[j,i]) labels.append(int(veccolor[j,i])) labcol.append(my_cmap(norm(veccolor[j,i]))) flowvals.append(fvals[j,i]) if not mapv is None: lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:]) lat00 = yvals[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = np.min(xvals[::yfreq,::xfreq]) xlon = np.max(xvals[::yfreq,::xfreq]) nlat = np.min(yvals[::yfreq,::xfreq]) xlat = np.max(yvals[::yfreq,::xfreq]) lon2 = xvals[dy/2,dx/2] lat2 = yvals[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() drawcountry = True if drawcountry: m.drawcountries() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') if veclength != 'auto': wind = np.sqrt(uvals**2 + vvals**2) uvals = np.where(wind == 0., 0., uvals) vvals = np.where(wind == 0., 0., vvals) uvals = np.float(veclength)*uvals/wind vvals = np.float(veclength)*vvals/wind vecolorvals = veccolor[::yfreq,::xfreq] vecolorvals = np.where(vecolorvals < vcolmin, vcolmin - 1, vecolorvals) vecolorvals = np.where(vecolorvals > vcolmax, vcolmax + 1, vecolorvals) # plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ # uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], veccolor[::yfreq,::xfreq], \ # cmap=plt.get_cmap(vcolor), pivot='middle') plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], \ uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], vecolorvals, \ cmap=plt.get_cmap(vcolor), pivot='middle') cbar = plt.colorbar() vN = vecoln.split('@')[1] vU = vecoln.split('@')[2] cbar.set_label(vN + ' (' + units_lunits(vU) + ')') if basinid: for i in range(len(xlabpos)): plt.text(xlabpos[i]+0.5*ddx, ylabpos[i]+0.95*ddy, labels[i], \ color=labcol[i], fontdict=fontcharac) # Sea-flow for i in range(len(xlabpos)): if flowvals[i] == 97: plt.plot(xlabpos[i], ylabpos[i], 'x', color=labcol[i]) elif flowvals[i] == 98: plt.plot(xlabpos[i], ylabpos[i], '*', color=labcol[i]) elif flowvals[i] == 99: plt.plot(xlabpos[i], ylabpos[i], 'h', color=labcol[i]) plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ xy=(0.80,-0.15), xycoords='axes fraction', color='black') plt.title(graphtit.replace('|',' ').replace('&','\&')) output_kind(kfig, figname, True) return def plot_river_pattern(xvals, yvals, fvals, veccolor, descid, desclon, desclat, \ stn, stlon, stlat, drawcountry, drawbasinid, colordescid, coloriver, colornoriver, \ colorst, stsize, graphtit, mapv, kfig, figname): """ Function to plot vectors xvals= values for the x-axis yvals= values for the y-axis fvals= values for the flow (1-8: N, NE, E, ..., NW; 97: sub-basin; 98: small to sea; 99: large to sea) veccolor= basin id descid= description id , desclon, desclat, mapv= 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 drawcountry= whether country lines should be plotted or not drawbasinid= whether basins id should be plotted or not graphtit= title of the graph ('|', for spaces) kfig= kind of figure figname= name of the figure """ fname = 'plot_river_pattern' dx=xvals.shape[1] dy=xvals.shape[0] xmin = np.min(xvals) xmax = np.max(xvals) ymin = np.min(yvals) ymax = np.max(yvals) # flow direction plt.rc('text', usetex=True) ddx = np.abs(xvals[dy/2+1,dx/2] - xvals[dy/2,dx/2]) ddy = np.abs(yvals[dy/2,dx/2+1] - yvals[dy/2,dx/2]) fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3} xlabpos = [] ylabpos = [] labels = [] labcol = [] flow = [] flowvals = [] xtrack = [] ytrack = [] colortrack = [] trackid = [] stnlL = [] stlonlL = [] stlatlL = [] # Setting up colors for each label # From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar my_cmap = plt.cm.get_cmap(colordescid) vcolmin = np.min(descid) vcolmax = np.max(descid) norm = mpl.colors.Normalize(vcolmin, vcolmax) # Vector angles lengthtrac = np.float(xvals[0,1] - xvals[0,0]) lengthtrac2 = lengthtrac*np.sqrt(2.) for j in range(0,dy): for i in range(0,dx): if veccolor[j,i] != '--': # ibeg = xvals[j,i]-lengthtrac/2. # jbeg = yvals[j,i]-lengthtrac/2. ibeg = xvals[j,i] jbeg = yvals[j,i] labels.append(int(veccolor[j,i])) flowvals.append(fvals[j,i]) angle = (fvals[j,i] - 1)*np.pi/4 if gen.searchInlist([2,4,6,8], fvals[j,i]): iend = ibeg + lengthtrac2*np.sin(angle) jend = jbeg + lengthtrac2*np.cos(angle) elif gen.searchInlist([1,3,5,7], fvals[j,i]): iend = ibeg + lengthtrac*np.sin(angle) jend = jbeg + lengthtrac*np.cos(angle) else: ibeg = xvals[j,i] jbeg = yvals[j,i] iend = None jend = None xlabpos.append(ibeg) ylabpos.append(jbeg) xtrack.append(ibeg) xtrack.append(iend) xtrack.append(None) ytrack.append(jbeg) ytrack.append(jend) ytrack.append(None) if len(desclon.shape) == 2: difflonlat = np.sqrt((desclon-xvals[j,i])**2 + (desclat-yvals[j,i])**2) mindiffLl = np.min(difflonlat) ilatlon = gen.index_mat(difflonlat, mindiffLl) else: ilatlon = np.zeros((2), dtype=int) difflon = np.abs(desclon - xvals[j,i]) mindiffl = np.min(difflon) ilatlon[1] = gen.index_vec(difflon,mindiffl) difflat = np.abs(desclat - yvals[j,i]) mindiffL = np.min(difflat) ilatlon[0] = gen.index_vec(difflat,mindiffL) if descid.mask[ilatlon[0],ilatlon[1]]: labcol.append(colornoriver) colortrack.append(colornoriver) else: if veccolor[j,i] != 6: labcol.append(colornoriver) colortrack.append(colornoriver) # print 'Lluis: veccol:', veccolor[j,i], 'mindiffl', mindiffl,'mindiffL:',mindiffL else: # print 'Lluis right! mindiffl', mindiffl,'mindiffL:',mindiffL labcol.append(coloriver) colortrack.append(my_cmap(norm(descid[ilatlon[0],ilatlon[1]]))) trackid.append(descid[ilatlon[0],ilatlon[1]]) totvals = len(flowvals) if not mapv is None: lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:]) lat00 = yvals[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] nlon = np.min(xvals) xlon = np.max(xvals) nlat = np.min(yvals) xlat = np.max(yvals) lon2 = xvals[dy/2,dx/2] lat2 = yvals[dy/2,dx/2] print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() drawcountry = True if drawcountry: m.drawcountries() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') j = 0 for i in range(len(xlabpos)): plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i]) j = j + 3 # Sea-flow for i in range(len(xlabpos)): if flowvals[i] == 97: plt.plot(xlabpos[i], ylabpos[i], 'x', color=colortrack[i]) elif flowvals[i] == 98: plt.plot(xlabpos[i], ylabpos[i], '*', color=colortrack[i]) elif flowvals[i] == 99: plt.plot(xlabpos[i], ylabpos[i], 'h', color=colortrack[i]) plt.xlim(xmin,xmax) plt.ylim(ymin,ymax) if drawbasinid: for i in range(len(xlabpos)): plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, labels[i], \ color=labcol[i], fontdict=fontcharac) # plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, trackid[i], \ # color=colortrack[i], fontdict=fontcharac) # plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, str(flowvals[i]), \ # color=colortrack[i], fontdict=fontcharac) for ist in range(len(stlon)): if stlon[ist] > xmin and stlon[ist] < xmax and stlat[ist] > ymin and stlat[ist] < ymax: stname = '*' for ic in range(len(stn[ist,:])): stname = stname + stn[ist,ic] # plt.annotate(stname, xy=(stlon[ist],stlat[ist]), color=colorst, ) plt.text(stlon[ist], stlat[ist], stname, color=colorst, fontsize=stsize) # cbar = plt.colorbar() plt.title(graphtit.replace('&','\&')) output_kind(kfig, figname, True) return def var_3desc(ovar): """ Function to provide std_name, long_name and units from an object variable ovar= object variable """ fname = 'var_desc' varattrs = ovar.ncattrs() if searchInlist(varattrs,'std_name'): stdn = ovar.getncattr('std_name') else: vvalues = variables_values(ovar._name) stdn = vvalues[1] if searchInlist(varattrs,'long_name'): lonn = ovar.getncattr('long_name') else: lonn = vvalues[4].replace('|',' ') if searchInlist(varattrs,'units'): un = ovar.getncattr('units') else: un = vvalues[5] return stdn, lonn, un def plot_river_desc(lons, lats, rns, rss, rus, ros, bcolor, ucolor, uuts, mapv, \ graphtit, kfig, lloc, figname): """ Function to plot rivers from 'river_desc.nc' ORCDHIEE lons= values for the x-axis lats= values for the y-axis rns= list of the name of the rivers to plot rss= dictionary with the lon,lats of the subbasins of each river rus= dictionary with the lon,lats of the upstream values of each river ros= dictionary with the lon,lats of the outflow points of each river bcolor= color of the lines for the subbasins ucolor= bar color for the upstreams uuts= units of the upstream mapv= 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 graphtit= title of the graph ('|', for spaces) lloc= location of the legend (0, automatic) 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' kfig= kind of figure figname= name of the figure """ fname = 'plot_basins' colvalues, linekinds, pointkinds, lwidths, psizes = ColorsLinesPointsStyles( \ len(rns), bcolor, ['-'], ['.'], [1.], [1.], None) dx=lons.shape[1] dy=lats.shape[0] plt.rc('text', usetex=True) nlon = np.min(lons) xlon = np.max(lons) nlat = np.min(lats) xlat = np.max(lats) dlon = xlon - nlon dlat = xlat - nlat # Making bigger the area to map nlon = np.min(lons)-dlon*0.1 xlon = np.max(lons)+dlon*0.1 nlat = np.min(lats)-dlat*0.1 xlat = np.max(lats)+dlat*0.1 plt.xlim(nlon,xlon) plt.ylim(nlat,xlat) if not mapv is None: lon00 = np.where(lons[:] < 0., 360. + lons[:], lons[:]) lat00 = lats[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] lon2 = (nlon + xlon)/2. lat2 = (nlat + xlat)/2. print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') umin = 10.e20 umax = -umin for rn in rns: nrns = np.min(rus[rn]) xrns = np.max(rus[rn]) if umin > nrns: umin = nrns if umax < xrns: umax = xrns i = 0 for rn in rns: plt.pcolormesh(lons, lats, rus[rn], vmin=umin, vmax=umax, \ cmap=plt.get_cmap(ucolor), norm=matplotlib.colors.LogNorm()) # plt.pcolormesh(lons, lats, rss[rn], cmap=plt.get_cmap(ucolor)) if rn == rns[0]: cbar = plt.colorbar() riverarea = rss[rn] riverarea.mask = ma.nomask rarea = np.where(riverarea == 999999999, -5, riverarea) rarea = np.where(rarea > 0, 1, 0) plt.contour(lons, lats, rarea, levels = [0.], colors=colvalues[i], \ linewidths=1.5) oflow = ros[rn] plt.plot(oflow[0], oflow[1], 'h', color=colvalues[i], label=rn) # plt.annotate(rn, xy=(0., i*0.05), xycoords='axes fraction', color=colvalues[i]) i = i+1 # Attempts to plot sub-basins # Gradient subbasins (not working) # gradriver = np.zeros((dy,dx), dtype=int) # for j in range(dy-2): # for i in range(dx-2): # gradriver[j,i] = riverarea[j,i+1]-riverarea[j,i]+riverarea[j+1,i]- \ # riverarea[j,i] # garea = np.where(gradriver != 0, 1, gradriver) # print garea # plt.contour(lons, lats, garea, levels = [0.], colors='red', linewidths=1.) # Not working # subbasins = np.sort(np.unique(rarea)) # for subb in subbasins: # if subb > 0.: # sarea = np.where(rarea != subb, 0, rarea) # plt.contour(lons, lats, sarea, levels = [1.], colors='gray', linewidths=1.) cbar.set_label('upstream (' + units_lunits(uuts) + ')') plt.legend(loc=lloc) # plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')', \ # xy=(0.80,-0.15), xycoords='axes fraction', color='black') plt.title(graphtit.replace('|',' ').replace('&','\&')) output_kind(kfig, figname, True) return def plot_vertical_lev(vertz, vertp, zlog, dzlog, plog, dplog, gtit, kfig, lloc): """ plotting vertical levels distribution plot_vertical_lev(vertz, gtit, kfig, lloc) vertz= distribution of vertical heights [z] vertp= distribution of vertical pressures [Pa] zlog: to use logarithmic scale on the height axis ('true/false') dzlog: to use logarithmic scale on the difference of height between levels axis ('true/false') plog: to use logarithmic scale on the pressure axis ('true/false') pzlog: to use logarithmic scale on the difference of pressure between levels axis ('true/false') gtit= title of the graph ('!' for spaces) kfig= kind of figure (jpg, pdf, png) lloc= location of the legend (0, automatic) 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' kfig= kind of figure """ fname = 'plot_vertical_lev' figname = 'vertical_lev' if vertz == 'h': print fname + '_____________________________________________________________' print plot_vertical_lev.__doc__ quit() Nlev = len(vertz) dvertz = vertz[1:Nlev] - vertz[0:Nlev-1] # From: http://stackoverflow.com/questions/14762181/adding-a-y-axis-label-to-secondary-y-axis-in-matplotlib plt.rc('text', usetex=True) # Height plot ## if vertp is not None: print plt.subplots.__doc__ fig, (ax1, ax3) = plt.subplots(nrows=2, ncols=1, sharex=True) else: fig, ax1 = plt.subplots() ax2 = ax1.twinx() plt.xlim(0,Nlev) l1 = ax1.plot(range(Nlev), vertz, 'r-x', label='height') l2 = ax2.plot(range(1,Nlev), dvertz, 'b-x', label='dheight') # ax1.set_xlabel('level (\#)') ax1.set_ylabel('height (m)', color='r') ax1.set_ylim(1,np.max(vertz)) ax1.grid() ax2.set_ylabel('difference between levels (m)', color='b') if vertp is not None: ax1.set_title('height') if zlog: ax1.set_yscale('log') if dzlog: ax2.set_yscale('log') # Pressure plot ## if vertp is not None: Nlev = len(vertp) dvertp = vertp[0:Nlev-1] - vertp[1:Nlev] ax4 = ax3.twinx() l3 = ax3.plot(range(Nlev), vertp, 'r-o', label='pressure') l4 = ax4.plot(range(1,Nlev), dvertp, 'b-o', label='dpressure') ax3.set_xlim(0, Nlev) ax3.set_xlabel('level (\#)') ax3.set_ylabel('pressure (Pa)', color='r') ax3.set_ylim(np.min(vertp), np.max(vertp)) ax3.grid() ax4.set_ylabel('difference between levels (Pa)', color='b') if plog: ax3.set_yscale('log') if dplog: ax4.set_yscale('log') if vertp is not None: fig.suptitle(gtit) ax3.set_title('pressure') else: plt.title(gtit) output_kind(kfig, figname, True) return def plot_subbasin(subname, lons, lats, rsf, rLlf, rof, rcolors, mapv, basinwidth, \ drawsubid, graphtit, kfig, lloc, figname): """ Function to plot rivers from 'river_desc.nc' ORCDHIEE subname= name of the subbasin lons= values for the x-axis lats= values for the y-axis rsf= list of the name of the sub-flows to plot rLlf= dictionary with the lat,lon of the subflows rof= dictionary with the outflows of the subflows rcolors= base of colors of the lines for the subflows (first level) mapv= 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 basinwidth = with of the lines of the basin drawsubid= wehther sub-flow ids should be plot or not graphtit= title of the graph ('|', for spaces) lloc= location of the legend (0, automatic) 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' kfig= kind of figure figname= name of the figure """ fname = 'plot_subbasin' dx=lons.shape[1] dy=lats.shape[0] lengthtrac = np.float(lons[0,1] - lons[0,0]) lengthtrac2 = lengthtrac*np.sqrt(2.) ddx = np.abs(lons[dy/2+1,dx/2] - lons[dy/2,dx/2]) ddy = np.abs(lats[dy/2,dx/2+1] - lats[dy/2,dx/2]) fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 10} subhr = gen.chainSnumHierarchy(rsf) # Getting colors Nsubs = len(rsf) maxNgrids = -1 collev = {} if rcolors is None: # 1st level colors from colorsauto for isub in range(subhr.xlevs[1]): collev[gen.num_chainSnum(isub+1)] = gen.colhex_dec(colorsauto[isub]) else: if len(rcolors) != subhr.Llev1: print errormsg print ' ' + fname + ': number of first subbasins:', subhr.Llev1, \ 'and number of provided colors:', len(rcolors), 'differ !!' print ' provided colors:', rcolors quit(-1) for isub in range(subhr.xlevs[1]): collev[gen.num_chainSnum(isub+1)] = rcolors[isub] # All the colors should finish white... endcol = [1., 1., 1.] # endcol = [0.9, 0.9, 0.9] colorsub = {} xlabpos = [] ylabpos = [] xtrack = [] ytrack = [] colortrack = [] outflows = [] for isub in rsf: levs = gen.num_split(isub) Nlevs = len(levs) if Nlevs == 1: begcol = [0., 0., 0.] else: begcol = collev[levs[1]] for ilev in range(2,Nlevs): # We do not want a last white value... if Nlevs < subhr.Nlevs: maxNcol = subhr.xlevs[ilev]+1 else: maxNcol = subhr.xlevs[ilev]+1 colors = gen.color_deg([begcol, endcol], maxNcol, 'std') begcol = colors[1:4,gen.chainSnum_num(levs[ilev])-1] colorsub[isub] = tuple(begcol) # Avoiding too clear colors cols = colorsub[isub] if (cols[0] >= 0.95)*(cols[1] >= 0.95)*(cols[2] >= 0.95): print warnmsg print ' ' + fname + ': avoiding to clear color:', colorsub[isub], '!!' maxcol = np.max(cols) imaxcol = gen.index_vec(cols, maxcol) clear = np.ones((3), dtype=np.float)*0.9 clear[imaxcol] = 0.95 print ' keeping it clear:', clear colorsub[isub] = tuple(clear) Llfs = rLlf[isub] outfs = rof[isub] for igrid in range(len(outfs)): Llf = Llfs[igrid] ibeg = Llf[1] jbeg = Llf[0] outflows.append(outfs[igrid]) angle = (outfs[igrid] - 1)*np.pi/4 if gen.searchInlist([2,4,6,8], outfs[igrid]): iend = ibeg + lengthtrac2*np.sin(angle) jend = jbeg + lengthtrac2*np.cos(angle) elif gen.searchInlist([1,3,5,7], outfs[igrid]): iend = ibeg + lengthtrac*np.sin(angle) jend = jbeg + lengthtrac*np.cos(angle) else: ibeg = xvals[j,i] jbeg = yvals[j,i] iend = None jend = None xtrack.append(ibeg) xtrack.append(iend) xtrack.append(None) ytrack.append(jbeg) ytrack.append(jend) ytrack.append(None) colortrack.append(colorsub[isub]) xlabpos.append(ibeg) ylabpos.append(jbeg) plt.rc('text', usetex=True) nlon = np.min(xlabpos) xlon = np.max(xlabpos) nlat = np.min(ylabpos) xlat = np.max(ylabpos) dlon = xlon - nlon dlat = xlat - nlat # Making bigger the area to map # nlon = nlon-dlon-lengthtrac2 # xlon = xlon+dlon+lengthtrac2 # nlat = nlat-dlat-lengthtrac2 # xlat = xlat+dlat+lengthtrac2 nlon = nlon-2.*lengthtrac*dlon/np.abs(dlon) xlon = xlon+2.*lengthtrac*dlon/np.abs(dlon) nlat = nlat-2.*lengthtrac*dlat/np.abs(dlat) xlat = xlat+2.*lengthtrac*dlat/np.abs(dlat) plt.xlim(nlon,xlon) plt.ylim(nlat,xlat) if not mapv is None: lon00 = np.where(lons[:] < 0., 360. + lons[:], lons[:]) lat00 = lats[:] map_proj=mapv.split(',')[0] map_res=mapv.split(',')[1] lon2 = (nlon + xlon)/2. lat2 = (nlat + xlat)/2. print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \ xlon, ',', xlat if map_proj == 'cyl': m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \ urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) elif map_proj == 'lcc': m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \ llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res) else: print errormsg print ' ' + fname + ": projection '" + map_proj + "' not ready!!" print ' projections available: cyl, lcc' quit(-1) m.drawcoastlines() meridians = gen.pretty_int(nlon,xlon,5) m.drawmeridians(meridians,labels=[True,False,False,True],color="black") parallels = gen.pretty_int(nlat,xlat,5) m.drawparallels(parallels,labels=[False,True,True,False],color="black") plt.xlabel('W-E') plt.ylabel('S-N') j = 0 for i in range(len(xtrack)/3): plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i], linewidth=basinwidth) j = j + 3 # Sea-flow for i in range(len(outflows)): if outflows[i] == 97: plt.plot(xlabpos[i], ylabpos[i], 'x', color=colorsub[i]) elif outflows[i] == 98: plt.plot(xlabpos[i], ylabpos[i], '*', color=colorsub[i]) elif outflows[i] == 99: plt.plot(xlabpos[i], ylabpos[i], 'h', color=colorsub[i]) if drawsubid: for i in range(len(xlabpos)): plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, rsf[i], \ color=colorsub[rsf[i]], fontdict=fontcharac) plt.title(graphtit) output_kind(kfig, figname, True) return def plot_2lines(valsA, valsB, valsaxisA, valsaxisB, rangeA, rangeB, rangeaxis, \ axisvals, varns, varus, cols, wdths, styls, szmks, marks, gtitle, axisn, gloc, \ sloc, fign, figk, ifclose): """ Function to plot two lines in different axes (x/x2 or y/y2) valsA= values to be plotted on axis x or y valsB= values to be plotted on axis x2 or y2 valsaxisA= values at the common axis for valsA valsaxisB= values at the common axis for valsB rangeA= range of values for valsA rangeB= range of values for valsB rangeaxis= range of values for the common axis axisvals= which is the axis to plot the values ('x' or 'y') varns= names of the variables in the plot varus= units of the variables cols= list with color names of the lines for the variables wdths= list with widths of the lines for the variables styls= list with the styles of the lines szmks= list with the size of the markers of the lines marks= list with the markers of the lines gtitle= title of the figure axisn= label in the figure of the common axis gloc= location of the legend gsiz= fontsize of the legend fign= name of the figure figk= kind of figure ifclose: Whether figure should be finished or not """ fname = 'plot_2lines' if axisvals == 'x': titX = axisn titX2 = axisn titY = varns[0] + ' (' + varus[0] + ')' titY2 = varns[1] + ' (' + varus[1] + ')' xmin = rangeaxis[0] xmax = rangeaxis[1] xmin2 = rangeaxis[0] xmax2 = rangeaxis[1] ymin = rangeA[0] ymax = rangeA[1] ymin2 = rangeB[0] ymax2 = rangeB[1] else: titX = varns[0] + ' (' + varus[0] + ')' titX2 = varns[1] + ' (' + varus[1] + ')' titY = axisn titY2 = axisn xmin = rangeA[0] xmax = rangeA[1] xmin2 = rangeB[0] xmax2 = rangeB[1] ymin = rangeaxis[0] ymax = rangeaxis[1] ymin2 = rangeaxis[0] ymax2 = rangeaxis[1] # From: http://matplotlib.org/examples/api/two_scales.html fig, ax1 = plt.subplots() plt.rc('text', usetex=True) # Plotting lines il=0 if axisvals == 'x': lA, = ax1.plot(valsaxisA, valsA, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax1.set_xlabel(titX, color='black') ax1.set_ylabel(titY, color=cols[il]) ax1.set_xlim(xmin,xmax) ax1.set_ylim(ymin,ymax) ax2 = ax1.twinx() il = il + 1 lB, = ax2.plot(valsaxisB, valsB, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax2.set_ylabel(titY2, color=cols[il]) ax2.set_ylim(ymin2,ymax2) ax2.set_xlim(xmin,xmax) titleloc = (0.5,1.) else: lA, = ax1.plot(valsA, valsaxisA, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax1.set_xlabel(titX, color=cols[il]) ax1.set_ylabel(titY, color='black') ax1.set_xlim(xmin,xmax) ax1.set_ylim(ymin,ymax) ax2 = ax1.twiny() il = il + 1 lB, = ax2.plot(valsB, valsaxisB, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax2.set_xlabel(titX2, color=cols[il]) ax2.set_xlim(xmin2,xmax2) ax2.set_ylim(ymin,ymax) titleloc = (0.5,1.075) graphtit = gen.latex_text(gtitle) plt.title(graphtit,position=titleloc) # plt.legend([lA, lB], loc=gloc, prop={'size':gsize}) output_kind(figk, fign, ifclose) return def plot_2lines_time(valsA, valsB, valsaxisA, valsaxisB, rangeA, rangeB, valstaxis, \ labelstaxis, taxis, varns, varus, cols, wdths, styls, szmks, marks, gtitle, axisn, \ gloc, gsize, fign, figk, ifclose): """ Function to plot two time-lines in different axes (x/x2 or y/y2) valsA= values to be plotted on axis x or y valsB= values to be plotted on axis x2 or y2 valsaxisA= values at the common axis for valsA valsaxisB= values at the common axis for valsB rangeA= range of values for valsA rangeB= range of values for valsB valstaxis= values for the tick in the time-axis labelstaxis= labels for the tick in the time-axis taxis= which is the axis for the time values ('x' or 'y') varns= names of the variables in the plot varus= units of the variables cols= list with color names of the lines for the variables wdths= list with widths of the lines for the variables styls= list with the styles of the lines szmks= list with the size of the markers of the lines marks= list with the markers of the lines gtitle= title of the figure axisn= label in the figure of the time axis gloc= location of the legend gsiz= fontsize of the legend fign= name of the figure figk= kind of figure ifclose: Whether figure should be finished or not """ fname = 'plot_2lines_time' if taxis == 'x': titX = axisn titX2 = axisn titY = varns[0] + ' (' + varus[0] + ')' titY2 = varns[1] + ' (' + varus[1] + ')' xmin = np.min(valstaxis) xmax = np.max(valstaxis) ymin = rangeA[0] ymax = rangeA[1] ymin2 = rangeB[0] ymax2 = rangeB[1] else: titX = varns[0] + ' (' + varus[0] + ')' titX2 = varns[1] + ' (' + varus[1] + ')' titY = axisn titY2 = axisn xmin = rangeA[0] xmax = rangeA[1] xmin2 = rangeB[0] xmax2 = rangeB[1] ymin = np.min(valstaxis) ymax = np.max(valstaxis) # From: http://matplotlib.org/examples/api/two_scales.html fig, ax1 = plt.subplots() plt.rc('text', usetex=True) # Plotting lines il=0 if taxis == 'x': lA, = ax1.plot(valsaxisA, valsA, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax1.set_xlabel(titX, color='black') ax1.set_ylabel(titY, color=cols[il]) ax1.set_xlim(xmin,xmax) ax1.set_ylim(ymin,ymax) ax1.set_xticks(valstaxis) ax1.set_xticklabels(labelstaxis) ax2 = ax1.twinx() il = il + 1 lB, = ax2.plot(valsaxisB, valsB, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax2.set_ylabel(titY2, color=cols[il]) ax2.set_ylim(ymin2,ymax2) ax2.set_xlim(xmin,xmax) ax2.set_xticks(valstaxis) ax2.set_xticklabels(labelstaxis) titleloc = (0.5,1.) else: lA, = ax1.plot(valsA, valsaxisA, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax1.set_xlabel(titX, color=cols[il]) ax1.set_ylabel(titY, color='black') ax1.set_xlim(xmin,xmax) ax1.set_ylim(ymin,ymax) ax1.set_yticks(valstaxis) ax1.set_yticklabels(labelstaxis) ax2 = ax1.twiny() il = il + 1 lB, = ax2.plot(valsB, valsaxisB, styls[il], marker=marks[il], \ linewidth=wdths[il], markersize=szmks[il], label=varns[il], color=cols[il]) ax2.set_xlabel(titX2, color=cols[il]) ax2.set_xlim(xmin2,xmax2) ax2.set_ylim(ymin,ymax) ax2.set_yticks(valstaxis) ax2.set_yticklabels(labelstaxis) titleloc = (0.5,1.075) graphtit = gen.latex_text(gtitle) plt.title(graphtit,position=titleloc) # plt.legend([lA, lB], loc=gloc, prop={'size':gsize}) output_kind(figk, fign, ifclose) return