From 0e315ee7a0d43deb5854ad06b3b17c848369b3b1 Mon Sep 17 00:00:00 2001 From: Francesc Alted Date: Tue, 13 Sep 2016 14:35:15 +0200 Subject: [PATCH 1/4] First attempt at making part of reflexible Python 2/3 compliant --- reflexible/conv2netcdf4/flexpart_read.py | 24 +- .../conv2netcdf4/fortflex/build_FortFlex.sh | 2 +- reflexible/conv2netcdf4/grid_read.py | 205 +- reflexible/data_structures.py | 4 +- reflexible/fprof/prof.py | 4 +- reflexible/legacy/mapping.py | 3239 +++++++++-------- reflexible/plotting.py | 8 +- reflexible/tests/test_plotting.py | 14 +- 8 files changed, 1768 insertions(+), 1732 deletions(-) diff --git a/reflexible/conv2netcdf4/flexpart_read.py b/reflexible/conv2netcdf4/flexpart_read.py index 0d2d81e..0a27dc2 100644 --- a/reflexible/conv2netcdf4/flexpart_read.py +++ b/reflexible/conv2netcdf4/flexpart_read.py @@ -1,5 +1,7 @@ #### Functions for reading FLEXPART output ##### +from __future__ import print_function + import os import re import datetime @@ -162,7 +164,7 @@ def read_command_old(path, headerrows): ================ ================================================================= """ - lines = file(path, 'r').readlines() + lines = open(path, 'r').readlines() command_vals = [i.strip() for i in lines[headerrows:]] # clean line ends COMMAND_KEYS = ( 'SIM_DIR', @@ -251,7 +253,7 @@ def getfile_lines(infile): reverts to beginning of file.""" if isinstance(infile, str): - return file(infile, 'r').readlines() + return open(infile, 'r').readlines() else: infile.seek(0) return infile.readlines() @@ -375,12 +377,12 @@ def read_trajectories(H, trajfile='trajectories.txt', if isinstance(H, str): try: - alltraj = file(H, 'r').readlines() + alltraj = open(H, 'r').readlines() except: raise IOError('Could not open file: %s' % H) else: path = H.path - alltraj = file(os.path.join(path, trajfile), 'r').readlines() + alltraj = open(os.path.join(path, trajfile), 'r').readlines() try: ibdate, ibtime, model, version = alltraj[0].strip().split()[:4] @@ -505,7 +507,7 @@ def curtain_agltoasl(H, curtain_agl, coords, below_gl=0.0): xp = H.outheight - H.outheight[0] casl = np.zeros((len(H.asl_axis), len(coords))) - for i in xrange(len(coords)): + for i in range(len(coords)): casl[:, i] = np.interp(H.asl_axis, xp + gl[i], curtain_agl[:, i], left=below_gl) return casl @@ -543,7 +545,7 @@ def read_agespectrum(filename, part=False, ndays=20): """ - f = file(filename, 'r').readlines() + f = open(filename, 'r').readlines() line1 = f[0].strip().split() if part: @@ -666,7 +668,7 @@ def gridarea(H): ylata = outlat0 + (float(iy) + 0.5) * dyout ylatp = ylata + 0.5 * dyout ylatm = ylata - 0.5 * dyout - if (ylatm < 0 and ylatp > 0): + if ylatm < 0 and ylatp > 0: hzone = dyout * r_earth * pih else: # cosfact = cosfunc(ylata) @@ -792,10 +794,10 @@ def read_header(pathname, **kwargs): OPS.update(kwargs) if OPS.verbose: - print "Reading Header with:\n" + print("Reading Header with:\n") for o in OPS: - print "%s ==> %s" % (o, OPS[o]) + print("%s ==> %s" % (o, OPS[o])) # Define utility functions for reading binary file skip = lambda n = 8: bf.seek(n, 1) @@ -831,7 +833,7 @@ def read_header(pathname, **kwargs): raise IOError("No DATEFILE: {0}".format(datefile)) else: try: - fd = file(datefile, 'r').readlines() + fd = open(datefile, 'r').readlines() except: raise IOError("Could not read datefile: {0}".format(datefile)) @@ -1180,7 +1182,7 @@ def _readV6(bf, h): if bf: # bf.read('i') - print h['numpoint'] + print(h['numpoint']) for i in range(h['numpoint']): # r2=getbin('i') i1 = getbin('i') diff --git a/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh b/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh index a06f106..9f3ce4e 100755 --- a/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh +++ b/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh @@ -3,5 +3,5 @@ rm -f FortFlex.pyf f2py -m FortFlex -h FortFlex.pyf FortFlex.f f2py -c --fcompiler=gfortran FortFlex.pyf FortFlex.f -mv FortFlex.so .. +mv *.so .. diff --git a/reflexible/conv2netcdf4/grid_read.py b/reflexible/conv2netcdf4/grid_read.py index aff18b4..5a0e8c8 100644 --- a/reflexible/conv2netcdf4/grid_read.py +++ b/reflexible/conv2netcdf4/grid_read.py @@ -1,5 +1,7 @@ ########### Grid Reading Routines ############### +from __future__ import print_function + import itertools import datetime import os @@ -11,7 +13,6 @@ from .helpers import _shout - def _readgrid_noFF(H, **kwargs): """ accepts a header dictionary as input, returns dictionary of Grid values %=========================================== @@ -43,10 +44,10 @@ def _readgrid_noFF(H, **kwargs): # try: # from FortFlex import readgrid, sumgrid # useFortFlex = 1 - # print 'using FortFlex' + # print('using FortFlex') # except: # useFortFlex = 0 - # print 'Cannot find FortFlex.so' + # print('Cannot find FortFlex.so') useFortFlex = 0 if 'date' in kwargs.keys(): @@ -111,21 +112,21 @@ def _readgrid_noFF(H, **kwargs): 'grid_time_nest_', 'footprint_total_nest'] # local functions for reading binary data and creating grid dump -# skip = lambda n=8 : f.read(n) -# getbin = lambda fmt,n=1 : struct.unpack(fmt*n,f.read(struct.calcsize(fmt)))[0] + # skip = lambda n=8 : f.read(n) + # getbin = lambda fmt,n=1 : struct.unpack(fmt*n,f.read(struct.calcsize(fmt)))[0] # Utility functions - skip = lambda n = 8: f2.seek(n, 1) - getbin = lambda dtype, n = 1: f2.read(dtype, (n,)) + skip = lambda n=8: f2.seek(n, 1) + getbin = lambda dtype, n=1: f2.read(dtype, (n,)) def getdump(n, fmt='f'): """ function to get the dump values for the sparse format """ - # print n + # print(n) skip() # Dfmt=[fmt]*n -# a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt] + # a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt] a = f2.read(fmt, n) -# dumplist=[a[j][0] for j in range(len(a))] + # dumplist=[a[j][0] for j in range(len(a))] # dumplist=[a[j] for j in range(len(a))] return a # dumplist @@ -148,10 +149,10 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): kz = n / (numxgrid * numygrid) jy = (n - kz * numxgrid * numygrid) / numxgrid ix = n - numxgrid * numygrid * kz - numxgrid * jy - # print "n ==> ix,jy,kz,k,nage" - # print "%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage) - # print grd.shape - # print grd[0,0,0,0,0] + # print("n ==> ix,jy,kz,k,nage") + # print("%s ==> %s,%s,%s,%s,%s") % (n,ix,jy,kz,k,nage) + # print(grd.shape) + # print(grd[0,0,0,0,0]) grd[ix, jy, kz - 1, k, nage] = abs(dmp_r[ir]) return grd # flipud(grd.transpose()) @@ -170,12 +171,12 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): dtype=np.float32) zplot = np.empty((numxgrid, numygrid, numzgrid, numpoint, nageclass)) - #-------------------------------------------------- + # -------------------------------------------------- # Loop over all times, given in field H['dates'] - #-------------------------------------------------- + # -------------------------------------------------- for ks in range(nspec_ret, nspec_ret + 1): - # print 'processing: ' + str(H['dates'][date_i]) + ' ' + str(ks) - # print 'processing: ' + str(date) + ' ' + str(ks) + # print('processing: ' + str(H['dates'][date_i]) + ' ' + str(ks)) + # print('processing: ' + str(date) + ' ' + str(ks)) specs = '_' + str(ks).zfill(3) fpmax = -999. if date is None and time_ret is None: @@ -191,7 +192,7 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): else: get_dates = available_dates - print 'getting grid for: ', get_dates + print('getting grid for: ', get_dates) for date_i in range(len(get_dates)): if unit != 4: @@ -203,15 +204,15 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): filename = os.path.join(H['pathname'], prefix[(unit) + (nested * 4)]) - # print 'reading: ' + filename + # print('reading: ' + filename) if os.path.exists(filename): - # print nspec,numxgrid,numygrid,numzgrid,nageclass,scaledepo,scaleconc,decaycons - # print ks - # print date + # print(nspec,numxgrid,numygrid,numzgrid,nageclass,scaledepo,scaleconc,decaycons) + # print(ks) + # print(date) if useFortFlex == 1: - # print 'Using FortFLEX' + # print('Using FortFLEX') ########################################################### # FORTRAN WRAPPER CODE - only works on linux # # @@ -232,8 +233,8 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): scaleconc, decayconstant) # contribution[:,:,:] = concgrid[:,:,:,:,0] - print np.min(concgrid) - print np.max(concgrid) + print(np.min(concgrid)) + print(np.max(concgrid)) # altitude = 50000 zplot = sumgrid(zplot, concgrid, @@ -247,12 +248,12 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): # nspec,nageclass),np.float) datagrid = np.zeros( (numxgrid, numygrid, numzgrid, 1, 1), np.float) - # f = file(filename, 'rb') - # print filename + # f = open(filename, 'rb') + # print(filename) f2 = reflexible.conv2netcdf4.BinaryFile(filename, order='fortran') skip(4) G['itime'] = getbin('i') - print H['available_dates'][date_i] + print(H['available_dates'][date_i]) # Read Wet Depostion skip() @@ -278,36 +279,37 @@ def dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage): skip() cnt_r = getbin('i') dmp_r = getdump(cnt_r) - # print dmp_i, cnt_r, dmp_r, datagrid, ks-1, nage + # print(dmp_i, cnt_r, dmp_r, datagrid, ks-1, nage) concgrid = dumpgrid(dmp_i, cnt_r, dmp_r, datagrid, ks - 1, nage) # G[H['ibtime']].append(concgrid) G[H['ibtime']] = concgrid f2.close() fail = 0 else: - print "\n\n INPUT ERROR: Could not find file: %s" % filename + print("\n\n INPUT ERROR: Could not find file: %s" % filename) raise IOError('No file: %s' % filename) if useFortFlex == 1: G = zplot return G + readgridBF = _readgrid_noFF def _readgridBF(H, filename): """ Read grid using BinaryFile class""" # Utility functions - skip = lambda n = 8: f2.seek(n, 1) - getbin = lambda dtype, n = 1: f2.read(dtype, (n,)) + skip = lambda n=8: f2.seek(n, 1) + getbin = lambda dtype, n=1: f2.read(dtype, (n,)) def getdump(n, fmt='f'): """ function to get the dump values for the sparse format """ skip() # Dfmt=[fmt]*n -# a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt] + # a=[struct.unpack(ft,f.read(struct.calcsize(ft))) for ft in Dfmt] a = f2.read(fmt, n) -# dumplist=[a[j][0] for j in range(len(a))] + # dumplist=[a[j][0] for j in range(len(a))] # dumplist=[a[j] for j in range(len(a))] return a # dumplist @@ -339,11 +341,11 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): ix = n - H.numxgrid * H.numygrid * kz - H.numxgrid * jy grd[ix, jy, kz - 1, k, nage] = abs(dmp_r[ir]) -# -# print "n ==> ix,jy,kz,k,nage" -# print "%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage) -# print grd.shape -# print grd[0,0,0,0,0] + # + # print("n ==> ix,jy,kz,k,nage") + # print("%s ==> %s,%s,%s,%s,%s" % (n,ix,jy,kz,k,nage)) + # print(grd.shape) + # print(grd[0,0,0,0,0]) else: if dmp_r[ir] * fact > 0: @@ -362,15 +364,17 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): try: # Using Pyximport for compiling on-the flight. # See: http://docs.cython.org/src/userguide/source_files_and_compilation.html#pyximport - import pyximport; pyximport.install() + import pyximport; + pyximport.install() from pflexcy import dumpdatagrid, dumpdepogrid - print 'using pflexcy' + print('using pflexcy') except: - print """WARNING: Using PURE Python to readgrid, execution will be slow. - Try compiling the FortFlex module or the pflexcy module - for your machine. For more information see the - reflexible/f2py_build directory or use cython with pflexcy.pyx - """ + print( + """WARNING: Using PURE Python to readgrid, execution will be slow. + Try compiling the FortFlex module or the pflexcy module + for your machine. For more information see the + reflexible/f2py_build directory or use cython with pflexcy.pyx + """) dumpdatagrid = _dumpgrid dumpdepogrid = _dumpgrid @@ -381,8 +385,8 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): drygrid = np.zeros((H.numxgrid, H.numygrid, H.numpointspec, 1), np.float) datagrid = np.zeros( (H.numxgrid, H.numygrid, H.numzgrid, H.numpointspec, nage), np.float) - # f = file(filename,'rb') - # print filename + # f = open(filename,'rb') + # print(filename) f2 = reflexible.conv2netcdf4.BinaryFile(filename, order='fortran') # read data: skip(4) @@ -399,7 +403,7 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): cnt_r = getbin('i') dmp_r = getdump(cnt_r) if dmp_r.any(): - # print dmp_r, dmp_i + # print(dmp_r, dmp_i) wetgrid = dumpdepogrid( dmp_i, cnt_r, dmp_r, wetgrid, ks, na, H.numxgrid, H.numygrid) @@ -409,7 +413,7 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): cnt_r = getbin('i') dmp_r = getdump(cnt_r) if dmp_r.any(): - # print dmp_r, dmp_i + # print(dmp_r, dmp_i) drygrid = dumpdepogrid( dmp_i, cnt_r, dmp_r, drygrid, ks, na, H.numxgrid, H.numygrid) @@ -421,11 +425,11 @@ def _dumpgrid(dmp_i, cnt_r, dmp_r, grd, k, nage, nx, ny): skip() cnt_r = getbin('i') dmp_r = getdump(cnt_r) - # print len(dmp_r),len(dmp_i) - # print cnt_r,cnt_i - # print dmp_i - # print type(dmp_i),type(cnt_r),type(dmp_r),type(datagrid),type(ks),type(na) - # print type(H.numxgrid),type(H.numygrid) + # print(len(dmp_r),len(dmp_i)) + # print(cnt_r,cnt_i) + # print(dmp_i) + # print(type(dmp_i),type(cnt_r),type(dmp_r),type(datagrid),type(ks),type(na)) + # print(type(H.numxgrid),type(H.numygrid)) datagrid = dumpdatagrid( dmp_i, cnt_r, dmp_r, datagrid, ks, na, H.numxgrid, H.numygrid) @@ -450,7 +454,7 @@ def _read_headerFF(pathname, h=None, try: from .FortFlex import readheader except: - print "Error with FortFlex.readheader, use read_header" + print("Error with FortFlex.readheader, use read_header") headervars = ['numxgrid', 'numygrid', 'numzgrid', 'outlon0', 'outlat0', 'compoint', 'dxout', 'dyout', 'outheight', 'ibdate', @@ -464,14 +468,16 @@ def _read_headerFF(pathname, h=None, h = reflexible.conv2netcdf4.Structure() if verbose: - print """Reading Header with: - maxpoint : %s - maxspec : %s - maxageclass : %s - nxmax : %s - nymax : %s - nzmax : %s - """ % (maxpoint, maxspec, maxageclass, nxmax, nymax, nzmax) + print( + """Reading Header with: + maxpoint : %s + maxspec : %s + maxageclass : %s + nxmax : %s + nymax : %s + nzmax : %s + """ % (maxpoint, maxspec, maxageclass, nxmax, nymax, nzmax) + ) (numxgrid, numygrid, numzgrid, outlon0, outlat0, dxout, dyout, outheight, ibdate, ibtime, loutstep, nspec, nageclass, lage, ireleasestart, @@ -621,7 +627,7 @@ def readgridV8(H, **kwargs): if time_ret[0] < 0: # if forward == False: - # get all dates for calculating footprint. + # get all dates for calculating footprint. time_ret = np.arange(len(H.available_dates)) # else: # raise ValueError("Must enter a positive time_ret for forward runs") @@ -650,7 +656,7 @@ def readgridV8(H, **kwargs): # assign grid dates for indexing fd fd.grid_dates = get_dates[:] - print 'getting grid for: ', get_dates + print('getting grid for: ', get_dates) # Some pre-definitions fail = 0 # set filename prefix @@ -676,7 +682,7 @@ def readgridV8(H, **kwargs): from nilu.pflexpart.FortFlex import readgrid, sumgrid useFortFlex = True print('Using old nilu.pflexpart.FortFlex') - # print 'using nilu.pflexpart FortFlex' + # print('using nilu.pflexpart FortFlex') except: useFortFlex = False print('Cannot load FortFlex, reverting to BinaryFile.') @@ -686,20 +692,20 @@ def readgridV8(H, **kwargs): print('Using BinaryFile') # reserve output fields - print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint + print(H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint) # ------------------------------------------------- # add the requests to the fd object to be returned OPS.unit = unit - #-------------------------------------------------- + # -------------------------------------------------- # Loop over all times, given in field H['available_dates'] - #-------------------------------------------------- + # -------------------------------------------------- for date_i in range(len(get_dates)): datestring = get_dates[date_i] - print datestring + print(datestring) for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1):A FLEXDATA[(s, datestring)] = fdc = reflexible.conv2netcdf4.FDC() spec_fid = '_' + str(s + 1).zfill(3) @@ -712,7 +718,7 @@ def readgridV8(H, **kwargs): else: # grid total footprint - print "Total footprint" + print("Total footprint") filename = os.path.join( H['pathname'], prefix[(unit_i) + (H.nested * 5)] + spec_fid) @@ -720,15 +726,15 @@ def readgridV8(H, **kwargs): if os.path.exists(filename): H.filename = filename - # print 'reading: ' + filename + # print('reading: ' + filename) if OPS.verbose: - print 'with values:' + print('with values:') inputvars = ['filename', 'numxgrid', 'numygrid', 'zdims', 'numpoint', 'nageclass', 'scaledepo', 'scaleconc', 'decayconstant', 'numpointspec'] for v in inputvars: - print v, " ==> ", H[v] + print(v, " ==> ", H[v]) if OPS.BinaryFile: print("Reading {0} with BinaryFile".format(filename)) @@ -899,7 +905,7 @@ def readgridV6(H, **kwargs): # assign grid dates for indexing fd fd.grid_dates = get_dates[:] - print 'getting grid for: ', get_dates + print('getting grid for: ', get_dates) # Some pre-definitions fail = 0 # set filename prefix @@ -918,7 +924,7 @@ def readgridV6(H, **kwargs): from .FortFlex import sumgrid useFortFlex = True OPS.useFortFlex = useFortFlex - print 'using FortFlex VERSION 6' + print('using FortFlex VERSION 6') except: # get the original module (no memory allocation) useFortFlex = False @@ -927,16 +933,16 @@ def readgridV6(H, **kwargs): readgrid = _readgridBF OPS.BinaryFile = True - # cheat: use key2var function to get values from header dict, H - # need to change this to using H.xxxxx - # headervars = ['nested','nspec','numxgrid','numygrid','numzgrid','nageclass',\ - # 'dates','pathname','decayconstant','numpoint','numpointspec', - # 'area','Heightnn','lage'] - # for k in headervars: - # key2var(H,k) + # cheat: use key2var function to get values from header dict, H + # need to change this to using H.xxxxx + # headervars = ['nested','nspec','numxgrid','numygrid','numzgrid','nageclass',\ + # 'dates','pathname','decayconstant','numpoint','numpointspec', + # 'area','Heightnn','lage'] + # for k in headervars: + # key2var(H,k) - # reserve output fields - print H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint + # reserve output fields + print(H.numxgrid, H.numygrid, H.numzgrid, OPS.nspec_ret, OPS.pspec_ret, OPS.age_ret, len(get_dates), H.numpoint) # ------------------------------------------------- @@ -951,13 +957,13 @@ def readgridV6(H, **kwargs): # grid=np.empty((numxgrid,numygrid,numzgrid,nspec_ret,pspec_ret,age_ret,len(get_dates))) # zplot=np.empty((numxgrid,numygrid,numzgrid,numpointspec)) - #-------------------------------------------------- + # -------------------------------------------------- # Loop over all times, given in field H['available_dates'] - #-------------------------------------------------- + # -------------------------------------------------- for date_i in range(len(get_dates)): datestring = get_dates[date_i] - print datestring + print(datestring) for s in nspec_ret: # range(OPS.nspec_ret,OPS.nspec_ret+1): FLEXDATA[(s, datestring)] = fdc = reflexible.conv2netcdf4.FDC() # spec_fid = '_'+str(s+1).zfill(3) @@ -969,25 +975,25 @@ def readgridV6(H, **kwargs): else: # grid total footprint - print "Total footprint" + print("Total footprint") filename = os.path.join(H['pathname'], prefix[(unit_i) + (H.nested * 5)]) H.zdims = 1 if os.path.exists(filename): H.filename = filename - # print 'reading: ' + filename + # print('reading: ' + filename) if OPS.verbose: - print 'with values:' + print('with values:') inputvars = ['filename', 'numxgrid', 'numygrid', 'zdims', 'numpoint', 'nageclass', 'scaledepo', 'scaleconc', 'decayconstant', 'numpointspec'] for v in inputvars: - print v, " ==> ", H[v] + print(v, " ==> ", H[v]) if OPS.BinaryFile: - # print 'Using BinaryFile' + # print('Using BinaryFile') gridT, wetgrid, drygrid, itime = _readgridBF(H, filename) else: gridT, wetgrid, drygrid, itime = readgrid( @@ -1044,12 +1050,11 @@ def readgridV6(H, **kwargs): def monthly_footprints(H): - footprints = np.zeros((H.ny, H.nx, len(H.C.keys()))) for i, key in enumerate(H.C): footprints[:, :, i] = H.C[key].slabs[0] - #footprints = np.average(footprints, axis=2) + # footprints = np.average(footprints, axis=2) return footprints @@ -1090,7 +1095,7 @@ def fill_grids(H, nspec=0, FD=None): """ -# assert H.direction == 'backward', "fill_backward is only valid for backward runs" + # assert H.direction == 'backward', "fill_backward is only valid for backward runs" # make sure npsec is iterable if isinstance(nspec, int): species = [nspec] @@ -1117,7 +1122,7 @@ def fill_grids(H, nspec=0, FD=None): c.spec_i = s # read data grids and attribute/sum sensitivity - print species + print(species) for s in species: for d in FD.grid_dates: @@ -1216,7 +1221,7 @@ def get_slabs(H, G, index=None, normAreaHeight=True, scale=1.0): area = H['area'] Slabs = reflexible.conv2netcdf4.Structure() grid_shape = G.shape - nageclass = 0 # XXX assume nageclass equal to 0 + nageclass = 0 # XXX assume nageclass equal to 0 if len(grid_shape) is 5: if index is None: if H.direction is 'forward': diff --git a/reflexible/data_structures.py b/reflexible/data_structures.py index 2ee9428..b28b9de 100755 --- a/reflexible/data_structures.py +++ b/reflexible/data_structures.py @@ -724,9 +724,9 @@ def __init__(self, **options): ## below here, not actual COMMAND input 'HEADER': """**********************************************\n\n\n Input file for FLEXPART\n\n*********************************************\n\n""", 'FLEXPART_VER': [10, '''FLEXPART VERSION Used to define format of COM MAND File'''], - 'SIM_START': [dt.datetime(2000, 01, 01, 00, 00, 00), + 'SIM_START': [dt.datetime(2000, 1, 1, 00, 00, 00), '''Beginning date and time of simulation. Must be given in format YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR, MI is MINUTE and SS is SECOND. Current version utilizes UTC.'''], - 'SIM_END': [dt.datetime(2000, 02, 01, 00, 00, 00), + 'SIM_END': [dt.datetime(2000, 2, 1, 00, 00, 00), '''Ending date and time of simulation. Same format as 2'''], 'AGECLASSES': [[86400 * 30], '''list of ageclasses (seconds) in the simulation'''], 'RELEASE_SECONDS': [86400, '''duration of the releases in seconds'''] diff --git a/reflexible/fprof/prof.py b/reflexible/fprof/prof.py index 23885cb..57d49e3 100644 --- a/reflexible/fprof/prof.py +++ b/reflexible/fprof/prof.py @@ -14,7 +14,7 @@ import time import cProfile import pstats -import StringIO +import io # Global variable for disabling all the profiles @@ -58,7 +58,7 @@ def cprof(explain='', nlines=20): pr.enable() yield pr.disable() - s = StringIO.StringIO() + s = io.StringIO() #sortby = ('time', 'calls') sortby = ('cumulative',) ps = pstats.Stats(pr, stream=s).sort_stats(*sortby) diff --git a/reflexible/legacy/mapping.py b/reflexible/legacy/mapping.py index 30c0244..be892a5 100644 --- a/reflexible/legacy/mapping.py +++ b/reflexible/legacy/mapping.py @@ -2,27 +2,29 @@ """ Matplotlib Basemap Tool Suite """ # -*- coding: utf-8 -*- -#John F Burkhart - 2010 -#Licence : this code is released under the matplotlib license +# John F Burkhart - 2010 +# Licence : this code is released under the matplotlib license + +from __future__ import print_function __author__ = "John F Burkhart " __version__ = "0.02" - import numpy as np import math import copy import datetime as dt import pdb -#from matplotlib import interactive, use -#use("Agg") -#interactive(False) +# from matplotlib import interactive, use +# use("Agg") +# interactive(False) import matplotlib as mpl from matplotlib import colors, cm from matplotlib.collections import LineCollection -#mpl.use('Agg') +# mpl.use('Agg') import matplotlib.pyplot as plt + try: from mpl_toolkits.basemap import Basemap, shiftgrid except: @@ -34,17 +36,18 @@ # !!NEED TO FIX!! # from matplotlib.ticker import NullFormatter + class Structure(dict): def __getattr__(self, attr): return self[attr] + def __setattr__(self, attr, value): self[attr] = value - def set_with_dict(self,D): + def set_with_dict(self, D): """ set attributes with a dict """ for k in D.keys(): - self.__setattr__(k,D[k]) - + self.__setattr__(k, D[k]) class KML_File: @@ -54,81 +57,86 @@ class KML_File: modifications, John Burkhart \n"\ - "\n"\ - "\n"\ - "\n") + "\n" \ + "\n" \ + "\n" \ + "\n") file.close() def close(self): - file = open(self.filepath,"a") + file = open(self.filepath, "a") file.write( - "\n"\ - "") + "\n" \ + "") file.close() def open_folder(self, name): - file = open(self.filepath,"a") + file = open(self.filepath, "a") file.write( - "\n"\ - " " + name + "\n") + "\n" \ + " " + name + "\n") file.close() def close_folder(self): - file = open(self.filepath,"a") + file = open(self.filepath, "a") file.write( - "\n") + "\n") file.close() - def add_placemarker(self, latitude, longitude, altitude = 0.0,\ - description = " ", name = " ", \ - range = 6000, tilt = 45, heading = 0): + def add_placemarker(self, latitude, longitude, altitude=0.0, \ + description=" ", name=" ", \ + range=6000, tilt=45, heading=0): "adds the point to a kml file" - file = open(self.filepath,"a") + file = open(self.filepath, "a") if not np.iterable(longitude) and not np.iterable(latitude): point = True - print 'no lines' + print + 'no lines' else: point = False - print 'lines!' + print + 'lines!' if point: file.write( - "\n"\ - " " + description + "\n"\ - " " + name + "\n"\ - " #normalPlaceMarker" + - " \n"\ - " " + str(longitude) + "\n"\ - " " + str(latitude) + "\n"\ - " " + str(range) + "\n"\ - " " + str(tilt) + "\n"\ - " " + str(heading) + "\n"\ - " \n"\ - " 0\n"\ - " \n"\ - " 1\n"\ - " relativeToGround\n"\ - " " + str(longitude) + "," + str(latitude) +", " + str(altitude) + "\n"\ - " \n"\ - " \n") + "\n" \ + " " + description + "\n" \ + " " + name + "\n" \ + " #normalPlaceMarker" + + " \n" \ + " " + str(longitude) + "\n" \ + " " + str(latitude) + "\n" \ + " " + str( + range) + "\n" \ + " " + str(tilt) + "\n" \ + " " + str(heading) + "\n" \ + " \n" \ + " 0\n" \ + " \n" \ + " 1\n" \ + " relativeToGround\n" \ + " " + str( + longitude) + "," + str(latitude) + ", " + str(altitude) + "\n" \ + " \n" \ + " \n") file.close() else: @@ -136,38 +144,38 @@ def add_placemarker(self, latitude, longitude, altitude = 0.0,\ file.write( - "\n"\ - " " + description + "\n"\ - " " + name + "\n"\ - " #normalPlaceMarker" + - " 1\n"\ - " \n"\ - " 0\n"\ - " 0\n"\ - " \n") - for i,lat in enumerate(latitude): - file.write("%s,%s,0 " % (longitude[i],lat)) + "\n" \ + " " + description + "\n" \ + " " + name + "\n" \ + " #normalPlaceMarker" + + " 1\n" \ + " \n" \ + " 0\n" \ + " 0\n" \ + " \n") + for i, lat in enumerate(latitude): + file.write("%s,%s,0 " % (longitude[i], lat)) file.write("\n") file.write( - " \n"\ - " \n"\ - " \n") + " \n" \ + " \n" \ + " \n") file.close() - -def _val2var(val,var): +def _val2var(val, var): """ a helper function to convert a value to a variable. Hopefully it is not still used.""" - if isinstance(val,str): - cmd = "global %s; %s = '%s';" % (var,var,val) - if isinstance(val,str) == False: - cmd = "global %s; %s = %s;" % (var,var,val) + if isinstance(val, str): + cmd = "global %s; %s = '%s';" % (var, var, val) + if not isinstance(val, str): + cmd = "global %s; %s = %s;" % (var, var, val) - exec(cmfd) + exec(cmd) -def map_dynaMap((x,y)): + +def map_dynaMap(x, y): """Attempts to define the basemap instance settings based on lon/lat input. @@ -182,30 +190,35 @@ def map_dynaMap((x,y)): This function is not reliable yet. """ - ymin=min(y); ymax=max(y); y_edge=abs(ymax-ymin)*0.75 - xmin=min(x); xmax=max(x); x_edge=abs(xmax-xmin)*0.75 - - llcrnrlat= ymin-y_edge - if llcrnrlat<-90.: llcrnrlat=-90. - llcrnrlon= xmin-x_edge - if llcrnrlon<-180.: llcrnrlon=-180. - urcrnrlat= ymax+y_edge - if urcrnrlat>90.: urcrnrlat=90. - urcrnrlon= xmax+x_edge - if urcrnrlon>180.: urcrnrlon=180. - area_thresh=1000 - resolution='l' - projection='merc' - lat_0=y[int(len(y)/2)] - lat_1=lat_0-y_edge - lon_0=x[int(len(x)/2)] - #lon_1=110. - rsphere=(6378137.00,6356752.3142) - pd=y_edge - md=x_edge - return urcrnrlon,urcrnrlat,llcrnrlat,llcrnrlon,pd,md,lat_0,lat_1,lon_0,area_thresh,projection,resolution,rsphere - -def map_regions(map_region='default',projection=None,coords=None,m=None,MapPar=None): + ymin = min(y) + ymax = max(y) + y_edge = abs(ymax - ymin) * 0.75 + xmin = min(x) + xmax = max(x) + x_edge = abs(xmax - xmin) * 0.75 + + llcrnrlat = ymin - y_edge + if llcrnrlat < -90.: llcrnrlat = -90. + llcrnrlon = xmin - x_edge + if llcrnrlon < -180.: llcrnrlon = -180. + urcrnrlat = ymax + y_edge + if urcrnrlat > 90.: urcrnrlat = 90. + urcrnrlon = xmax + x_edge + if urcrnrlon > 180.: urcrnrlon = 180. + area_thresh = 1000 + resolution = 'l' + projection = 'merc' + lat_0 = y[int(len(y) / 2)] + lat_1 = lat_0 - y_edge + lon_0 = x[int(len(x) / 2)] + # lon_1=110. + rsphere = (6378137.00, 6356752.3142) + pd = y_edge + md = x_edge + return urcrnrlon, urcrnrlat, llcrnrlat, llcrnrlon, pd, md, lat_0, lat_1, lon_0, area_thresh, projection, resolution, rsphere + + +def map_regions(map_region='default', projection=None, coords=None, m=None, MapPar=None): """ This is a core function, used throughout this module. It is called by the :func:`get_FIGURE` function. If you want to create a new @@ -269,662 +282,663 @@ def map_regions(map_region='default',projection=None,coords=None,m=None,MapPar=N """ - #set some default values - MP=Structure() + # set some default values + MP = Structure() if MapPar: MP.set_with_dict(MapPar) MapPar = MP else: - MapPar=Structure() - FigPar=Structure() + MapPar = Structure() + FigPar = Structure() - MapPar.anchor='C' - FigPar.figsize=(8,7) #w,h tuple + MapPar.anchor = 'C' + FigPar.figsize = (8, 7) # w,h tuple # rect = l,b,w,h - FigPar.axlocs=[0.05,0.01,.8,.9] + FigPar.axlocs = [0.05, 0.01, .8, .9] - if map_region==0: + if map_region == 0: if coords == None: - print "ERROR: for region 0 must provide coordinate tuple (x,y)" - raise IOError - MapPar.urcrnrlon,MapPar.urcrnrlat,MapPar.llcrnrlat,MapPar.llcrnrlon,\ - pd,md,MapPar.lat_0,MapPar.lat_1,MapPar.lon_0,\ - MapPar.area_thresh,MapPar.projection,MapPar.resolution,\ - MapPar.rsphere=map_dynaMap(coords) - - if map_region == None or map_region in ['NorthernHemisphere','default']: + print + "ERROR: for region 0 must provide coordinate tuple (x,y)" + raise IOError() + MapPar.urcrnrlon, MapPar.urcrnrlat, MapPar.llcrnrlat, MapPar.llcrnrlon, \ + pd, md, MapPar.lat_0, MapPar.lat_1, MapPar.lon_0, \ + MapPar.area_thresh, MapPar.projection, MapPar.resolution, \ + MapPar.rsphere = map_dynaMap(*coords) + + if map_region == None or map_region in ['NorthernHemisphere', 'default']: "NorthernHemisphere is default" - MapPar.projection='cyl' - MapPar.llcrnrlat=0 - MapPar.urcrnrlat=90 - MapPar.llcrnrlon=-180 - MapPar.urcrnrlon=180 - MapPar.resolution='c' - #MapPar.anchor='W' - FigPar.figsize=(8,3) #w,h tuple - FigPar.axlocs=[0.1,0.1,.7,.8] + MapPar.projection = 'cyl' + MapPar.llcrnrlat = 0 + MapPar.urcrnrlat = 90 + MapPar.llcrnrlon = -180 + MapPar.urcrnrlon = 180 + MapPar.resolution = 'c' + # MapPar.anchor='W' + FigPar.figsize = (8, 3) # w,h tuple + FigPar.axlocs = [0.1, 0.1, .7, .8] if map_region == 'EarthMBTFPQ': "Most of the Earth" - MapPar.projection='mbtfpq' - #MapPar.llcrnrlat=-40 - #MapPar.urcrnrlat=85 - #MapPar.llcrnrlon=-180 - #MapPar.urcrnrlon=180 - MapPar.resolution='c' - MapPar.lon_0=0 - #MapPar.anchor='W' - FigPar.figsize=(8,7) #w,h tuple - FigPar.axlocs=[0.05,0.01,.8,.96] + MapPar.projection = 'mbtfpq' + # MapPar.llcrnrlat=-40 + # MapPar.urcrnrlat=85 + # MapPar.llcrnrlon=-180 + # MapPar.urcrnrlon=180 + MapPar.resolution = 'c' + MapPar.lon_0 = 0 + # MapPar.anchor='W' + FigPar.figsize = (8, 7) # w,h tuple + FigPar.axlocs = [0.05, 0.01, .8, .96] if map_region == 'EarthMerc': "Most of the Earth" - MapPar.projection='merc' - MapPar.llcrnrlat=-40 - MapPar.urcrnrlat=85 - MapPar.llcrnrlon=-180 - MapPar.urcrnrlon=180 - MapPar.resolution='c' - MapPar.lat_ts=20 - #MapPar.anchor='W' - FigPar.figsize=(8,7) #w,h tuple - FigPar.axlocs=[0.05,0.01,.8,.96] + MapPar.projection = 'merc' + MapPar.llcrnrlat = -40 + MapPar.urcrnrlat = 85 + MapPar.llcrnrlon = -180 + MapPar.urcrnrlon = 180 + MapPar.resolution = 'c' + MapPar.lat_ts = 20 + # MapPar.anchor='W' + FigPar.figsize = (8, 7) # w,h tuple + FigPar.axlocs = [0.05, 0.01, .8, .96] if map_region == 'NH': "NorthernHemisphere is default" - MapPar.projection='merc' - MapPar.llcrnrlat=-40 - MapPar.urcrnrlat=85 - MapPar.llcrnrlon=-180 - MapPar.urcrnrlon=180 - MapPar.resolution='c' - MapPar.lat_ts=20 - #MapPar.anchor='W' - FigPar.figsize=(8,7) #w,h tuple - FigPar.axlocs=[0.05,0.01,.8,.96] - - elif map_region in ['Globe','Earth']: + MapPar.projection = 'merc' + MapPar.llcrnrlat = -40 + MapPar.urcrnrlat = 85 + MapPar.llcrnrlon = -180 + MapPar.urcrnrlon = 180 + MapPar.resolution = 'c' + MapPar.lat_ts = 20 + # MapPar.anchor='W' + FigPar.figsize = (8, 7) # w,h tuple + FigPar.axlocs = [0.05, 0.01, .8, .96] + + elif map_region in ['Globe', 'Earth']: """ Whole earth, mercator """ - MapPar.projection='merc' - MapPar.llcrnrlat=-75 - MapPar.urcrnrlat=85 - MapPar.llcrnrlon=-180 - MapPar.urcrnrlon=180 - MapPar.resolution='c' - #MapPar.anchor='W' - #FigPar.figsize=(8,3) #w,h tuple - FigPar.axlocs=[0.05,0.05,.7,.8] + MapPar.projection = 'merc' + MapPar.llcrnrlat = -75 + MapPar.urcrnrlat = 85 + MapPar.llcrnrlon = -180 + MapPar.urcrnrlon = 180 + MapPar.resolution = 'c' + # MapPar.anchor='W' + # FigPar.figsize=(8,3) #w,h tuple + FigPar.axlocs = [0.05, 0.05, .7, .8] elif map_region == 'WholeEarth': " Use Mollweid " - MapPar.projection='moll' - MapPar.lon_0=0 - MapPar.resolution='c' - FigPar.figsize=(6,6) - FigPar.axlocs=[0.1,0.05,0.7,.85] + MapPar.projection = 'moll' + MapPar.lon_0 = 0 + MapPar.resolution = 'c' + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] - elif map_region==1 and projection=='lcc' or map_region=="NorwegianSea": + elif map_region == 1 and projection == 'lcc' or map_region == "NorwegianSea": "Norwegian Sea" - MapPar.llcrnrlat=55. - MapPar.llcrnrlon=-5. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=50. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=60. - MapPar.lon_0=15. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='SOLAS' or map_region=="solas": - MapPar.llcrnrlat=45. - MapPar.llcrnrlon=-125. - MapPar.urcrnrlat=78. - MapPar.urcrnrlon=-40. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.axlocs=[0.1,0.05,0.7,.85] - - elif map_region=='obuoy': - MapPar.llcrnrlat=50. - MapPar.llcrnrlon=-179. - MapPar.urcrnrlat=80. - MapPar.urcrnrlon=-119. - MapPar.area_thresh=1000. - MapPar.resolution='i' - MapPar.projection='merc' - MapPar.lat_0=70. - MapPar.lon_0=-150. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='NORTHSEA' or map_region=="NorthSea": - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-55. - MapPar.urcrnrlat=78. - MapPar.urcrnrlon=50. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.axlocs=[0.1,0.05,0.7,.85] - - elif map_region=='Asia': - MapPar.llcrnrlat=0. - MapPar.llcrnrlon=60. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=180. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=20. - MapPar.lon_0=100. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='Japan': - MapPar.llcrnrlat=25. - MapPar.llcrnrlon=125. - MapPar.urcrnrlat=45. - MapPar.urcrnrlon=160. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=35. - MapPar.lon_0=140. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(8,7) #w,h tuple + MapPar.llcrnrlat = 55. + MapPar.llcrnrlon = -5. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = 50. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 60. + MapPar.lon_0 = 15. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'SOLAS' or map_region == "solas": + MapPar.llcrnrlat = 45. + MapPar.llcrnrlon = -125. + MapPar.urcrnrlat = 78. + MapPar.urcrnrlon = -40. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] + + elif map_region == 'obuoy': + MapPar.llcrnrlat = 50. + MapPar.llcrnrlon = -179. + MapPar.urcrnrlat = 80. + MapPar.urcrnrlon = -119. + MapPar.area_thresh = 1000. + MapPar.resolution = 'i' + MapPar.projection = 'merc' + MapPar.lat_0 = 70. + MapPar.lon_0 = -150. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'NORTHSEA' or map_region == "NorthSea": + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -55. + MapPar.urcrnrlat = 78. + MapPar.urcrnrlon = 50. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] + + elif map_region == 'Asia': + MapPar.llcrnrlat = 0. + MapPar.llcrnrlon = 60. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = 180. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 20. + MapPar.lon_0 = 100. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'Japan': + MapPar.llcrnrlat = 25. + MapPar.llcrnrlon = 125. + MapPar.urcrnrlat = 45. + MapPar.urcrnrlon = 160. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 35. + MapPar.lon_0 = 140. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (8, 7) # w,h tuple # rect = l,b,w,h - FigPar.axlocs=[0.05,0.01,0.8,.99] - - elif map_region=='JapanOld': - MapPar.llcrnrlat=25. - MapPar.llcrnrlon=125. - MapPar.urcrnrlat=45. - MapPar.urcrnrlon=160. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=35. - MapPar.lon_0=140. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='Pacific_old': + FigPar.axlocs = [0.05, 0.01, 0.8, .99] + + elif map_region == 'JapanOld': + MapPar.llcrnrlat = 25. + MapPar.llcrnrlon = 125. + MapPar.urcrnrlat = 45. + MapPar.urcrnrlon = 160. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 35. + MapPar.lon_0 = 140. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'Pacific_old': """Larger Scale Overview of the North Atlantic LAEA""" - MapPar.llcrnrlat=-30. - MapPar.llcrnrlon=130. - MapPar.urcrnrlat=60. - MapPar.urcrnrlon=-80. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='merc' - MapPar.lat_0=40. - MapPar.lat_1=40. - MapPar.lon_0=-160. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,6) - FigPar.axlocs=[0.08,0.05,0.7,.85] - - - elif map_region=='NorthAtlantic': - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=77. - MapPar.urcrnrlon=50. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='NorthAtlanticWide': + MapPar.llcrnrlat = -30. + MapPar.llcrnrlon = 130. + MapPar.urcrnrlat = 60. + MapPar.urcrnrlon = -80. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'merc' + MapPar.lat_0 = 40. + MapPar.lat_1 = 40. + MapPar.lon_0 = -160. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.08, 0.05, 0.7, .85] + + + elif map_region == 'NorthAtlantic': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 77. + MapPar.urcrnrlon = 50. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'NorthAtlanticWide': """Larger Scale Overview of the North Atlantic LAEA""" - MapPar.llcrnrlat=50. - MapPar.llcrnrlon=-120. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=40. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=40. - MapPar.lat_1=20. - MapPar.lon_0=-40. - MapPar.rsphere=(6378137.00,6356752.3142) - - - elif map_region=='EYJ2': - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-25. - MapPar.urcrnrlat=65. - MapPar.urcrnrlon=45. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=15. - MapPar.lon_0=-100. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(5,4) - - elif map_region=='CALNEX': - MapPar.llcrnrlat=33. - MapPar.llcrnrlon=-125. - MapPar.urcrnrlat=39. - MapPar.urcrnrlon=-115. - MapPar.area_thresh=1000. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=35. - MapPar.lon_0=-120. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,6) - FigPar.axlocs=[0.05,0.05,0.7,.85] - - elif map_region=='EYJ': - MapPar.llcrnrlat=32. - MapPar.llcrnrlon=-125. - MapPar.urcrnrlat=40. - MapPar.urcrnrlon=-117. - MapPar.area_thresh=100. - MapPar.resolution='h' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-120. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,6) - FigPar.axlocs=[0.05,0.05,0.8,.85] - - elif map_region=='SAVAAv2': + MapPar.llcrnrlat = 50. + MapPar.llcrnrlon = -120. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 40. + MapPar.lat_1 = 20. + MapPar.lon_0 = -40. + MapPar.rsphere = (6378137.00, 6356752.3142) + + + elif map_region == 'EYJ2': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -25. + MapPar.urcrnrlat = 65. + MapPar.urcrnrlon = 45. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 15. + MapPar.lon_0 = -100. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (5, 4) + + elif map_region == 'CALNEX': + MapPar.llcrnrlat = 33. + MapPar.llcrnrlon = -125. + MapPar.urcrnrlat = 39. + MapPar.urcrnrlon = -115. + MapPar.area_thresh = 1000. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 35. + MapPar.lon_0 = -120. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.05, 0.05, 0.7, .85] + + elif map_region == 'EYJ': + MapPar.llcrnrlat = 32. + MapPar.llcrnrlon = -125. + MapPar.urcrnrlat = 40. + MapPar.urcrnrlon = -117. + MapPar.area_thresh = 100. + MapPar.resolution = 'h' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -120. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.05, 0.05, 0.8, .85] + + elif map_region == 'SAVAAv2': "for Volcanoe eruption, requires a m instance passed in the MapPar Structure" - MapPar.projection="geos" - MapPar.resolution='l' - MapPar.lon_0=6. - MapPar.llcrnrx=0 - MapPar.llcrnry=0 - #In this case a Basemap instance is required - #in advance - MapPar.urcrnrx=MapPar.m.urcrnrx/2. - MapPar.urcrnry=MapPar.m.urcrnry/2. - - - elif map_region=='VAUUAV_MODIS': + MapPar.projection = "geos" + MapPar.resolution = 'l' + MapPar.lon_0 = 6. + MapPar.llcrnrx = 0 + MapPar.llcrnry = 0 + # In this case a Basemap instance is required + # in advance + MapPar.urcrnrx = MapPar.m.urcrnrx / 2. + MapPar.urcrnry = MapPar.m.urcrnry / 2. + + + elif map_region == 'VAUUAV_MODIS': "Zoom over Kongsfjorden" - MapPar.llcrnrlat=78.8 - MapPar.llcrnrlon=11.6 - MapPar.urcrnrlat=79.2 - MapPar.urcrnrlon=12.6 - MapPar.area_thresh=10. - MapPar.resolution='h' - MapPar.projection='laea' - MapPar.lat_ts=79. - MapPar.lon_0=12. - MapPar.lat_0=79. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='VAUUAV09': + MapPar.llcrnrlat = 78.8 + MapPar.llcrnrlon = 11.6 + MapPar.urcrnrlat = 79.2 + MapPar.urcrnrlon = 12.6 + MapPar.area_thresh = 10. + MapPar.resolution = 'h' + MapPar.projection = 'laea' + MapPar.lat_ts = 79. + MapPar.lon_0 = 12. + MapPar.lat_0 = 79. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'VAUUAV09': "Zoom over Kongsfjorden" - MapPar.llcrnrlat=78.87 - MapPar.llcrnrlon=11.7 - MapPar.urcrnrlat=78.98 - MapPar.urcrnrlon=12.5 - MapPar.area_thresh=10. - MapPar.resolution='h' - MapPar.projection='lcc' - MapPar.lat_1=75. - MapPar.lon_0=20. - MapPar.rsphere=(6378137.00,6356752.3142) + MapPar.llcrnrlat = 78.87 + MapPar.llcrnrlon = 11.7 + MapPar.urcrnrlat = 78.98 + MapPar.urcrnrlon = 12.5 + MapPar.area_thresh = 10. + MapPar.resolution = 'h' + MapPar.projection = 'lcc' + MapPar.lat_1 = 75. + MapPar.lon_0 = 20. + MapPar.rsphere = (6378137.00, 6356752.3142) pd = 2 md = 5 - elif map_region=='alaska': + elif map_region == 'alaska': "Over alaska" - MapPar.llcrnrlat=30. - MapPar.llcrnrlon=-185. - MapPar.urcrnrlat=60. - MapPar.urcrnrlon=-60. - MapPar.area_thresh=1000. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_0=50. - MapPar.lon_0=-135. - MapPar.rsphere=(6378137.00,6356752.3142) - - - elif map_region=='GreenlandPoly': + MapPar.llcrnrlat = 30. + MapPar.llcrnrlon = -185. + MapPar.urcrnrlat = 60. + MapPar.urcrnrlon = -60. + MapPar.area_thresh = 1000. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_0 = 50. + MapPar.lon_0 = -135. + MapPar.rsphere = (6378137.00, 6356752.3142) + + + elif map_region == 'GreenlandPoly': "Over Greenland" - MapPar.llcrnrlat=50. - MapPar.llcrnrlon=-55. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=-21. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='poly' - MapPar.lat_0=70. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='CICCI': + MapPar.llcrnrlat = 50. + MapPar.llcrnrlon = -55. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = -21. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'poly' + MapPar.lat_0 = 70. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'CICCI': "Zoom over Kongsfjorden" - MapPar.llcrnrlat=78.7 - MapPar.llcrnrlon=11 - MapPar.urcrnrlat=79.8 - MapPar.urcrnrlon=14.5 - MapPar.area_thresh=10. - MapPar.resolution='h' - MapPar.projection='lcc' - MapPar.lat_1=75. - MapPar.lon_0=20. - MapPar.rsphere=(6378137.00,6356752.3142) + MapPar.llcrnrlat = 78.7 + MapPar.llcrnrlon = 11 + MapPar.urcrnrlat = 79.8 + MapPar.urcrnrlon = 14.5 + MapPar.area_thresh = 10. + MapPar.resolution = 'h' + MapPar.projection = 'lcc' + MapPar.lat_1 = 75. + MapPar.lon_0 = 20. + MapPar.rsphere = (6378137.00, 6356752.3142) pd = 2 md = 5 - elif map_region=='Svalbard': + elif map_region == 'Svalbard': "Over Svalbard" MapPar.update({'area_thresh': 10.0, - 'lat_0': 80.0, - 'llcrnrlat': 72.0, - 'llcrnrlon': -5.0, - 'lon_0': 20.0, - 'projection': 'lcc', - 'resolution': 'h', - 'rsphere': (6378137.0, 6356752.3142), - 'urcrnrlat': 80.0, - 'urcrnrlon': 45.0}) - - elif map_region=='stads-ortho': + 'lat_0': 80.0, + 'llcrnrlat': 72.0, + 'llcrnrlon': -5.0, + 'lon_0': 20.0, + 'projection': 'lcc', + 'resolution': 'h', + 'rsphere': (6378137.0, 6356752.3142), + 'urcrnrlat': 80.0, + 'urcrnrlon': 45.0}) + + elif map_region == 'stads-ortho': "Stads nested output" MapPar.update({ - 'lat_0': 70.0, - 'lon_0': 10.0, - 'projection': 'ortho', - 'resolution': 'l' - }) + 'lat_0': 70.0, + 'lon_0': 10.0, + 'projection': 'ortho', + 'resolution': 'l' + }) - elif map_region=='stads-nest': + elif map_region == 'stads-nest': "Stads nested output" MapPar.update({ - 'area_thresh' : 1000., - 'lat_0': 60.0, - 'llcrnrlat': 70.0, - 'llcrnrlon': -25.0, - 'lon_0': 10.0, - 'projection': 'merc', - 'resolution': 'i', - 'rsphere': (6378137.0, 6356752.3142), - 'urcrnrlat': 85.0, - 'urcrnrlon': 35.0}) - - elif map_region=='GreenlandLaea': + 'area_thresh': 1000., + 'lat_0': 60.0, + 'llcrnrlat': 70.0, + 'llcrnrlon': -25.0, + 'lon_0': 10.0, + 'projection': 'merc', + 'resolution': 'i', + 'rsphere': (6378137.0, 6356752.3142), + 'urcrnrlat': 85.0, + 'urcrnrlon': 35.0}) + + elif map_region == 'GreenlandLaea': "Over Greenland" - MapPar.llcrnrlat=57. - MapPar.llcrnrlon=-55. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=-3. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='laea' - MapPar.lat_0=70. - MapPar.lon_0=-42. - MapPar.rsphere=(6378137.00,6356752.3142) - - - elif map_region=='CentralGreenland': + MapPar.llcrnrlat = 57. + MapPar.llcrnrlon = -55. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = -3. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'laea' + MapPar.lat_0 = 70. + MapPar.lon_0 = -42. + MapPar.rsphere = (6378137.00, 6356752.3142) + + + elif map_region == 'CentralGreenland': "Centered around Summit, Greenland" - MapPar.llcrnrlat=68. - MapPar.llcrnrlon=-55. - MapPar.urcrnrlat=74. - MapPar.urcrnrlon=-15 - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=70. - MapPar.lon_0=-40. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='GEOSummit': + MapPar.llcrnrlat = 68. + MapPar.llcrnrlon = -55. + MapPar.urcrnrlat = 74. + MapPar.urcrnrlon = -15 + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 70. + MapPar.lon_0 = -40. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'GEOSummit': "Zoom near Summit, Greenland" - MapPar.llcrnrlat=72. - MapPar.llcrnrlon=-41 - MapPar.urcrnrlat=73. - MapPar.urcrnrlon=-37. - MapPar.area_thresh=1. - MapPar.resolution='h' - MapPar.projection='lcc' - MapPar.lat_1=70. - MapPar.lon_0=-40. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='GEOSummitZoom': + MapPar.llcrnrlat = 72. + MapPar.llcrnrlon = -41 + MapPar.urcrnrlat = 73. + MapPar.urcrnrlon = -37. + MapPar.area_thresh = 1. + MapPar.resolution = 'h' + MapPar.projection = 'lcc' + MapPar.lat_1 = 70. + MapPar.lon_0 = -40. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'GEOSummitZoom': "Zoom near Summit, Greenland" - MapPar.llcrnrlat=72.55 - MapPar.llcrnrlon=-38.52 - MapPar.urcrnrlat=72.6 - MapPar.urcrnrlon=-38.4 - MapPar.area_thresh=10. - MapPar.resolution='h' - MapPar.projection='lcc' - MapPar.lat_1=70. - MapPar.lon_0=-38. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region==2 and projection=='lcc': + MapPar.llcrnrlat = 72.55 + MapPar.llcrnrlon = -38.52 + MapPar.urcrnrlat = 72.6 + MapPar.urcrnrlon = -38.4 + MapPar.area_thresh = 10. + MapPar.resolution = 'h' + MapPar.projection = 'lcc' + MapPar.lat_1 = 70. + MapPar.lon_0 = -38. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 2 and projection == 'lcc': "Whole North Atlantic" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=50. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='north-atlantic': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = 50. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'north-atlantic': "Whole North Atlantic" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=75. - MapPar.urcrnrlon=50. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='finland': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 75. + MapPar.urcrnrlon = 50. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'finland': "North Norway" - MapPar.llcrnrlat=59. - MapPar.llcrnrlon=15. - MapPar.urcrnrlat=71. - MapPar.urcrnrlon=40. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=65. - MapPar.lon_0=25. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='Europe': + MapPar.llcrnrlat = 59. + MapPar.llcrnrlon = 15. + MapPar.urcrnrlat = 71. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 65. + MapPar.lon_0 = 25. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'Europe': "Europe Mercator" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-20. - MapPar.urcrnrlat=68. - MapPar.urcrnrlon=40. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='merc' - MapPar.lat_1=60. - MapPar.lon_0=10. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,6) - FigPar.axlocs=[0.08,0.05,0.7,.85] - - elif map_region=='Europe_ortho': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -20. + MapPar.urcrnrlat = 68. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'merc' + MapPar.lat_1 = 60. + MapPar.lon_0 = 10. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.08, 0.05, 0.7, .85] + + elif map_region == 'Europe_ortho': "Europe ORTHO" - #MapPar.llcrnrlat=35. - #MapPar.llcrnrlon=-20. - #MapPar.urcrnrlat=68. - #MapPar.urcrnrlon=40. - MapPar.area_thresh=10000. - MapPar.resolution='c' - MapPar.projection='ortho' - MapPar.lat_0=60. - MapPar.lon_0=10. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,6) - FigPar.axlocs=[0.05,0.05,0.8,.85] - - elif map_region==3 and projection=='lcc': + # MapPar.llcrnrlat=35. + # MapPar.llcrnrlon=-20. + # MapPar.urcrnrlat=68. + # MapPar.urcrnrlon=40. + MapPar.area_thresh = 10000. + MapPar.resolution = 'c' + MapPar.projection = 'ortho' + MapPar.lat_0 = 60. + MapPar.lon_0 = 10. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.05, 0.05, 0.8, .85] + + elif map_region == 3 and projection == 'lcc': "North Norway" - MapPar.llcrnrlat=65. - MapPar.llcrnrlon=10. - MapPar.urcrnrlat=72. - MapPar.urcrnrlon=30. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=60. - MapPar.lon_0=15. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region==4 and projection=='lcc': + MapPar.llcrnrlat = 65. + MapPar.llcrnrlon = 10. + MapPar.urcrnrlat = 72. + MapPar.urcrnrlon = 30. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 60. + MapPar.lon_0 = 15. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 4 and projection == 'lcc': "Greenland Sea, Fram Strait" - MapPar.llcrnrlat=72. - MapPar.llcrnrlon=-20. - MapPar.urcrnrlat=80. - MapPar.urcrnrlon=30. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_1=70. - MapPar.lon_0=0. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region==4 and projection=='aea': + MapPar.llcrnrlat = 72. + MapPar.llcrnrlon = -20. + MapPar.urcrnrlat = 80. + MapPar.urcrnrlon = 30. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_1 = 70. + MapPar.lon_0 = 0. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 4 and projection == 'aea': "Greenland Sea, Fram Strait" - MapPar.llcrnrlat=72. - MapPar.llcrnrlon=-20. - MapPar.urcrnrlat=80. - MapPar.urcrnrlon=30. - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='lcc' - MapPar.lat_0=70. - MapPar.lon_0=0. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region==5 and projection=='lcc': + MapPar.llcrnrlat = 72. + MapPar.llcrnrlon = -20. + MapPar.urcrnrlat = 80. + MapPar.urcrnrlon = 30. + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'lcc' + MapPar.lat_0 = 70. + MapPar.lon_0 = 0. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 5 and projection == 'lcc': "Whole North Atlantic and Norwegian, Greenland Seas" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=81. - MapPar.urcrnrlon=40. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='lcc' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - FigPar.figsize=(6,5) - FigPar.axlocs=[0.08,0.05,0.7,.8] - - elif map_region=='Pacific': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 81. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'lcc' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + FigPar.figsize = (6, 5) + FigPar.axlocs = [0.08, 0.05, 0.7, .8] + + elif map_region == 'Pacific': """National Geospatial-Intelligence Agency: DMANC1 North America""" - MapPar.llcrnrlat=5. - MapPar.llcrnrlon=-200. - MapPar.urcrnrlat=55. - MapPar.urcrnrlon=-110. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=40. - MapPar.lat_1=20. - MapPar.lon_0=-170. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='DMANC1': + MapPar.llcrnrlat = 5. + MapPar.llcrnrlon = -200. + MapPar.urcrnrlat = 55. + MapPar.urcrnrlon = -110. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 40. + MapPar.lat_1 = 20. + MapPar.lon_0 = -170. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'DMANC1': """National Geospatial-Intelligence Agency: DMANC1 North America""" - MapPar.llcrnrlat=5. - MapPar.llcrnrlon=-130. - MapPar.urcrnrlat=55. - MapPar.urcrnrlon=-10. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=40. - MapPar.lat_1=20. - MapPar.lon_0=-90. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='NorthPole': + MapPar.llcrnrlat = 5. + MapPar.llcrnrlon = -130. + MapPar.urcrnrlat = 55. + MapPar.urcrnrlon = -10. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 40. + MapPar.lat_1 = 20. + MapPar.lon_0 = -90. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'NorthPole': "Whole North Atlantic and Norwegian, Greenland Seas" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=81. - MapPar.urcrnrlon=40. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='npstere' - MapPar.lat_1=5. - MapPar.lon_0=-80. - MapPar.rsphere=(6378137.00,6356752.3142) - MapPar.boundinglat=65. - - elif map_region=='npolar': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 81. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'npstere' + MapPar.lat_1 = 5. + MapPar.lon_0 = -80. + MapPar.rsphere = (6378137.00, 6356752.3142) + MapPar.boundinglat = 65. + + elif map_region == 'npolar': "Same as POLARCAT" - MapPar.area_thresh=1000. - MapPar.resolution='i' - MapPar.projection='npstere' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - MapPar.boundinglat=25. - MapPar.anchor='W' - FigPar.figsize=(6,6) - FigPar.axlocs=[0.1,0.05,0.7,.85] - #FigPar.axlocs = [0.1,0.05,0.9,0.8] - - - elif map_region=='npolar_70': + MapPar.area_thresh = 1000. + MapPar.resolution = 'i' + MapPar.projection = 'npstere' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + MapPar.boundinglat = 25. + MapPar.anchor = 'W' + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] + # FigPar.axlocs = [0.1,0.05,0.9,0.8] + + + elif map_region == 'npolar_70': "Zoom on polar project 70N" - MapPar.area_thresh=1000. - MapPar.resolution='i' - MapPar.projection='npstere' - MapPar.lat_1=5. - MapPar.lon_0=-80. - MapPar.rsphere=(6378137.00,6356752.3142) - MapPar.boundinglat=65. - MapPar.anchor='W' - FigPar.figsize=(6,6) - FigPar.axlocs=[0.1,0.05,0.7,.85] - - - elif map_region=='POLARCAT': + MapPar.area_thresh = 1000. + MapPar.resolution = 'i' + MapPar.projection = 'npstere' + MapPar.lat_1 = 5. + MapPar.lon_0 = -80. + MapPar.rsphere = (6378137.00, 6356752.3142) + MapPar.boundinglat = 65. + MapPar.anchor = 'W' + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] + + + elif map_region == 'POLARCAT': "Whole North Atlantic and Norwegian, Greenland Seas" - MapPar.llcrnrlat=35. - MapPar.llcrnrlon=-95. - MapPar.urcrnrlat=81. - MapPar.urcrnrlon=40. - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='npstere' - MapPar.lat_1=5. - MapPar.lon_0=-20. - MapPar.rsphere=(6378137.00,6356752.3142) - MapPar.boundinglat=45. - MapPar.anchor='W' - FigPar.figsize=(6,6) - FigPar.axlocs=[0.1,0.05,0.7,.85] - #FigPar.axlocs = [0.1,0.05,0.9,0.8] - - #elif map_region=='NorthernHemisphere': + MapPar.llcrnrlat = 35. + MapPar.llcrnrlon = -95. + MapPar.urcrnrlat = 81. + MapPar.urcrnrlon = 40. + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'npstere' + MapPar.lat_1 = 5. + MapPar.lon_0 = -20. + MapPar.rsphere = (6378137.00, 6356752.3142) + MapPar.boundinglat = 45. + MapPar.anchor = 'W' + FigPar.figsize = (6, 6) + FigPar.axlocs = [0.1, 0.05, 0.7, .85] + # FigPar.axlocs = [0.1,0.05,0.9,0.8] + + # elif map_region=='NorthernHemisphere': # """Northern Hemisphere""" # MapPar.llcrnrlat=0. # MapPar.llcrnrlon=-180. @@ -938,100 +952,101 @@ def map_regions(map_region='default',projection=None,coords=None,m=None,MapPar=N # MapPar.lon_0=-110. # MapPar.lon_1=110. # MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='met.no:c_map1': + + elif map_region == 'met.no:c_map1': """Met.no Sea Ice Chart: c_map1 Norwegian / Greenland Seas""" - MapPar.llcrnrlat=57.6 - MapPar.llcrnrlon=-27.3 - MapPar.urcrnrlat=68.27 - MapPar.urcrnrlon=86.5 - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=90. - MapPar.lat_1=75. - MapPar.lon_0=0. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='met.no:cust1': + MapPar.llcrnrlat = 57.6 + MapPar.llcrnrlon = -27.3 + MapPar.urcrnrlat = 68.27 + MapPar.urcrnrlon = 86.5 + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 90. + MapPar.lat_1 = 75. + MapPar.lon_0 = 0. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'met.no:cust1': """Met.no Sea Ice Chart: c_map1 Norwegian / Greenland Seas""" - MapPar.llcrnrlat=74.74 - MapPar.llcrnrlon=-14.24 - MapPar.urcrnrlat=80.44 - MapPar.urcrnrlon=10.95 - MapPar.area_thresh=1000. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=90. - MapPar.lat_1=75. - MapPar.lon_0=0. - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='GRACE': + MapPar.llcrnrlat = 74.74 + MapPar.llcrnrlon = -14.24 + MapPar.urcrnrlat = 80.44 + MapPar.urcrnrlon = 10.95 + MapPar.area_thresh = 1000. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 90. + MapPar.lat_1 = 75. + MapPar.lon_0 = 0. + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'GRACE': """Based around GRACE Campaigns""" - MapPar.llcrnrlat=60.0 - MapPar.llcrnrlon=-70.0 - MapPar.urcrnrlat=75.0 - MapPar.urcrnrlon=-30.00 - MapPar.area_thresh=100. - MapPar.resolution='i' - MapPar.projection='aea' - MapPar.lat_0=60. - MapPar.lat_1=60. - MapPar.lon_0=-55. - MapPar.rsphere=(6378137.00,6356752.3142) - pd=2 - md=5 - - elif map_region=='falcon': + MapPar.llcrnrlat = 60.0 + MapPar.llcrnrlon = -70.0 + MapPar.urcrnrlat = 75.0 + MapPar.urcrnrlon = -30.00 + MapPar.area_thresh = 100. + MapPar.resolution = 'i' + MapPar.projection = 'aea' + MapPar.lat_0 = 60. + MapPar.lat_1 = 60. + MapPar.lon_0 = -55. + MapPar.rsphere = (6378137.00, 6356752.3142) + pd = 2 + md = 5 + + elif map_region == 'falcon': """Based around falcon flights""" - MapPar.llcrnrlat=48.0 - MapPar.llcrnrlon=-15.0 - MapPar.urcrnrlat=54.0 - MapPar.urcrnrlon=-5.00 - MapPar.area_thresh=100. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=50. - MapPar.lat_1=50. - MapPar.lon_0=-10. - MapPar.rsphere=(6378137.00,6356752.3142) - pd=2 - md=5 - - elif map_region=='EUCAARI': + MapPar.llcrnrlat = 48.0 + MapPar.llcrnrlon = -15.0 + MapPar.urcrnrlat = 54.0 + MapPar.urcrnrlon = -5.00 + MapPar.area_thresh = 100. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 50. + MapPar.lat_1 = 50. + MapPar.lon_0 = -10. + MapPar.rsphere = (6378137.00, 6356752.3142) + pd = 2 + md = 5 + + elif map_region == 'EUCAARI': """EUCAARI Domain""" - MapPar.llcrnrlat=38.0 - MapPar.llcrnrlon=-15.0 - MapPar.urcrnrlat=68.0 - MapPar.urcrnrlon=30.00 - MapPar.area_thresh=100. - MapPar.resolution='l' - MapPar.projection='aea' - MapPar.lat_0=50. - MapPar.lat_1=50. - MapPar.lon_0=7.5 - MapPar.rsphere=(6378137.00,6356752.3142) - - elif map_region=='auto' and coords != None: + MapPar.llcrnrlat = 38.0 + MapPar.llcrnrlon = -15.0 + MapPar.urcrnrlat = 68.0 + MapPar.urcrnrlon = 30.00 + MapPar.area_thresh = 100. + MapPar.resolution = 'l' + MapPar.projection = 'aea' + MapPar.lat_0 = 50. + MapPar.lat_1 = 50. + MapPar.lon_0 = 7.5 + MapPar.rsphere = (6378137.00, 6356752.3142) + + elif map_region == 'auto' and coords != None: """Based around falcon flights""" - MapPar.urcrnrlon,MapPar.urcrnrlat,MapPar.llcrnrlat,MapPar.llcrnrlon,\ - pd,md,MapPar.lat_0,MapPar.lat_1,MapPar.lon_0,\ - MapPar.area_thresh,MapPar.projection,\ - MapPar.resolution,MapPar.rsphere = map_dynaMap(coords) - - #MapPar.urcrnrlon,MapPar.urcrnrlat,MapPar.llcrnrlat,MapPar.llcrnrlon,\ - # pd,md,MapPar.lat_0,MapPar.lat_1,MapPar.lon_0,\ - #MapPar.area_thresh,MapPar.projection,MapPar.resolution,\ - #MapPar.rsphere=map_dynaMap(([-180,180],[-90,90])) - #MapPar.boundinglat=50. - try: del MapPar['m'] #in case - except: pass - return MapPar,FigPar - - - -def draw_grid(m,xdiv=10.,ydiv=5.,location=[1,0,0,1],\ + MapPar.urcrnrlon, MapPar.urcrnrlat, MapPar.llcrnrlat, MapPar.llcrnrlon, \ + pd, md, MapPar.lat_0, MapPar.lat_1, MapPar.lon_0, \ + MapPar.area_thresh, MapPar.projection, \ + MapPar.resolution, MapPar.rsphere = map_dynaMap(*coords) + + # MapPar.urcrnrlon,MapPar.urcrnrlat,MapPar.llcrnrlat,MapPar.llcrnrlon,\ + # pd,md,MapPar.lat_0,MapPar.lat_1,MapPar.lon_0,\ + # MapPar.area_thresh,MapPar.projection,MapPar.resolution,\ + # MapPar.rsphere=map_dynaMap(([-180,180],[-90,90])) + # MapPar.boundinglat=50. + try: + del MapPar['m'] # in case + except: + pass + return MapPar, FigPar + + +def draw_grid(m, xdiv=10., ydiv=5., location=[1, 0, 0, 1], \ linewidth=0.5, color='k'): """ draw parallels and meridians on a map. @@ -1052,66 +1067,66 @@ def draw_grid(m,xdiv=10.,ydiv=5.,location=[1,0,0,1],\ """ - #Label properties + # Label properties p_leg = mpl.font_manager.FontProperties(size='8') - - pd_options = [0.0001,0.001,0.01, 0.1, 0.2, 0.5, 1, 2.5, 5, 10, 20] - md_options = [0.0001,0.001,0.01, 0.1, 0.2, 0.5, 1, 2.5, 5, 10, 20] - + + pd_options = [0.0001, 0.001, 0.01, 0.1, 0.2, 0.5, 1, 2.5, 5, 10, 20] + md_options = [0.0001, 0.001, 0.01, 0.1, 0.2, 0.5, 1, 2.5, 5, 10, 20] + xdiff = np.abs(m.urcrnrlon - m.llcrnrlon) ydiff = np.abs(m.urcrnrlat - m.llcrnrlat) - - if m.projection in ['npstere','spstere']: + + if m.projection in ['npstere', 'spstere']: ydiff = 90. - np.abs(m.urcrnrlat) maxlat = 90. else: maxlat = None - + # setup map to have meridians: if xdiff > xdiv: - md = np.round((xdiff / xdiv )) + md = np.round((xdiff / xdiv)) else: - md = np.round((xdiff / xdiv ),1) - + md = np.round((xdiff / xdiv), 1) + md = md_options[(np.abs(np.array(md_options) - md)).argmin()] - #print md - if m.projection in ['npstere','spstere']: - meridians = np.arange(-180,181,20) + # print(md) + if m.projection in ['npstere', 'spstere']: + meridians = np.arange(-180, 181, 20) else: - meridians = np.arange(m.llcrnrlon,m.urcrnrlon,md) - MD = m.drawmeridians(meridians,labels=location, - linewidth=linewidth,color=color,fontproperties=p_leg) - - + meridians = np.arange(m.llcrnrlon, m.urcrnrlon, md) + MD = m.drawmeridians(meridians, labels=location, + linewidth=linewidth, color=color, fontproperties=p_leg) + # setup map to have parallels. if ydiff > ydiv: - pd = np.round((ydiff / ydiv )) + pd = np.round((ydiff / ydiv)) else: - pd = np.round((ydiff / ydiv ), 2) + pd = np.round((ydiff / ydiv), 2) if not maxlat: maxlat = m.urcrnrlat - - pd = pd_options[(np.abs(np.array(pd_options) - pd)).argmin()] - #print pd - if m.projection in ['npstere','spstere']: - parallels = np.arange(30,91,10) + + pd = pd_options[(np.abs(np.array(pd_options) - pd)).argmin()] + # print(pd) + if m.projection in ['npstere', 'spstere']: + parallels = np.arange(30, 91, 10) else: if pd > 0.01: - parallels = np.arange(np.round(m.llcrnrlat),maxlat,pd) + parallels = np.arange(np.round(m.llcrnrlat), maxlat, pd) else: - parallels = np.arange(m.llcrnrlat,maxlat,pd) + parallels = np.arange(m.llcrnrlat, maxlat, pd) + + print(parallels, m.llcrnrlat, maxlat, pd) + MP = m.drawparallels(parallels, labels=location, + linewidth=linewidth, color=color, fontproperties=p_leg) + return MP, MD - print parallels, m.llcrnrlat, maxlat, pd - MP = m.drawparallels(parallels,labels=location, - linewidth=linewidth,color=color,fontproperties=p_leg) - return MP,MD -def plot_ECMWF(nc,variable,time,level,map_region='NorthAmerica'): +def plot_ECMWF(nc, variable, time, level, map_region='NorthAmerica'): """ uses :function:`plot_grid`_ to plot a netcdf file downloaded from ECMWF MARS data store. """ - if isinstance(nc,str): + if isinstance(nc, str): nc = NetCDFFile(nc) alldata = nc.variables[variable][:] @@ -1120,21 +1135,19 @@ def plot_ECMWF(nc,variable,time,level,map_region='NorthAmerica'): lats = nc.variables['latitude'][:] lat2 = lats[::-1] lons = nc.variables['longitude'][:] - lon1 = np.arange(-180,180,360./len(lons)) - data = alldata[time,level,:,:] - data = data[::-1,:] #flip lats - #shift the grid to go from -180 to 180 - d2,lon2 = shiftgrid(-180,data,lon1) - plot_grid((lon2,lat2,d2),map_region=map_region) - + lon1 = np.arange(-180, 180, 360. / len(lons)) + data = alldata[time, level, :, :] + data = data[::-1, :] # flip lats + # shift the grid to go from -180 to 180 + d2, lon2 = shiftgrid(-180, data, lon1) + plot_grid((lon2, lat2, d2), map_region=map_region) -def get_OpenDAP(analysistime,date_range=None,level_range=None, - URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd%s/gfs_hd_%sz/" , - datalist=['prmslmsl','ugrdprs','vgrdprs'], +def get_OpenDAP(analysistime, date_range=None, level_range=None, + URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd%s/gfs_hd_%sz/", + datalist=['prmslmsl', 'ugrdprs', 'vgrdprs'], ): - - """ Open an OpenDAP connection and return requested data from list + """ Open an OpenDAP connection and return requested data from list Written for use with GFS/NCEP data, but could be modified """ @@ -1149,24 +1162,24 @@ def get_OpenDAP(analysistime,date_range=None,level_range=None, YYYYMMDD = atime[:8] HH = atime[8:10] elif isinstance(analysistime, datetime.datime): - YYYYMMDD = dt.datetime.strftime(analysistime-dt.timedelta(2),'%Y%m%d') + YYYYMMDD = dt.datetime.strftime(analysistime - dt.timedelta(2), '%Y%m%d') HH = str(analysistime.hour).zfill(2) - #YYYYMMDDHH1 = '2010042100' - #YYYYMMDDHH2 = '2010042212' + # YYYYMMDDHH1 = '2010042100' + # YYYYMMDDHH2 = '2010042212' - #YYYY = YYYYMMDDHH1[0:4] - #if YYYY != YYYYMMDDHH2[0:4]: - # raise ValueError,'dates must be in same year' + # YYYY = YYYYMMDDHH1[0:4] + # if YYYY != YYYYMMDDHH2[0:4]: + # raise ValueError('dates must be in same year') # set OpenDAP server URL. - URLbase = URLbase % (YYYYMMDD,HH) + URLbase = URLbase % (YYYYMMDD, HH) try: print("trying to retrieve:\n{0}".format(URLbase)) - #use the basemap NetCDFFile to retrieve OpenDAP + # use the basemap NetCDFFile to retrieve OpenDAP data = mpl_NetCDFFile(URLbase) except: - raise IOError, 'opendap server not providing the requested data' + raise IOError('opendap server not providing the requested data') D['extracted_levels'] = [] D['extracted_dates'] = [] lev_units = data.variables['lev'].units @@ -1178,37 +1191,37 @@ def get_OpenDAP(analysistime,date_range=None,level_range=None, times = data.variables['time'][:] time_units = data.variables['time'].units # convert numeric time values to datetime objects. - fdates = num2date(times,time_units) + fdates = num2date(times, time_units) # put times in YYYYMMDDHH format. dates = [fdate.strftime('%Y%m%d%H%M%S') for fdate in fdates] # find indices bounding desired times. - #ntime1 = dates.index(YYYYMMDDHH1) - #ntime2 = dates.index(YYYYMMDDHH2) - #ntime1 = dates.index(date_range[0]) - #ntime2 = dates.index(date_range[-1]) - #print 'ntime1,ntime2:',ntime1,ntime2 - #if ntime1 >= ntime2: - # raise ValueError,'date2 must be greater than date1' + # ntime1 = dates.index(YYYYMMDDHH1) + # ntime2 = dates.index(YYYYMMDDHH2) + # ntime1 = dates.index(date_range[0]) + # ntime2 = dates.index(date_range[-1]) + # print('ntime1,ntime2:',ntime1,ntime2) + # if ntime1 >= ntime2: + # raise ValueError('date2 must be greater than date1') # get requested data from DODs: for dset in datalist: - lrange,trange = False,False + lrange, trange = False, False D[dset] = {} D[dset]['long_name'] = data.variables[dset].long_name D[dset]['dimensions'] = data.variables[dset].dimensions D[dset]['dtype'] = data.variables[dset].dtype D[dset]['shape'] = data.variables[dset].shape D[dset]['missing_value'] = data.variables[dset].missing_value - #now make some decisions based on data shape + # now make some decisions based on data shape if 'time' in D[dset]['dimensions']: trange = True - if isinstance(date_range,int): + if isinstance(date_range, int): nt0 = dates.index(date_range) nt1 = nt0 + 1 - elif isinstance(date_range,tuple): + elif isinstance(date_range, tuple): nt0 = dates.index(date_range[0]) nt1 = dates.index(date_range[-1]) + 1 - elif isinstance(date_range,list): + elif isinstance(date_range, list): nt0 = dates.index(date_range[0]) nt1 = dates.index(date_range[-1]) + 1 else: @@ -1217,10 +1230,10 @@ def get_OpenDAP(analysistime,date_range=None,level_range=None, if 'lev' in D[dset]['dimensions']: lrange = True - if isinstance(level_range,int): + if isinstance(level_range, int): l0 = level_range l1 = l0 + 1 - elif isinstance(level_range,list): + elif isinstance(level_range, list): l0 = level_range[0] l1 = level_range[-1] else: @@ -1229,18 +1242,15 @@ def get_OpenDAP(analysistime,date_range=None,level_range=None, if lrange: if trange: - D[dset]['data'] = data.variables[dset][nt0:nt1,l0:l1,:,:] + D[dset]['data'] = data.variables[dset][nt0:nt1, l0:l1, :, :] D['extracted_dates'] = dates[nt0:nt1] D['extracted_levels'] = levels[l0:l1] else: if trange: - D[dset]['data'] = data.variables[dset][nt0:nt1,:,:] + D[dset]['data'] = data.variables[dset][nt0:nt1, :, :] D['dates'] = dates[nt0:nt1] - - - D['lev_units'] = lev_units D['levels'] = levels D['dt'] = fdates @@ -1252,19 +1262,19 @@ def get_OpenDAP(analysistime,date_range=None,level_range=None, # get sea level pressure and 10-m wind data. - #slpdata = data.variables['prmslmsl'] - #udata = data.variables['ugrdprs'] - #vdata = data.variables['vgrdprs'] + # slpdata = data.variables['prmslmsl'] + # udata = data.variables['ugrdprs'] + # vdata = data.variables['vgrdprs'] # mult slp by 0.01 to put in units of millibars. - #slpin = 0.01*slpdata[ntime1:ntime2+1,:,:] - #uin = udata[ntime1:ntime2+1,0,:,:] - #vin = vdata[ntime1:ntime2+1,0,:,:] + # slpin = 0.01*slpdata[ntime1:ntime2+1,:,:] + # uin = udata[ntime1:ntime2+1,0,:,:] + # vin = vdata[ntime1:ntime2+1,0,:,:] def default_netcdf(nco_filename, - lon0=-179.5,lat0=-89.5, - nx=720,ny=360, - dx=0.5,dy=0.5, + lon0=-179.5, lat0=-89.5, + nx=720, ny=360, + dx=0.5, dy=0.5, title="", institution="Norwegian Institute for Air Research", source="FLEXPART Model Data", @@ -1304,7 +1314,7 @@ def default_netcdf(nco_filename, """ - nco = NetCDFFile(nco_filename,'w') + nco = NetCDFFile(nco_filename, 'w') nco.author = author nco.createdate = dt.datetime.now().ctime() nco.contact = contact @@ -1316,25 +1326,26 @@ def default_netcdf(nco_filename, nco.references = references nco.history = history - nco.createDimension('lon',nx) - nco.createDimension('lat',ny) - nco.createVariable('lat','d',('lat',)) - nco.createVariable('lon','d',('lon',)) - lon = np.arange(lon0,lon0+(nx*dx),dx) - lat = np.arange(lat0,lat0+(ny*dy),dy) + nco.createDimension('lon', nx) + nco.createDimension('lat', ny) + nco.createVariable('lat', 'd', ('lat',)) + nco.createVariable('lon', 'd', ('lon',)) + lon = np.arange(lon0, lon0 + (nx * dx), dx) + lat = np.arange(lat0, lat0 + (ny * dy), dy) nco.variables['lat'][:] = lat nco.variables['lat'].long_name = 'latitude' - nco.variables['lat'].units= 'degrees_north' + nco.variables['lat'].units = 'degrees_north' nco.variables['lon'][:] = lon nco.variables['lon'].long_name = 'longitude' - nco.variables['lon'].units= 'degrees_east' + nco.variables['lon'].units = 'degrees_east' - #nco.createVariable('data','d',('lon','lat')) + # nco.createVariable('data','d',('lon','lat')) return nco -def grid_to_netcdf(H,D,nco_filename,spc,time): + +def grid_to_netcdf(H, D, nco_filename, spc, time): """INCOMPLETE, NOT YETWORKING Dump a grid, returned from :module:`rf.get_grid` into a netcdf file. @@ -1343,9 +1354,12 @@ def grid_to_netcdf(H,D,nco_filename,spc,time): nco.sync() nco.close() """ - print "Not working" + print + "Not working" return + + """ @@ -1384,11 +1398,10 @@ def grid_to_netcdf(H,D,nco_filename,spc,time): """ -def get_FIGURE(fig=None,ax=None,m=None,map_region=None, - projection=None,coords=None,getm=True, - MapPar=None,FigPar=None, +def get_FIGURE(fig=None, ax=None, m=None, map_region=None, + projection=None, coords=None, getm=True, + MapPar=None, FigPar=None, image=None): - """ This is a core function, used throughout this module. It is called by various functions. The idea is that I create a :class:`Structure` that @@ -1433,26 +1446,26 @@ def get_FIGURE(fig=None,ax=None,m=None,map_region=None, """ FIGURE = Structure() - + if getm: if m == None: if image: - fig,m = get_base_image(image,map_region=map_region, - projection=projection, - coords=coords, - MapPar=MapPar, - FigPar=FigPar, - ) + fig, m = get_base_image(image, map_region=map_region, + projection=projection, + coords=coords, + MapPar=MapPar, + FigPar=FigPar, + ) else: - - fig,m = get_base1(map_region=map_region, - projection=projection, - coords=coords, - MapPar=MapPar, - FigPar=FigPar, - fig=fig, - ) - + + fig, m = get_base1(map_region=map_region, + projection=projection, + coords=coords, + MapPar=MapPar, + FigPar=FigPar, + fig=fig, + ) + FIGURE.fig = fig FIGURE.m = m FIGURE.ax = fig.gca() @@ -1460,28 +1473,30 @@ def get_FIGURE(fig=None,ax=None,m=None,map_region=None, FIGURE.m = None else: FIGURE.m = m - + if fig == None: FIGURE.fig = plt.figure() - fig=FIGURE.fig + fig = FIGURE.fig else: FIGURE.fig = fig - + if ax == None: FIGURE.ax = fig.gca() ax = FIGURE.ax else: FIGURE.ax = ax - + FIGURE.indices = Structure() FIGURE.indices.texts = len(FIGURE.ax.texts) FIGURE.indices.images = len(FIGURE.ax.images) FIGURE.indices.collections = len(FIGURE.ax.collections) FIGURE.indices.lines = len(FIGURE.ax.lines) - print "Using figure: %s" % FIGURE.fig.number + print + "Using figure: %s" % FIGURE.fig.number return FIGURE + def get_base1(map_region=1, projection=None, figname=None, @@ -1502,42 +1517,46 @@ def get_base1(map_region=1, ## Use map_regions function to define ## input paramters for Basemap - MapPar_sd,FigPar_sd=map_regions(map_region=map_region, - projection=projection, - coords=coords, - MapPar=MapPar) + MapPar_sd, FigPar_sd = map_regions(map_region=map_region, + projection=projection, + coords=coords, + MapPar=MapPar) if MapPar: MapPar_sd.set_with_dict(MapPar) if FigPar: FigPar_sd.set_with_dict(FigPar) ## create the figure. - if fig==None: + if fig == None: axlocs = FigPar_sd.pop('axlocs') - fig=plt.figure(**FigPar_sd) + fig = plt.figure(**FigPar_sd) ax = fig.add_axes(axlocs) else: - ax=fig.gca() - - print Basemap - m=Basemap(**MapPar_sd) - print 'getting base1' - print m + ax = fig.gca() + + print + Basemap + m = Basemap(**MapPar_sd) + print + 'getting base1' + print + m plt.axes(ax) ## make sure axes ax are current ## draw coastlines and political boundaries. m.drawcoastlines(linewidth=0.8) m.drawcountries(linewidth=0.2) m.drawstates(linewidth=0.2) if drawlsmask: - m.drawlsmask(ocean_color='#008EBA',zorder=0) - #m.fillcontinents(zorder=0.) + m.drawlsmask(ocean_color='#008EBA', zorder=0) + # m.fillcontinents(zorder=0.) ## draw parallels and meridians. ## use draw_grid function - MP,MD = draw_grid(m) + MP, MD = draw_grid(m) - if figname!=None: + if figname != None: plt.savefig(figname) - return fig,m + return fig, m + def get_base2(**kwargs): """ Warps NASA Blue Marble Image version @@ -1577,40 +1596,43 @@ def get_base2(**kwargs): # read in jpeg image to rgba array of normalized floats. pilImage = Image.open('land_shallow_topo_2048.jpg') rgba = pil_to_array(pilImage) - rgba = rgba.astype(P.Float32)/255. # convert to normalized floats. + rgba = rgba.astype(P.Float32) / 255. # convert to normalized floats. # define lat/lon grid that image spans (projection='cyl'). - nlons = rgba.shape[1]; nlats = rgba.shape[0] - delta = 360./float(nlons) - lons = P.arange(-180.+0.5*delta,180.,delta) - lats = P.arange(-90.+0.5*delta,90.,delta) + nlons = rgba.shape[1]; + nlats = rgba.shape[0] + delta = 360. / float(nlons) + lons = P.arange(-180. + 0.5 * delta, 180., delta) + lats = P.arange(-90. + 0.5 * delta, 90., delta) # create new figure - fig=P.figure(1,figsize=(8,6)) + fig = P.figure(1, figsize=(8, 6)) # define Lambert Conformal basemap for North America. - mr={'map_region':reg,'coords':coord} - mp,figpar=map_regions(map_region=reg,coords=coord) + mr = {'map_region': reg, 'coords': coord} + mp, figpar = map_regions(map_region=reg, coords=coord) m = Basemap(**mp) - ax = fig.add_axes([0.1,0.1,0.7,0.7]) #need to change back to [0.1,0.1,0.7,0.7] + ax = fig.add_axes([0.1, 0.1, 0.7, 0.7]) # need to change back to [0.1,0.1,0.7,0.7] P.axes(ax) # make the original axes current again # transform to nx x ny regularly spaced native projection grid # nx and ny chosen to have roughly the same horizontal res as original image. - dx = 2.*P.pi*m.rmajor/float(nlons) - nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 - rgba_warped = P.zeros((ny,nx,4),P.Float64) + dx = 2. * P.pi * m.rmajor / float(nlons) + nx = int((m.xmax - m.xmin) / dx) + 1; + ny = int((m.ymax - m.ymin) / dx) + 1 + rgba_warped = P.zeros((ny, nx, 4), P.Float64) # interpolate rgba values from proj='cyl' (geographic coords) to 'lcc' for k in range(4): - rgba_warped[:,:,k] = m.transform_scalar(rgba[:,:,k],lons,lats,nx,ny) + rgba_warped[:, :, k] = m.transform_scalar(rgba[:, :, k], lons, lats, nx, ny) # plot warped rgba image. im = m.imshow(rgba_warped) # draw coastlines. - m.drawcoastlines(linewidth=0.5,color='0.5') + m.drawcoastlines(linewidth=0.5, color='0.5') # draw parallels and meridians. - draw_grid(m,linewidth=0.5, color='0.5') - #draw() - return fig,m + draw_grid(m, linewidth=0.5, color='0.5') + # draw() + return fig, m + def get_base3(**kwargs): """ Custom version for sea ice. """ @@ -1641,30 +1663,30 @@ def get_base3(**kwargs): from matplotlib.toolkits.basemap import Basemap from matplotlib.image import pil_to_array from pylab import title, colorbar, show, close, \ - axes, cm, load, arange, figure, \ - text, savefig, setp, draw, clf, cla + axes, cm, load, arange, figure, \ + text, savefig, setp, draw, clf, cla # read in jpeg image to rgba array of normalized floats. - if basefile!=None: + if basefile != None: pilImage = Image.open(basefile) rgba = pil_to_array(pilImage) - rgba = rgba.astype(np.float32)/255. # convert to normalized floats. + rgba = rgba.astype(np.float32) / 255. # convert to normalized floats. interactive(False) ## create the figure. - if fig==None: - fig=figure(1,\ - dpi=100, - frameon=False, - facecolor=None, - edgecolor=None, - figsize=(12,9.11) - ) - #fig=figure(1,figsize=(57.625,43.75)) - #Use map_regions function to define input paramters for Basemap - mp,figpar=map_regions(map_region) - m=Basemap(**mp) - ax = fig.add_axes([0,0,1,1],frameon=False) #need to change back to [0.1,0.1,0.7,0.7] + if fig == None: + fig = figure(1, \ + dpi=100, + frameon=False, + facecolor=None, + edgecolor=None, + figsize=(12, 9.11) + ) + # fig=figure(1,figsize=(57.625,43.75)) + # Use map_regions function to define input paramters for Basemap + mp, figpar = map_regions(map_region) + m = Basemap(**mp) + ax = fig.add_axes([0, 0, 1, 1], frameon=False) # need to change back to [0.1,0.1,0.7,0.7] axes(ax) # make the original axes current again ## draw coastlines and political boundaries. m.drawcoastlines(linewidth=.5) @@ -1673,13 +1695,12 @@ def get_base3(**kwargs): ## draw parallels and meridians. draw_grid(m, linewidth=0.5) - if figname!=None: + if figname != None: savefig(figname) - return fig,m - + return fig, m -def get_base_image(imagefile,**kwargs): +def get_base_image(imagefile, **kwargs): """ Warps NASA Blue Marble Image version Create basemap figure for plotting on top of returns a figure and a basemap instance. @@ -1723,58 +1744,61 @@ def get_base_image(imagefile,**kwargs): # read in jpeg image to rgba array of normalized floats. pilImage = Image.open(imagefile) rgba = pil_to_array(pilImage) - rgba = rgba.astype(np.float32)/255. # convert to normalized floats. + rgba = rgba.astype(np.float32) / 255. # convert to normalized floats. # define lat/lon grid that image spans (projection='cyl'). - nlons = rgba.shape[1]; nlats = rgba.shape[0] - delta = 360./float(nlons) - lons = np.arange(-180.+0.5*delta,180.,delta) - lats = np.arange(-90.+0.5*delta,90.,delta) + nlons = rgba.shape[1]; + nlats = rgba.shape[0] + delta = 360. / float(nlons) + lons = np.arange(-180. + 0.5 * delta, 180., delta) + lats = np.arange(-90. + 0.5 * delta, 90., delta) # create new figure - fig=plt.figure(1,figsize=(8,6)) + fig = plt.figure(1, figsize=(8, 6)) # define Lambert Conformal basemap for North America. - mr={'map_region':reg,'coords':coord} - mp,figpar=map_regions(map_region=reg,coords=coord) + mr = {'map_region': reg, 'coords': coord} + mp, figpar = map_regions(map_region=reg, coords=coord) m = Basemap(**mp) - ax = fig.add_axes([0.1,0.1,0.7,0.7]) #need to change back to [0.1,0.1,0.7,0.7] + ax = fig.add_axes([0.1, 0.1, 0.7, 0.7]) # need to change back to [0.1,0.1,0.7,0.7] plt.axes(ax) # make the original axes current again # transform to nx x ny regularly spaced native projection grid # nx and ny chosen to have roughly the same horizontal res as original image. - dx = 2.*np.pi*m.rmajor/float(nlons) - nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 - rgba_warped = np.zeros((ny,nx,4),np.float64) + dx = 2. * np.pi * m.rmajor / float(nlons) + nx = int((m.xmax - m.xmin) / dx) + 1; + ny = int((m.ymax - m.ymin) / dx) + 1 + rgba_warped = np.zeros((ny, nx, 4), np.float64) # interpolate rgba values from proj='cyl' (geographic coords) to 'lcc' try: for k in range(4): - rgba_warped[:,:,k] = m.transform_scalar(rgba[:,:,k],lons,lats,nx,ny) + rgba_warped[:, :, k] = m.transform_scalar(rgba[:, :, k], lons, lats, nx, ny) except: rgba_warped = rgba - print 'problem with transform_scalar' + print + 'problem with transform_scalar' # plot warped rgba image. im = m.imshow(rgba_warped) # draw coastlines. - m.drawcoastlines(linewidth=0.5,color='0.5') + m.drawcoastlines(linewidth=0.5, color='0.5') # draw parallels and meridians. - draw_grid(m,linewidth=0.5, color='0.5') - #draw() - return fig,m + draw_grid(m, linewidth=0.5, color='0.5') + # draw() + return fig, m -def set_plotkwargs(kwargs,plot_kwargs=None): +def set_plotkwargs(kwargs, plot_kwargs=None): """ Internal function to check for valid plotting keywords and extracts them into a new kwarg dict for the plotting routine """ - valid_kwargs = ['alpha','animated','antialiased','axes','clip_box',\ - 'clip_on ','clip_path ','color','contains','dash_capstyle',\ - 'dash_joinstyle','dashes','data','drawstyle','figure ',\ - 'label','linestyle','linewidth ','lod','marker', \ - 'markeredgecolor','markeredgewidth ','markerfacecolor ',\ - 'markersize','picker ','snap ','solid_capstyle',\ - 'solid_joinstyle ','transform','url ',\ - 'visible','xdata','ydata','zorder'] + valid_kwargs = ['alpha', 'animated', 'antialiased', 'axes', 'clip_box', \ + 'clip_on ', 'clip_path ', 'color', 'contains', 'dash_capstyle', \ + 'dash_joinstyle', 'dashes', 'data', 'drawstyle', 'figure ', \ + 'label', 'linestyle', 'linewidth ', 'lod', 'marker', \ + 'markeredgecolor', 'markeredgewidth ', 'markerfacecolor ', \ + 'markersize', 'picker ', 'snap ', 'solid_capstyle', \ + 'solid_joinstyle ', 'transform', 'url ', \ + 'visible', 'xdata', 'ydata', 'zorder'] if plot_kwargs is None: plot_kwargs = {} for arg in kwargs.keys(): @@ -1782,9 +1806,10 @@ def set_plotkwargs(kwargs,plot_kwargs=None): plot_kwargs[arg] = kwargs[arg] return plot_kwargs -def plot_track(lon,lat, - figname=None,FIGURE=None,overlay=False,cbar2=False, - zlevel=None,zsize=None,base=1,marker='o', + +def plot_track(lon, lat, + figname=None, FIGURE=None, overlay=False, cbar2=False, + zlevel=None, zsize=None, base=1, marker='o', plotargs=None, map_region=None, projection=None, @@ -1825,12 +1850,11 @@ def plot_track(lon,lat, """ ## make tick lables smaller and set parameters for colorbar - mpl.rcParams['xtick.labelsize']=8.5 - mpl.rcParams['ytick.labelsize']=8.5 + mpl.rcParams['xtick.labelsize'] = 8.5 + mpl.rcParams['ytick.labelsize'] = 8.5 p_cax = mpl.font_manager.FontProperties(size='7') - - #check that lon/lat are iterable + # check that lon/lat are iterable if not np.iterable(lon): lon = [lon] if not np.iterable(lat): @@ -1838,119 +1862,120 @@ def plot_track(lon,lat, scatter = True if zsize == None: - zsize = np.ones(len(lon))*1000 + zsize = np.ones(len(lon)) * 1000 ## Default plot_kwargs plot_kwargs = {} if plotargs != None: - plot_kwargs = set_plotkwargs(plotargs,plot_kwargs) + plot_kwargs = set_plotkwargs(plotargs, plot_kwargs) else: - plotargs = {'zorder':10, - 'alpha':0.35, - 'edgecolor':None, - #'marker':marker #should get rid of marker + plotargs = {'zorder': 10, + 'alpha': 0.35, + 'edgecolor': None, + # 'marker':marker #should get rid of marker } - plot_kwargs = set_plotkwargs(plotargs,plot_kwargs) + plot_kwargs = set_plotkwargs(plotargs, plot_kwargs) - if FIGURE==None: - FIGURE = get_FIGURE(map_region=map_region,projection=projection,coords=coords) + if FIGURE == None: + FIGURE = get_FIGURE(map_region=map_region, projection=projection, coords=coords) ##Get fig if exists - fig=FIGURE.fig - m=FIGURE.m - ax=FIGURE.ax - nullfmt = NullFormatter() - if ax==None: - az=plt.gca() + fig = FIGURE.fig + m = FIGURE.m + ax = FIGURE.ax + nullfmt = NullFormatter() + if ax == None: + az = plt.gca() pos = az.get_position() l, b, w, h = getattr(pos, 'bounds', pos) - ax = fig.add_axes([l,b,w,h],frameon=False) - if overlay==False: + ax = fig.add_axes([l, b, w, h], frameon=False) + if overlay == False: del ax.collections[FIGURE.indices.collections:] # Make sure the axes and figure are current plt.figure(fig.number) plt.axes(ax) - #PRINT CRUISE TRACK - cx,cy = m(lon,lat) - if zlevel!=None: + # PRINT CRUISE TRACK + cx, cy = m(lon, lat) + if zlevel != None: if scatter: ## Set plotting kwargs - #if plotargs == None: + # if plotargs == None: # plotargs = {'zorder':10, # 'alpha':0.35, # 'edgecolor':None # } - plot_kwargs = set_plotkwargs(plotargs,plot_kwargs) - #from pylab import linspace, cm, colorbar,scatter, cla, axes - #m.scatter(cx,cy,25*linspace(0,1,(len(cx))),zlevel,cmap=cm.spectral,marker='o',faceted=False,zorder=10) + plot_kwargs = set_plotkwargs(plotargs, plot_kwargs) + # from pylab import linspace, cm, colorbar,scatter, cla, axes + # m.scatter(cx,cy,25*linspace(0,1,(len(cx))),zlevel,cmap=cm.spectral,marker='o',faceted=False,zorder=10) cmap = plt.get_cmap('jet') - #c=m.scatter(cx,cy,zsize,zlevel,cmap=cmap,marker=marker, + # c=m.scatter(cx,cy,zsize,zlevel,cmap=cmap,marker=marker, # **plot_kwargs) - c=m.scatter(cx,cy,zsize,zlevel,cmap=cmap, - **plot_kwargs) - #pos = az.get_position() - #l, b, w, h = getattr(pos, 'bounds', pos) - #cax = axes([l+w+0.075, b, 0.05, h]) - #c.set_alpha(0.01) - #m.scatter has no color bar, so create a ghost 'scatter' instance + c = m.scatter(cx, cy, zsize, zlevel, cmap=cmap, + **plot_kwargs) + # pos = az.get_position() + # l, b, w, h = getattr(pos, 'bounds', pos) + # cax = axes([l+w+0.075, b, 0.05, h]) + # c.set_alpha(0.01) + # m.scatter has no color bar, so create a ghost 'scatter' instance pos = ax.get_position() l, b, w, h = getattr(pos, 'bounds', pos) jnkfig = plt.figure() - jnkax = jnkfig.add_axes([l,b,w,h],frameon=False) - jnkax.scatter(cx,cy,zsize,zlevel,cmap=cmap,marker=marker, + jnkax = jnkfig.add_axes([l, b, w, h], frameon=False) + jnkax.scatter(cx, cy, zsize, zlevel, cmap=cmap, marker=marker, **plot_kwargs) plt.figure(fig.number) if cbar2 is True: - cax = plt.axes([l+w+0.12, b, 0.02, h-0.035]) + cax = plt.axes([l + w + 0.12, b, 0.02, h - 0.035]) else: - cax = plt.axes([l+w+0.03, b, 0.025, h-0.035]) - #cax = plt.axes([l+w+0.03, b, 0.02, h]) - plt.colorbar(cax=cax) # draw colorbar + cax = plt.axes([l + w + 0.03, b, 0.025, h - 0.035]) + # cax = plt.axes([l+w+0.03, b, 0.02, h]) + plt.colorbar(cax=cax) # draw colorbar if TEX: cax.set_title(r'textit{%s}' % units) else: cax.set_title('%s' % units, fontproperties=p_cax) - - #delete the ghost instance + + # delete the ghost instance plt.close(jnkfig.number) plt.axes(ax) else: - #raise IOError('Vidit is going to add this method') + # raise IOError('Vidit is going to add this method') print('trying now') if np.rank(cx) == 1: points = np.array([cx, cy]).T.reshape(-1, 1, 2) segments = np.concatenate([points[:-1], points[1:]], axis=1) - - lc = LineCollection(segments, cmap=plt.get_cmap('jet'),norm=plt.Normalize(zlevel.min(),zlevel.max())) + + lc = LineCollection(segments, cmap=plt.get_cmap('jet'), norm=plt.Normalize(zlevel.min(), zlevel.max())) lc.set_array(zlevel.flatten()) lc.set_linewidth(3) plt.gca().add_collection(lc) elif np.rank(cx) == 2: - for i in xrange(cx.shape[0]): - points = np.array([cx[i,:],cy[i,:]]).T.reshape(-1, 1, 2) + for i in range(cx.shape[0]): + points = np.array([cx[i, :], cy[i, :]]).T.reshape(-1, 1, 2) segments = np.concatenate([points[:-1], points[1:]], axis=1) - - lc = LineCollection(segments, cmap=plt.get_cmap('jet'),norm=plt.Normalize(zlevel.min(),zlevel.max())) + + lc = LineCollection(segments, cmap=plt.get_cmap('jet'), + norm=plt.Normalize(zlevel.min(), zlevel.max())) lc.set_array(zlevel.flatten()) lc.set_linewidth(3) plt.gca().add_collection(lc) else: - raise InputError("input array shape zlevel cannot be greater than rank-2") - - + raise InputError("input array shape zlevel cannot be greater than rank-2") + + else: if scatter: - m.scatter(cx,cy,**plot_kwargs) + m.scatter(cx, cy, **plot_kwargs) else: - m.plot(cx,cy,**plot_kwargs) - print 'boring' + m.plot(cx, cy, **plot_kwargs) + print('boring') plt.axes(ax) - ax.xaxis.set_major_formatter( nullfmt ) - ax.yaxis.set_major_formatter( nullfmt ) - plt.setp(ax, xticks=[],yticks=[]) - #ax.axesPatch.set_alpha(0.0) + ax.xaxis.set_major_formatter(nullfmt) + ax.yaxis.set_major_formatter(nullfmt) + plt.setp(ax, xticks=[], yticks=[]) + # ax.axesPatch.set_alpha(0.0) if figname: plt.savefig(figname) FIGURE.fig = fig @@ -1958,8 +1983,9 @@ def plot_track(lon,lat, FIGURE.ax = ax return FIGURE -def plot_grid(D,map_region='POLARCAT',dres=0.5, - transform=True,figname=None,fillcontinents=False, + +def plot_grid(D, map_region='POLARCAT', dres=0.5, + transform=True, figname=None, fillcontinents=False, points=False): """ plot an array over a basemap. The required argument "D" is either: @@ -1986,29 +2012,30 @@ def plot_grid(D,map_region='POLARCAT',dres=0.5, """ - print "length of D: %s" % len(D) - + print + "length of D: %s" % len(D) - if isinstance(D,np.ndarray): - #pdb.set_trace() + if isinstance(D, np.ndarray): + # pdb.set_trace() assert len(D.shape) == 2, "D grid must be 2d" - #if len(D)==1: - #assume full earth grid with dres - print 'received grid of shape:', D.shape + # if len(D)==1: + # assume full earth grid with dres + print + 'received grid of shape:', D.shape if D.shape[0] == 720: - lons = np.arange(-180,180,dres) + lons = np.arange(-180, 180, dres) elif D.shape[0] == 721: - lons = np.arange(-180,180.01,dres) + lons = np.arange(-180, 180.01, dres) if D.shape[1] == 360: - lats = np.arange(-90,90,dres) - lons = np.arange(-180,180.01,dres) + lats = np.arange(-90, 90, dres) + lons = np.arange(-180, 180.01, dres) elif D.shape[1] == 361: - lats = np.arange(-90,90.01,dres) - lons = np.arange(-180,180.01,dres) + lats = np.arange(-90, 90.01, dres) + lons = np.arange(-180, 180.01, dres) points = False z = D.T - elif len(D)==3: + elif len(D) == 3: points = True x = D[0] y = D[1] @@ -2019,115 +2046,116 @@ def plot_grid(D,map_region='POLARCAT',dres=0.5, lons = x lats = y - ## CHANGED THIS HERE FOR THE CODEX ## -## Be sure to set map_region=None ## - if isinstance(map_region,str): - print "getting basemap with map_region: %s" % map_region - fig,m = get_base1(map_region=map_region) + ## CHANGED THIS HERE FOR THE CODEX ## + ## Be sure to set map_region=None ## + if isinstance(map_region, str): + print + "getting basemap with map_region: %s" % map_region + fig, m = get_base1(map_region=map_region) else: fig = plt.figure() ax = fig.add_axes() a = 6378.273e3 ec = 0.081816153 - b = a*np.sqrt(1.-ec**2) - m = Basemap(projection='stere',lat_0=90,lon_0=-45,lat_ts=70,\ - llcrnrlat=33.92,llcrnrlon=279.96,\ - urcrnrlon=102.34,urcrnrlat=31.37,\ - rsphere=(a,b)) -# m = Basemap(width=12000000,height=8000000, -# resolution='l',projection='npstere',\ -# lat_ts=50,lat_0=50,lon_0=-107.)#Set up a basemap + b = a * np.sqrt(1. - ec ** 2) + m = Basemap(projection='stere', lat_0=90, lon_0=-45, lat_ts=70, \ + llcrnrlat=33.92, llcrnrlon=279.96, \ + urcrnrlon=102.34, urcrnrlat=31.37, \ + rsphere=(a, b)) + # m = Basemap(width=12000000,height=8000000, + # resolution='l',projection='npstere',\ + # lat_ts=50,lat_0=50,lon_0=-107.)#Set up a basemap m.drawcoastlines() if points == True: - # Plot individual data points - norm = colors.normalize(z.min(),z.max()) - xpt,ypt = m(x,y) + # Plot individual data points + norm = colors.normalize(z.min(), z.max()) + xpt, ypt = m(x, y) cmap = [cm.jet(norm(i)) for i in z] - m.scatter(xpt,ypt,s=z/100.,edgecolors='none',color=cmap) -# for i in range(len(y)): -# xpt,ypt = m(x[i],y[i]) -# cmap = cm.jet(norm(z[i])) -# #cmap = 'k' -# m.plot([xpt],[ypt],'.',color=cmap,markersize=2) -# #plt.plot(x[i],y[i],'.',color=cmap,markersize=2) + m.scatter(xpt, ypt, s=z / 100., edgecolors='none', color=cmap) + # for i in range(len(y)): + # xpt,ypt = m(x[i],y[i]) + # cmap = cm.jet(norm(z[i])) + # #cmap = 'k' + # m.plot([xpt],[ypt],'.',color=cmap,markersize=2) + # #plt.plot(x[i],y[i],'.',color=cmap,markersize=2) if points == False: - #transform Z data into projection - #transform to nx x ny regularly spaced native projection grid - dx = 2.*np.pi*m.rmajor/len(lons) - nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 + # transform Z data into projection + # transform to nx x ny regularly spaced native projection grid + dx = 2. * np.pi * m.rmajor / len(lons) + nx = int((m.xmax - m.xmin) / dx) + 1; + ny = int((m.ymax - m.ymin) / dx) + 1 if transform: # Need this if we choose lon,lat approach - Zt,xx,yy = m.transform_scalar(z,lons,lats,nx,ny,returnxy=True) + Zt, xx, yy = m.transform_scalar(z, lons, lats, nx, ny, returnxy=True) else: Zt = z - print m.projection + print + m.projection if 'moll' in m.projection: - x,y = m(xx,yy) - m.contourf(x,y,Zt) + x, y = m(xx, yy) + m.contourf(x, y, Zt) else: m.imshow(Zt) if fillcontinents: m.fillcontinents() plt.colorbar() - #plt.imshow(Z) + # plt.imshow(Z) if figname != None: - #plt.ylim([40,90]) - #plt.title('data locations on mercator grid') + # plt.ylim([40,90]) + # plt.title('data locations on mercator grid') plt.savefig(figname) else: plt.show() - return fig,m + return fig, m -def plot_imshow(x,y,z,\ - data_range=None,\ - units = 'ns m^2 / kg',\ - datainfo_str = ' ',\ - rel_i = 0, plottitle = None,\ - globe = False, lon_0 = -110, lat_0 = 45,\ +def plot_imshow(x, y, z, \ + data_range=None, \ + units='ns m^2 / kg', \ + datainfo_str=' ', \ + rel_i=0, plottitle=None, \ + globe=False, lon_0=-110, lat_0=45, \ map_region=None, projection=None, dropm=None, coords=None, overlay=False, transform=False, FIGURE=None): - """ plot x,y,z values output using matplotlib basemap and the imshow function. """ if FIGURE == None: - FIGURE = get_FIGURE(map_region=map_region,projection=projection,coords=coords) - + FIGURE = get_FIGURE(map_region=map_region, projection=projection, coords=coords) if dropm != None: try: del m plt.close('all') except: - print 'could not drop m' + print('could not drop m') - ## Extract grid and release info from header - ## (Try block for different versions) - #try: - #header2vars = ['outlon0','outlat0','dxout','dyout','numxgrid','numygrid','xpoint','ypoint'] - #for k in header2vars: - #key2var(H,k) - #except: - #header2vars = ['outlon0','outlat0','dxout','dyout','numxgrid','numygrid','xp1','yp1'] - #for k in header2vars: - #key2var(H,k) + ## Extract grid and release info from header + ## (Try block for different versions) + # try: + # header2vars = ['outlon0','outlat0','dxout','dyout','numxgrid','numygrid','xpoint','ypoint'] + # for k in header2vars: + # key2var(H,k) + # except: + # header2vars = ['outlon0','outlat0','dxout','dyout','numxgrid','numygrid','xp1','yp1'] + # for k in header2vars: + # key2var(H,k) - #releaselocation = (xpoint[rel_i],ypoint[rel_i]) + # releaselocation = (xpoint[rel_i],ypoint[rel_i]) ## make tick lables smaller - mpl.rcParams['xtick.labelsize']=6 - mpl.rcParams['ytick.labelsize']=6 + mpl.rcParams['xtick.labelsize'] = 6 + mpl.rcParams['ytick.labelsize'] = 6 fig = FIGURE.fig m = FIGURE.m @@ -2138,72 +2166,73 @@ def plot_imshow(x,y,z,\ plt.axes(ax) ### set up transformations for the data array - #if m.projection not in ['cyl','merc','mill']: - #lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ), dyout )[:-1] - #lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ), dxout )[:-1] - #data = data[:-1,:-1] - #else: - #lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ), dyout ) - #lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ), dxout ) - lats=y - lons=x - data=z + # if m.projection not in ['cyl','merc','mill']: + # lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ), dyout )[:-1] + # lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ), dxout )[:-1] + # data = data[:-1,:-1] + # else: + # lats = np.arange( outlat0, ( outlat0 + ( numygrid*dyout ) ), dyout ) + # lons = np.arange( outlon0, ( outlon0 + ( numxgrid*dxout ) ), dxout ) + lats = y + lons = x + data = z ## transform to nx x ny regularly spaced native projection grid if transform: - dx = 2.*np.pi*m.rmajor/len(lons) - nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 - topodat = m.transform_scalar(data,lons,lats,nx,ny) + dx = 2. * np.pi * m.rmajor / len(lons) + nx = int((m.xmax - m.xmin) / dx) + 1; + ny = int((m.ymax - m.ymin) / dx) + 1 + topodat = m.transform_scalar(data, lons, lats, nx, ny) else: topodat = data - ### get min/max range - #if data_range != None: - #dat_min = data_range[0] - #dat_max = data_range[1] - #else: - #dat_min = data.min() - #dat_max = data.max() + ### get min/max range + # if data_range != None: + # dat_min = data_range[0] + # dat_max = data_range[1] + # else: + # dat_min = data.min() + # dat_max = data.max() ### create logorithmic color scale ### method uses strings to get the magnitude of variable - #ss = "%1.2e" % (dat_max) - #dmx = int('%s%s' % (ss[-3],int(ss[-2:]))) + # ss = "%1.2e" % (dat_max) + # dmx = int('%s%s' % (ss[-3],int(ss[-2:]))) ##dmx=float(len(str(int(np.ceil(dat_max))))) - #ss = "%1.2e" % (.0000001*(dat_max-dat_min)) - #dmn = int('%s%s' % (ss[-3],int(ss[-2:]))) + # ss = "%1.2e" % (.0000001*(dat_max-dat_min)) + # dmn = int('%s%s' % (ss[-3],int(ss[-2:]))) ### create equally spaced range - #logspace = 10.**np.linspace(dmn, dmx, 100) - #clevs = [i for i in logspace] + # logspace = 10.**np.linspace(dmn, dmx, 100) + # clevs = [i for i in logspace] ## Get the current axes, and properties for use later pos = ax.get_position() l, b, w, h = pos.bounds ## draw land sea mask - #m.fillcontinents(zorder=0) - m.drawlsmask(ocean_color='white',zorder=0) + # m.fillcontinents(zorder=0) + m.drawlsmask(ocean_color='white', zorder=0) ### Plot Release Location - #xpt,ypt = m(releaselocation[0],releaselocation[1]) + # xpt,ypt = m(releaselocation[0],releaselocation[1]) ### Remove prior location point - #try: - #del ax.lines[-1] - #except: - #pass - #location, = m.plot([xpt],[ypt],'bx',linewidth=6,markersize=20,zorder=1000) + # try: + # del ax.lines[-1] + # except: + # pass + # location, = m.plot([xpt],[ypt],'bx',linewidth=6,markersize=20,zorder=1000) ### Plot the footprint - #if overlay == False: - #try: - #del ax.images[0] - #except: - #print 'could not delete image' + # if overlay == False: + # try: + # del ax.images[0] + # except: + # print('could not delete image') ## Set up the IMAGE ## cmapnames = ['jet', 'hsv', 'gist_ncar', 'gist_rainbow', 'cool', 'spectral'] colmap = plt.get_cmap('gist_ncar_r') - im = m.imshow(topodat,cmap=colmap) + im = m.imshow(topodat, cmap=colmap) ## OLDER METHODS ## im = m.imshow(topodat,cmap=cm.s3pcpn,vmin=clevs[0],vmax=clevs[-1]) ## im = m.contourf(nx,ny,topodat, locator=ticker.LogLocator()) @@ -2212,23 +2241,23 @@ def plot_imshow(x,y,z,\ ## make a copy of the image object, change ## colormap to linear version of the precip colormap. im2 = copy.copy(im) - #im2.set_cmap(cm.s3pcpn_l) + # im2.set_cmap(cm.s3pcpn_l) im2.set_cmap(colmap) ## create new axis for colorbar. - cax = plt.axes([l+w+0.03, b, 0.025, h-0.035]) + cax = plt.axes([l + w + 0.03, b, 0.025, h - 0.035]) ## using im2, not im (hack to prevent colors from being ## too compressed at the low end on the colorbar - results ## from highly nonuniform colormap) - cb = plt.colorbar(im2, cax)#, format='%3.2g') # draw colorbar + cb = plt.colorbar(im2, cax) # , format='%3.2g') # draw colorbar ## set colorbar label and ticks - #p_cax = mpl.font_manager.FontProperties(size='6') - #clabels = clevs[::10] ##clevs, by 10 steps - #clabels.append(clevs[-1]) ## add the last label - #cb.ax.set_yticks(np.linspace(0,1,len(clabels))) - #cb.ax.set_yticklabels(['%3.2g' % cl for cl in clabels], - #fontproperties=p_cax) - #cax.set_title('sensitivity\n(%s)' % units, - #fontproperties=p_cax) + # p_cax = mpl.font_manager.FontProperties(size='6') + # clabels = clevs[::10] ##clevs, by 10 steps + # clabels.append(clevs[-1]) ## add the last label + # cb.ax.set_yticks(np.linspace(0,1,len(clabels))) + # cb.ax.set_yticklabels(['%3.2g' % cl for cl in clabels], + # fontproperties=p_cax) + # cax.set_title('sensitivity\n(%s)' % units, + # fontproperties=p_cax) ## make the original axes current again plt.axes(ax) @@ -2246,14 +2275,14 @@ def plot_imshow(x,y,z,\ except: pass - plt.text(l,b+1000, + plt.text(l, b + 1000, datainfo_str, - fontsize = 10, - bbox = dict(boxstyle="round", - ec=(1., 0.5, 0.5), - fc=(1., 0.8, 0.8), - alpha=0.8 - ) + fontsize=10, + bbox=dict(boxstyle="round", + ec=(1., 0.5, 0.5), + fc=(1., 0.8, 0.8), + alpha=0.8 + ) ) FIGURE.ax = ax @@ -2262,22 +2291,21 @@ def plot_imshow(x,y,z,\ if plottitle != None: plt.title(plottitle) - #plt = plt + # plt = plt return FIGURE #### LateX ##### if TEX: - mpl.rc('font', **{'family':'sans-serif', - 'sans-serif':['Helvetica'] + mpl.rc('font', **{'family': 'sans-serif', + 'sans-serif': ['Helvetica'] } ) - mpl.rc('text',usetex=True) + mpl.rc('text', usetex=True) - -#### SOME SPECIFIC FUNCTIONS FOR PLOTTING SATELLITE GROUND TRACKS #### +#### SOME SPECIFIC FUNCTIONS FOR PLOTTING SATELLITE GROUND TRACKS #### def plot_sattracks(**kwargs): """ Plot satellite tracks. @@ -2331,7 +2359,7 @@ def plot_sattracks(**kwargs): overlay = False if 'plotargs' in kwargs.keys(): plotargs = kwargs['plotargs'] - plot_kwargs = set_plotkwargs(plotargs,plot_kwargs) + plot_kwargs = set_plotkwargs(plotargs, plot_kwargs) else: plotargs = None if 'print_times' in kwargs.keys(): @@ -2347,111 +2375,111 @@ def plot_sattracks(**kwargs): import mapping as mp import datetime from matplotlib.ticker import NullFormatter - from pylab import figure, axes, setp, title, text, savefig, show, close, gca - - #DEBUG - for k,v in kwargs.iteritems(): - print 'for %s ==> %s' % (k,v) - - #A few plot definitions - if FIGURE==None: - FIGURE = get_FIGURE(map_region=map_region,projection=projection,coords=coords) - - fig=FIGURE.fig - m=FIGURE.m - az=FIGURE.ax - nullfmt = NullFormatter() - -# if fig==None: -# fig,m=mp.get_base3(tstart=tstart,tstop=tstop,map_region=map_region) -# ax=gca() -# pos = ax.get_position() -# l, b, w, h = getattr(pos, 'bounds', pos) -# az = fig.add_axes([l,b,w,h],frameon=False) -# if az!=None and overlay==0: -# az.cla() + from pylab import figure, axes, setp, title, text, savefig, show, close, gca + + # DEBUG + for k, v in kwargs.iteritems(): + print('for %s ==> %s' % (k, v)) + + # A few plot definitions + if FIGURE == None: + FIGURE = get_FIGURE(map_region=map_region, projection=projection, coords=coords) + + fig = FIGURE.fig + m = FIGURE.m + az = FIGURE.ax + nullfmt = NullFormatter() + + # if fig==None: + # fig,m=mp.get_base3(tstart=tstart,tstop=tstop,map_region=map_region) + # ax=gca() + # pos = ax.get_position() + # l, b, w, h = getattr(pos, 'bounds', pos) + # az = fig.add_axes([l,b,w,h],frameon=False) + # if az!=None and overlay==0: + # az.cla() figure(fig.number) - #Get/Print Satellite Track - if overlay==False: + # Get/Print Satellite Track + if overlay == False: del az.lines[FIGURE.indices.lines:] del az.texts[FIGURE.indices.texts:] if satellite != None: - m,doy=mp.Sat_tracks(m,start_time=tstart,stop_time=tstop,satellite=satellite, - ifile=ifile,plotargs=plotargs,print_times=print_times,Sets=Sets) + m, doy = mp.Sat_tracks(m, start_time=tstart, stop_time=tstop, satellite=satellite, + ifile=ifile, plotargs=plotargs, print_times=print_times, Sets=Sets) else: - m,doy = None,None + m, doy = None, None - az.xaxis.set_major_formatter( nullfmt ) - az.yaxis.set_major_formatter( nullfmt ) - setp(az, xticks=[],yticks=[]) - #az.axesPatch.set_alpha(0.0) + az.xaxis.set_major_formatter(nullfmt) + az.yaxis.set_major_formatter(nullfmt) + setp(az, xticks=[], yticks=[]) + # az.axesPatch.set_alpha(0.0) # set title. - #pdb.set_trace() + # pdb.set_trace() tstart = doy[0] tstop = doy[1] if title != None: if tstop: - title('%s @ %s through %s ' % (satellite,tstart.strftime('%d-%b-%Y'), + title('%s @ %s through %s ' % (satellite, tstart.strftime('%d-%b-%Y'), tstop.strftime('%d-%b-%Y'))) else: - title('%s @ %s ' % (satellite,tstart.strftime('%d-%b-%Y')) ) + title('%s @ %s ' % (satellite, tstart.strftime('%d-%b-%Y'))) else: pass - #title('%s' % satellite) - if figname!=None: + # title('%s' % satellite) + if figname != None: savefig(figname) - #show() - #FIGURE=[fig,m,az] + # show() + # FIGURE=[fig,m,az] FIGURE.fig = fig FIGURE.m = m FIGURE.ax = az return FIGURE -def read_satdata(c,dat,D,start_time,stop_time=None,satellite='Calipso'): +def read_satdata(c, dat, D, start_time, stop_time=None, satellite='Calipso'): """ Reads text files for a few defined format satellite groundtrack data .. note:: Not ready """ import datetime - if stop_time==None: - stop_time=start_time+datetime.timedelta(1) - if satellite=='Calipso': + if stop_time == None: + stop_time = start_time + datetime.timedelta(1) + if satellite == 'Calipso': for l in c.readlines(): - l=l.strip().split() - if len(l)==7: + l = l.strip().split() + if len(l) == 7: dat.append(l) for d in dat: - t=' '.join(d[:4]) - t=t[:-4] - t=datetime.datetime.strptime(t,'%d %b %Y %H:%M:%S') - if t>start_time and t start_time and t < stop_time: + lat, lon, rng = d[4:] + D.append((t, lon, lat, rng)) + # pdb.set_trace() + if satellite in ['Terra', 'Aqua']: for l in c.readlines(): - l=l.strip().split() - if l!=[]: - if l[0]!='GMT': - if len(l)==3: - y,m,d=l[0].split('/') - if len(l)>3: - hr,mn,sc=l[0].split(':') - t=datetime.datetime(int(y),int(m),int(d),int(hr),int(mn),int(sc)) - #print start_time, t, stop_time - if t>start_time and t 3: + hr, mn, sc = l[0].split(':') + t = datetime.datetime(int(y), int(m), int(d), int(hr), int(mn), int(sc)) + # print(start_time, t, stop_time) + if t > start_time and t < stop_time: + # print('yes!') try: - lat,lon,rng=l[1:] - D.append((t,lon,lat,rng)) + lat, lon, rng = l[1:] + D.append((t, lon, lat, rng)) except: - lat,lon,hdg,ltlat,ltlon,rtlat,rtlon = l[1:] - D.append((t,lon,lat,hdg)) + lat, lon, hdg, ltlat, ltlon, rtlat, rtlon = l[1:] + D.append((t, lon, lat, hdg)) - t0 = datetime.datetime(2009,04,15) + t0 = datetime.datetime(2009, 4, 15) Sets = {} i = 0 @@ -2459,12 +2487,13 @@ def read_satdata(c,dat,D,start_time,stop_time=None,satellite='Calipso'): if t[0] - t0 < datetime.timedelta(minutes=60): Sets[i].append(t) else: - t0=t[0] + t0 = t[0] i += 1 - Sets[i]=[t] + Sets[i] = [t] return Sets + def read_LARC_predict(filepath): """ read path prediction file from the LARC website """ lines = open(filepath).readlines() @@ -2472,18 +2501,18 @@ def read_LARC_predict(filepath): for l in lines: s = l.strip().split() if len(s) == 3: - day = dt.datetime.strptime(s[0],'%Y/%m/%d') + day = dt.datetime.strptime(s[0], '%Y/%m/%d') continue if s: if 'GMT' in s[0]: pass else: - H,M,S =[int(i) for i in s[0].split(':')] - time = dt.datetime(day.year,day.month,day.day,H,M,S) + H, M, S = [int(i) for i in s[0].split(':')] + time = dt.datetime(day.year, day.month, day.day, H, M, S) lat = float(s[1]) lon = -1 * float(s[2]) hdg = float(s[3]) - out.append([time,lon,lat,hdg]) + out.append([time, lon, lat, hdg]) out = np.array(out) Sets = {} @@ -2497,17 +2526,18 @@ def read_LARC_predict(filepath): if t[0] - t0 < dt.timedelta(minutes=60): Sets[i].append(t) else: - t0=t[0] + t0 = t[0] i += 1 - Sets[i]=[t] + Sets[i] = [t] return Sets -def Sat_tracks(map,start_time=None,stop_time=None, - satellite='Calipso',ifile=None, + +def Sat_tracks(map, start_time=None, stop_time=None, + satellite='Calipso', ifile=None, plotargs=None, - get_Sets=False,Sets=None, - print_times=True,plot_frq=30): + get_Sets=False, Sets=None, + print_times=True, plot_frq=30): """ Plots CALIPSO/Terra Ground tracks. .. note:: @@ -2522,77 +2552,81 @@ def Sat_tracks(map,start_time=None,stop_time=None, import datetime, os, sys from pylab import text, title, floor, savefig from matplotlib import colors, cm - if ifile==None: - if satellite=='Calipso': + if ifile == None: + if satellite == 'Calipso': if 'lin' in sys.platform: - #SATELLITEDAT='/mnt/win/07_jfb/jfbin/pybin/data/CALIPSO_SpecialProduct_NorthAtlantic_20080408.txt' - SATELLITEDAT='/xnilu_wrk/jfb/DATASETS_TEMPLATES/CALIPSO/CALIPSO_SpecialProduct_PAM-ARCMIP_20110324.txt' - else: SATELLITEDAT='C:\07_jfb\jfbin\pybin\data\CALIPSO_SpecialProduct_NorthAtlantic_20080408.txt' - elif satellite=='Terra': + # SATELLITEDAT='/mnt/win/07_jfb/jfbin/pybin/data/CALIPSO_SpecialProduct_NorthAtlantic_20080408.txt' + SATELLITEDAT = '/xnilu_wrk/jfb/DATASETS_TEMPLATES/CALIPSO/CALIPSO_SpecialProduct_PAM-ARCMIP_20110324.txt' + else: + SATELLITEDAT = 'C:\07_jfb\jfbin\pybin\data\CALIPSO_SpecialProduct_NorthAtlantic_20080408.txt' + elif satellite == 'Terra': if 'lin' in sys.platform: - SATELLITEDAT='/mnt/win/07_jfb/jfbin/pybin/data/terra.dat' - else: SATELLITEDAT='C:\07_jfb\jfbin\pybin\data\terra.dat' + SATELLITEDAT = '/mnt/win/07_jfb/jfbin/pybin/data/terra.dat' + else: + SATELLITEDAT = 'C:\07_jfb\jfbin\pybin\data\terra.dat' else: - SATELLITEDAT=None; print 'Need satellite data file!!!'; AttributeError + SATELLITEDAT = None + print('Need satellite data file!!!') + AttributeError else: - SATELLITEDAT=ifile - print 'Using: '+ifile + SATELLITEDAT = ifile + print('Using: ' + ifile) if plotargs == None: - plot_kwargs = {'linewidth':1, + plot_kwargs = {'linewidth': 1, 'color': 'b', 'linestyle': '--'} else: plot_kwargs = set_plotkwargs(plotargs) - if start_time==None: - start_time=datetime.datetime(2010,05,11) - if stop_time==None: - stop_time=start_time+datetime.timedelta(30) + if start_time == None: + start_time = datetime.datetime(2010, 5, 11) + if stop_time == None: + stop_time = start_time + datetime.timedelta(30) if Sets == None: - c=file(SATELLITEDAT,'r') - dat,D=[],[] - #doy=start_time.strftime('%Y%m%d') - Sets=read_satdata(c,dat,D,start_time,stop_time,satellite) + c = file(SATELLITEDAT, 'r') + dat, D = [], [] + # doy=start_time.strftime('%Y%m%d') + Sets = read_satdata(c, dat, D, start_time, stop_time, satellite) if get_Sets: return Sets - #PRINT TRACK + # PRINT TRACK tinit = Sets[1][0][0] tend = Sets[24][0][0] - norm = colors.normalize(0,len(Sets)) - for k,D in enumerate(Sets.values()): + norm = colors.normalize(0, len(Sets)) + for k, D in enumerate(Sets.values()): cmap = cm.jet(norm(k)) - clon=[-1* float(i[1]) for i in D] - clat=[float(i[2]) for i in D] - cx,cy = map(clon,clat) - plot_kwargs.update({'color':cmap}) - map.plot(cx,cy,**plot_kwargs) - #Print legends + clon = [-1 * float(i[1]) for i in D] + clat = [float(i[2]) for i in D] + cx, cy = map(clon, clat) + plot_kwargs.update({'color': cmap}) + map.plot(cx, cy, **plot_kwargs) + # Print legends if 1: label = D[0][0].strftime('%Y%m%d') - #Print TRACK TIMES + # Print TRACK TIMES if print_times: - frq=plot_frq - tlon=[-1* float(i[1]) for i in D[::frq]]; - tlat=[float(i[2]) for i in D[::frq]]; - tx,ty=map(tlon,tlat) - map.plot(tx,ty,'r+',label=label,**plot_kwargs) - ty=[i+1000 for i in ty] #reposition for text placement - tx=[i+1000 for i in tx] - t_txt=[i[0].strftime('%H:%M') for i in D[::frq]] + frq = plot_frq + tlon = [-1 * float(i[1]) for i in D[::frq]] + tlat = [float(i[2]) for i in D[::frq]] + tx, ty = map(tlon, tlat) + map.plot(tx, ty, 'r+', label=label, **plot_kwargs) + ty = [i + 1000 for i in ty] # reposition for text placement + tx = [i + 1000 for i in tx] + t_txt = [i[0].strftime('%H:%M') for i in D[::frq]] for i in range(len(tx)): - if tx[i]map.llcrnrx and ty[i]>map.llcrnry: - text(tx[i],ty[i],t_txt[i],size=8) - #title('%s Track: %s' % (satellite,doy)) - - return map,[tinit,tend] + if tx[i] < map.urcrnrx and ty[i] < map.urcrnry and tx[i] > map.llcrnrx and ty[i] > map.llcrnry: + text(tx[i], ty[i], t_txt[i], size=8) + # title('%s Track: %s' % (satellite,doy)) + return map, [tinit, tend] -def greatCircleDistance(RLAT1,RLON1,RLAT2,RLON2): + +def greatCircleDistance(RLAT1, RLON1, RLAT2, RLON2): """ C C PROGRAM HISTORY LOG: @@ -2611,60 +2645,60 @@ def greatCircleDistance(RLAT1,RLON1,RLAT2,RLON2): C C """ - RERTH=6.3712E6 - PI=3.14159265358979 - DPR=180./PI - if abs(RLAT1-RLAT2) < 0.03 and abs(RLON1-RLON2) < 0.03: - DISTANCE=0. + RERTH = 6.3712E6 + PI = 3.14159265358979 + DPR = 180. / PI + if abs(RLAT1 - RLAT2) < 0.03 and abs(RLON1 - RLON2) < 0.03: + DISTANCE = 0. else: - CLAT1=math.cos(RLAT1/DPR) - SLAT1=math.sin(RLAT1/DPR) - CLAT2=math.cos(RLAT2/DPR) - SLAT2=math.sin(RLAT2/DPR) - CDLON=math.cos((RLON1-RLON2)/DPR) - CRD=SLAT1*SLAT2+CLAT1*CLAT2*CDLON - DISTANCE=RERTH*math.acos(CRD)/1000. + CLAT1 = math.cos(RLAT1 / DPR) + SLAT1 = math.sin(RLAT1 / DPR) + CLAT2 = math.cos(RLAT2 / DPR) + SLAT2 = math.sin(RLAT2 / DPR) + CDLON = math.cos((RLON1 - RLON2) / DPR) + CRD = SLAT1 * SLAT2 + CLAT1 * CLAT2 * CDLON + DISTANCE = RERTH * math.acos(CRD) / 1000. return DISTANCE + gcd = greatCircleDistance -def gridarea(xl,yl,xr,yr,method='gcd'): + +def gridarea(xl, yl, xr, yr, method='gcd'): """ calculates grid area as trapezoid using great circle distance """ - #yres = (yr-yl) - #lat = yl + yres - #R = 6371.009 - #if method == 'badc': - #A = (R**2)*(np.radians(xr)-np.radians(xl))*(np.sin(yr)-np.sin(yl)) - #return A - #if method == 'drmath': - #A = ((np.pi/180)*R**2)*abs((np.sin(yr)-np.sin(yl)))*abs(xr-xl) - #return A - #if method == 'radians': - ##Area is calculated using the latitude at the upper bound of the grid-cell, using Radians - ##lat = latitude of center of the grid cell. - #radians = (90.0 - (lat+yres/2))*3.141593/180.0 - ##Calculate cosines: - #cosines = np.cos(radians)-np.cos(radians + (yres*3.141593)/180.0) - ##Calculate area in square kilometers: - #area = (6371221.3*6371221.3*3.141593*cosines/360.0)*1.0e-6 - #return area + # yres = (yr-yl) + # lat = yl + yres + # R = 6371.009 + # if method == 'badc': + # A = (R**2)*(np.radians(xr)-np.radians(xl))*(np.sin(yr)-np.sin(yl)) + # return A + # if method == 'drmath': + # A = ((np.pi/180)*R**2)*abs((np.sin(yr)-np.sin(yl)))*abs(xr-xl) + # return A + # if method == 'radians': + ##Area is calculated using the latitude at the upper bound of the grid-cell, using Radians + ##lat = latitude of center of the grid cell. + # radians = (90.0 - (lat+yres/2))*3.141593/180.0 + ##Calculate cosines: + # cosines = np.cos(radians)-np.cos(radians + (yres*3.141593)/180.0) + ##Calculate area in square kilometers: + # area = (6371221.3*6371221.3*3.141593*cosines/360.0)*1.0e-6 + # return area if method == 'gcd': - dxa = gcd(yr,xl,yr,xr) - dxb = gcd(yl,xl,yl,xr) - dy = gcd(yl,xl,yr,xl) - #areafract = (dlon*dlat)/((xr-xl)/(yr-yl)) - #print dxa,dxb,dy - area = 0.5*dy*(dxa+dxb) + dxa = gcd(yr, xl, yr, xr) + dxb = gcd(yl, xl, yl, xr) + dy = gcd(yl, xl, yr, xl) + # areafract = (dlon*dlat)/((xr-xl)/(yr-yl)) + # print(dxa, dxb, dy) + area = 0.5 * dy * (dxa + dxb) return area - - __version__ = '1.0.1' -class GreatCircle(object): +class GreatCircle(object): """ formula for perfect sphere from Ed Williams' 'Aviation Formulary' (http://williams.best.vwh.net/avform.htm) @@ -2674,8 +2708,7 @@ class GreatCircle(object): Contact: Jeff Whitaker """ - - def __init__(self,rmajor,rminor,lon1,lat1,lon2,lat2): + def __init__(self, rmajor, rminor, lon1, lat1, lon2, lat2): """ Define a great circle by specifying: rmajor - radius of major axis of ellipsoid @@ -2690,40 +2723,39 @@ def __init__(self,rmajor,rminor,lon1,lat1,lon2,lat2): distance - distance along great circle in radians. lon1,lat1,lon2,lat2 - start and end points (in radians). """ - if rmajor==0 and rminor==0: - + if rmajor == 0 and rminor == 0: # WGS84 a = 6378137.0 b = 6356752.3142 - f = (a-b)/a - rmajor = (2*a+b)/3. - rminor = (2*a+b)/3. + f = (a - b) / a + rmajor = (2 * a + b) / 3. + rminor = (2 * a + b) / 3. # convert to radians from degrees. lat1 = math.radians(lat1) lon1 = math.radians(lon1) lat2 = math.radians(lat2) lon2 = math.radians(lon2) self.a = rmajor - self.f = (rmajor-rminor)/rmajor + self.f = (rmajor - rminor) / rmajor self.lat1 = lat1 self.lat2 = lat2 self.lon1 = lon1 self.lon2 = lon2 # distance along geodesic in meters. - d,a12,a21 = vinc_dist(self.f, self.a, lat1, lon1, lat2, lon2 ) + d, a12, a21 = vinc_dist(self.f, self.a, lat1, lon1, lat2, lon2) self.distance = d self.azimuth12 = a12 self.azimuth21 = a21 # great circle arc-length distance (in radians). - self.gcarclen = 2.*math.asin(math.sqrt((math.sin((lat1-lat2)/2))**2+\ - math.cos(lat1)*math.cos(lat2)*(math.sin((lon1-lon2)/2))**2)) + self.gcarclen = 2. * math.asin(math.sqrt((math.sin((lat1 - lat2) / 2)) ** 2 + \ + math.cos(lat1) * math.cos(lat2) * (math.sin((lon1 - lon2) / 2)) ** 2)) # check to see if points are antipodal (if so, route is undefined). if self.gcarclen == math.pi: self.antipodal = True else: self.antipodal = False - def points(self,npoints): + def points(self, npoints): """ compute arrays of npoints equally spaced intermediate points along the great circle. @@ -2739,32 +2771,33 @@ def points(self,npoints): """ # must ask for at least 2 points. if npoints <= 1: - raise ValueError,'npoints must be greater than 1' + raise ValueError('npoints must be greater than 1') elif npoints == 2: - return [math.degrees(self.lon1),math.degrees(self.lon2)],[math.degrees(self.lat1),math.degrees(self.lat2)] + return [math.degrees(self.lon1), math.degrees(self.lon2)], [math.degrees(self.lat1), + math.degrees(self.lat2)] # can't do it if endpoints are antipodal, since # route is undefined. if self.antipodal: - raise ValueError,'cannot compute intermediate points on a great circle whose endpoints are antipodal' + raise ValueError('cannot compute intermediate points on a great circle whose endpoints are antipodal') d = self.gcarclen - delta = 1.0/(npoints-1) - f = delta*np.arange(npoints) # f=0 is point 1, f=1 is point 2. - incdist = self.distance/(npoints-1) + delta = 1.0 / (npoints - 1) + f = delta * np.arange(npoints) # f=0 is point 1, f=1 is point 2. + incdist = self.distance / (npoints - 1) lat1 = self.lat1 lat2 = self.lat2 lon1 = self.lon1 lon2 = self.lon2 # perfect sphere, use great circle formula if self.f == 0.: - A = np.sin((1-f)*d)/math.sin(d) - B = np.sin(f*d)/math.sin(d) - x = A*math.cos(lat1)*math.cos(lon1)+B*math.cos(lat2)*math.cos(lon2) - y = A*math.cos(lat1)*math.sin(lon1)+B*math.cos(lat2)*math.sin(lon2) - z = A*math.sin(lat1) +B*math.sin(lat2) - lats=np.arctan2(z,np.sqrt(x**2+y**2)) - lons=np.arctan2(y,x) - lons = map(math.degrees,lons.tolist()) - lats = map(math.degrees,lats.tolist()) + A = np.sin((1 - f) * d) / math.sin(d) + B = np.sin(f * d) / math.sin(d) + x = A * math.cos(lat1) * math.cos(lon1) + B * math.cos(lat2) * math.cos(lon2) + y = A * math.cos(lat1) * math.sin(lon1) + B * math.cos(lat2) * math.sin(lon2) + z = A * math.sin(lat1) + B * math.sin(lat2) + lats = np.arctan2(z, np.sqrt(x ** 2 + y ** 2)) + lons = np.arctan2(y, x) + lons = map(math.degrees, lons.tolist()) + lats = map(math.degrees, lats.tolist()) # use ellipsoid formulas else: latpt = self.lat1 @@ -2772,15 +2805,18 @@ def points(self,npoints): azimuth = self.azimuth12 lons = [math.degrees(lonpt)] lats = [math.degrees(latpt)] - for n in range(npoints-2): - latptnew,lonptnew,alpha21=vinc_pt(self.f,self.a,latpt,lonpt,azimuth,incdist) - d,azimuth,a21=vinc_dist(self.f,self.a,latptnew,lonptnew,lat2,lon2) + for n in range(npoints - 2): + latptnew, lonptnew, alpha21 = vinc_pt(self.f, self.a, latpt, lonpt, azimuth, incdist) + d, azimuth, a21 = vinc_dist(self.f, self.a, latptnew, lonptnew, lat2, lon2) lats.append(math.degrees(latptnew)) lons.append(math.degrees(lonptnew)) - latpt = latptnew; lonpt = lonptnew + latpt = latptnew + lonpt = lonptnew lons.append(math.degrees(self.lon2)) lats.append(math.degrees(self.lat2)) - return lons,lats + return lons, lats + + # # --------------------------------------------------------------------- # | | @@ -2824,7 +2860,7 @@ def points(self,npoints): # | | # ---------------------------------------------------------------------- -def vinc_dist( f, a, phi1, lembda1, phi2, lembda2 ) : +def vinc_dist(f, a, phi1, lembda1, phi2, lembda2): """ Returns the distance between two geographic points on the ellipsoid @@ -2835,87 +2871,86 @@ def vinc_dist( f, a, phi1, lembda1, phi2, lembda2 ) : """ - if (abs( phi2 - phi1 ) < 1e-8) and ( abs( lembda2 - lembda1) < 1e-8 ) : + if (abs(phi2 - phi1) < 1e-8) and (abs(lembda2 - lembda1) < 1e-8): return 0.0, 0.0, 0.0 - two_pi = 2.0*math.pi + two_pi = 2.0 * math.pi b = a * (1.0 - f) - TanU1 = (1-f) * math.tan( phi1 ) - TanU2 = (1-f) * math.tan( phi2 ) + TanU1 = (1 - f) * math.tan(phi1) + TanU2 = (1 - f) * math.tan(phi2) U1 = math.atan(TanU1) U2 = math.atan(TanU2) lembda = lembda2 - lembda1 - last_lembda = -4000000.0 # an impossibe value + last_lembda = -4000000.0 # an impossibe value omega = lembda # Iterate the following equations, # until there is no significant change in lembda - while ( last_lembda < -3000000.0 or lembda != 0 and abs( (last_lembda - lembda)/lembda) > 1.0e-9 ) : + while (last_lembda < -3000000.0 or lembda != 0 and abs((last_lembda - lembda) / lembda) > 1.0e-9): + sqr_sin_sigma = pow(math.cos(U2) * math.sin(lembda), 2) + \ + pow((math.cos(U1) * math.sin(U2) - \ + math.sin(U1) * math.cos(U2) * math.cos(lembda)), 2) - sqr_sin_sigma = pow( math.cos(U2) * math.sin(lembda), 2) + \ - pow( (math.cos(U1) * math.sin(U2) - \ - math.sin(U1) * math.cos(U2) * math.cos(lembda) ), 2 ) - - Sin_sigma = math.sqrt( sqr_sin_sigma ) + Sin_sigma = math.sqrt(sqr_sin_sigma) Cos_sigma = math.sin(U1) * math.sin(U2) + math.cos(U1) * math.cos(U2) * math.cos(lembda) - sigma = math.atan2( Sin_sigma, Cos_sigma ) + sigma = math.atan2(Sin_sigma, Cos_sigma) Sin_alpha = math.cos(U1) * math.cos(U2) * math.sin(lembda) / math.sin(sigma) - alpha = math.asin( Sin_alpha ) + alpha = math.asin(Sin_alpha) - Cos2sigma_m = math.cos(sigma) - (2 * math.sin(U1) * math.sin(U2) / pow(math.cos(alpha), 2) ) + Cos2sigma_m = math.cos(sigma) - (2 * math.sin(U1) * math.sin(U2) / pow(math.cos(alpha), 2)) - C = (f/16) * pow(math.cos(alpha), 2) * (4 + f * (4 - 3 * pow(math.cos(alpha), 2))) + C = (f / 16) * pow(math.cos(alpha), 2) * (4 + f * (4 - 3 * pow(math.cos(alpha), 2))) last_lembda = lembda - lembda = omega + (1-C) * f * math.sin(alpha) * (sigma + C * math.sin(sigma) * \ - (Cos2sigma_m + C * math.cos(sigma) * (-1 + 2 * pow(Cos2sigma_m, 2) ))) - + lembda = omega + (1 - C) * f * math.sin(alpha) * (sigma + C * math.sin(sigma) * \ + (Cos2sigma_m + C * math.cos(sigma) * ( + -1 + 2 * pow(Cos2sigma_m, 2)))) - u2 = pow(math.cos(alpha),2) * (a*a-b*b) / (b*b) + u2 = pow(math.cos(alpha), 2) * (a * a - b * b) / (b * b) - A = 1 + (u2/16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))) + A = 1 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))) - B = (u2/1024) * (256 + u2 * (-128+ u2 * (74 - 47 * u2))) + B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) - delta_sigma = B * Sin_sigma * (Cos2sigma_m + (B/4) * \ - (Cos_sigma * (-1 + 2 * pow(Cos2sigma_m, 2) ) - \ - (B/6) * Cos2sigma_m * (-3 + 4 * sqr_sin_sigma) * \ - (-3 + 4 * pow(Cos2sigma_m,2 ) ))) + delta_sigma = B * Sin_sigma * (Cos2sigma_m + (B / 4) * \ + (Cos_sigma * (-1 + 2 * pow(Cos2sigma_m, 2)) - \ + (B / 6) * Cos2sigma_m * (-3 + 4 * sqr_sin_sigma) * \ + (-3 + 4 * pow(Cos2sigma_m, 2)))) s = b * A * (sigma - delta_sigma) - alpha12 = math.atan2( (math.cos(U2) * math.sin(lembda)), \ - (math.cos(U1) * math.sin(U2) - math.sin(U1) * math.cos(U2) * math.cos(lembda))) + alpha12 = math.atan2((math.cos(U2) * math.sin(lembda)), \ + (math.cos(U1) * math.sin(U2) - math.sin(U1) * math.cos(U2) * math.cos(lembda))) - alpha21 = math.atan2( (math.cos(U1) * math.sin(lembda)), \ - (-math.sin(U1) * math.cos(U2) + math.cos(U1) * math.sin(U2) * math.cos(lembda))) + alpha21 = math.atan2((math.cos(U1) * math.sin(lembda)), \ + (-math.sin(U1) * math.cos(U2) + math.cos(U1) * math.sin(U2) * math.cos(lembda))) - if ( alpha12 < 0.0 ) : - alpha12 = alpha12 + two_pi - if ( alpha12 > two_pi ) : + if (alpha12 < 0.0): + alpha12 = alpha12 + two_pi + if (alpha12 > two_pi): alpha12 = alpha12 - two_pi alpha21 = alpha21 + two_pi / 2.0 - if ( alpha21 < 0.0 ) : + if (alpha21 < 0.0): alpha21 = alpha21 + two_pi - if ( alpha21 > two_pi ) : + if (alpha21 > two_pi): alpha21 = alpha21 - two_pi - return s, alpha12, alpha21 + return s, alpha12, alpha21 # END of Vincenty's Inverse formulae -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- # Vincenty's Direct formulae | # Given: latitude and longitude of a point (phi1, lembda1) and | # the geodetic azimuth (alpha12) | @@ -2924,9 +2959,9 @@ def vinc_dist( f, a, phi1, lembda1, phi2, lembda2 ) : # Calculate: the latitude and longitude of the second point (phi2, lembda2) | # and the reverse azimuth (alpha21). | # | -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- -def vinc_pt( f, a, phi1, lembda1, alpha12, s ) : +def vinc_pt(f, a, phi1, lembda1, alpha12, s): """ Returns the lat and long of projected point and reverse azimuth @@ -2937,27 +2972,25 @@ def vinc_pt( f, a, phi1, lembda1, alpha12, s ) : """ + two_pi = 2.0 * math.pi - two_pi = 2.0*math.pi - - if ( alpha12 < 0.0 ) : + if (alpha12 < 0.0): alpha12 = alpha12 + two_pi - if ( alpha12 > two_pi ) : + if (alpha12 > two_pi): alpha12 = alpha12 - two_pi - b = a * (1.0 - f) - TanU1 = (1-f) * math.tan(phi1) - U1 = math.atan( TanU1 ) - sigma1 = math.atan2( TanU1, math.cos(alpha12) ) + TanU1 = (1 - f) * math.tan(phi1) + U1 = math.atan(TanU1) + sigma1 = math.atan2(TanU1, math.cos(alpha12)) Sinalpha = math.cos(U1) * math.sin(alpha12) cosalpha_sq = 1.0 - Sinalpha * Sinalpha - u2 = cosalpha_sq * (a * a - b * b ) / (b * b) + u2 = cosalpha_sq * (a * a - b * b) / (b * b) A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * \ - (320 - 175 * u2) ) ) - B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2) ) ) + (320 - 175 * u2))) + B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) # Starting with the approximation sigma = (s / (b * A)) @@ -2969,200 +3002,198 @@ def vinc_pt( f, a, phi1, lembda1, alpha12, s ) : # two_sigma_m , delta_sigma - while ( abs( (last_sigma - sigma) / sigma) > 1.0e-9 ) : - + while (abs((last_sigma - sigma) / sigma) > 1.0e-9): two_sigma_m = 2 * sigma1 + sigma - delta_sigma = B * math.sin(sigma) * ( math.cos(two_sigma_m) \ - + (B/4) * (math.cos(sigma) * \ - (-1 + 2 * math.pow( math.cos(two_sigma_m), 2 ) - \ - (B/6) * math.cos(two_sigma_m) * \ - (-3 + 4 * math.pow(math.sin(sigma), 2 )) * \ - (-3 + 4 * math.pow( math.cos (two_sigma_m), 2 ))))) + delta_sigma = B * math.sin(sigma) * (math.cos(two_sigma_m) \ + + (B / 4) * (math.cos(sigma) * \ + (-1 + 2 * math.pow(math.cos(two_sigma_m), 2) - \ + (B / 6) * math.cos(two_sigma_m) * \ + (-3 + 4 * math.pow(math.sin(sigma), 2)) * \ + (-3 + 4 * math.pow(math.cos(two_sigma_m), 2))))) last_sigma = sigma sigma = (s / (b * A)) + delta_sigma + phi2 = math.atan2((math.sin(U1) * math.cos(sigma) + math.cos(U1) * math.sin(sigma) * math.cos(alpha12)), \ + ((1 - f) * math.sqrt(math.pow(Sinalpha, 2) + \ + pow(math.sin(U1) * math.sin(sigma) - math.cos(U1) * math.cos( + sigma) * math.cos(alpha12), 2)))) - phi2 = math.atan2 ( (math.sin(U1) * math.cos(sigma) + math.cos(U1) * math.sin(sigma) * math.cos(alpha12) ), \ - ((1-f) * math.sqrt( math.pow(Sinalpha, 2) + \ - pow(math.sin(U1) * math.sin(sigma) - math.cos(U1) * math.cos(sigma) * math.cos(alpha12), 2)))) + lembda = math.atan2((math.sin(sigma) * math.sin(alpha12)), (math.cos(U1) * math.cos(sigma) - \ + math.sin(U1) * math.sin(sigma) * math.cos(alpha12))) + C = (f / 16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq)) - lembda = math.atan2( (math.sin(sigma) * math.sin(alpha12 )), (math.cos(U1) * math.cos(sigma) - \ - math.sin(U1) * math.sin(sigma) * math.cos(alpha12))) - - C = (f/16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq )) - - omega = lembda - (1-C) * f * Sinalpha * \ - (sigma + C * math.sin(sigma) * (math.cos(two_sigma_m) + \ - C * math.cos(sigma) * (-1 + 2 * math.pow(math.cos(two_sigma_m),2) ))) + omega = lembda - (1 - C) * f * Sinalpha * \ + (sigma + C * math.sin(sigma) * (math.cos(two_sigma_m) + \ + C * math.cos(sigma) * ( + -1 + 2 * math.pow(math.cos(two_sigma_m), 2)))) lembda2 = lembda1 + omega - alpha21 = math.atan2 ( Sinalpha, (-math.sin(U1) * math.sin(sigma) + \ - math.cos(U1) * math.cos(sigma) * math.cos(alpha12))) + alpha21 = math.atan2(Sinalpha, (-math.sin(U1) * math.sin(sigma) + \ + math.cos(U1) * math.cos(sigma) * math.cos(alpha12))) alpha21 = alpha21 + two_pi / 2.0 - if ( alpha21 < 0.0 ) : + if (alpha21 < 0.0): alpha21 = alpha21 + two_pi - if ( alpha21 > two_pi ) : + if (alpha21 > two_pi): alpha21 = alpha21 - two_pi - - return phi2, lembda2, alpha21 + return phi2, lembda2, alpha21 # END of Vincenty's Direct formulae -##--------------------------------------------------------------------------- -# Notes: -# -# * "The inverse formulae may give no solution over a line -# between two nearly antipodal points. This will occur when -# lembda ... is greater than pi in absolute value". (Vincenty, 1975) -# -# * In Vincenty (1975) L is used for the difference in longitude, -# however for consistency with other formulae in this Manual, -# omega is used here. -# -# * Variables specific to Vincenty's formulae are shown below, -# others common throughout the manual are shown in the Glossary. -# -# -# alpha = Azimuth of the geodesic at the equator -# U = Reduced latitude -# lembda = Difference in longitude on an auxiliary sphere (lembda1 & lembda2 -# are the geodetic longitudes of points 1 & 2) -# sigma = Angular distance on a sphere, from point 1 to point 2 -# sigma1 = Angular distance on a sphere, from the equator to point 1 -# sigma2 = Angular distance on a sphere, from the equator to point 2 -# sigma_m = Angular distance on a sphere, from the equator to the -# midpoint of the line from point 1 to point 2 -# u, A, B, C = Internal variables -# -# -# Sample Data -# -# Flinders Peak -# -37o57'03.72030" -# 144o25'29.52440" -# Buninyong -# -37o39'10.15610" -# 143o55'35.38390" -# Ellipsoidal Distance -# 54,972.271 m -# -# Forward Azimuth -# 306o52'05.37" -# -# Reverse Azimuth -# 127o10'25.07" -# -# -##******************************************************************* - -# Test driver - -#if __name__ == "__main__" : + ##--------------------------------------------------------------------------- + # Notes: + # + # * "The inverse formulae may give no solution over a line + # between two nearly antipodal points. This will occur when + # lembda ... is greater than pi in absolute value". (Vincenty, 1975) + # + # * In Vincenty (1975) L is used for the difference in longitude, + # however for consistency with other formulae in this Manual, + # omega is used here. + # + # * Variables specific to Vincenty's formulae are shown below, + # others common throughout the manual are shown in the Glossary. + # + # + # alpha = Azimuth of the geodesic at the equator + # U = Reduced latitude + # lembda = Difference in longitude on an auxiliary sphere (lembda1 & lembda2 + # are the geodetic longitudes of points 1 & 2) + # sigma = Angular distance on a sphere, from point 1 to point 2 + # sigma1 = Angular distance on a sphere, from the equator to point 1 + # sigma2 = Angular distance on a sphere, from the equator to point 2 + # sigma_m = Angular distance on a sphere, from the equator to the + # midpoint of the line from point 1 to point 2 + # u, A, B, C = Internal variables + # + # + # Sample Data + # + # Flinders Peak + # -37o57'03.72030" + # 144o25'29.52440" + # Buninyong + # -37o39'10.15610" + # 143o55'35.38390" + # Ellipsoidal Distance + # 54,972.271 m + # + # Forward Azimuth + # 306o52'05.37" + # + # Reverse Azimuth + # 127o10'25.07" + # + # + ##******************************************************************* + + # Test driver + + # if __name__ == "__main__" : ## WGS84 - #a = 6378137.0 - #b = 6356752.3142 - #f = (a-b)/a - - #print "\n Ellipsoidal major axis = %12.3f metres\n" % ( a ) - #print "\n Inverse flattening = %15.9f\n" % ( 1.0/f ) - - #print "\n Test Flinders Peak to Buninyon" - #print "\n ****************************** \n" - #phi1 = -(( 3.7203 / 60. + 57) / 60. + 37 ) - #lembda1 = ( 29.5244 / 60. + 25) / 60. + 144 - #print "\n Flinders Peak = %12.6f, %13.6f \n" % ( phi1, lembda1 ) - #deg = int(phi1) - #minn = int(abs( ( phi1 - deg) * 60.0 )) - #sec = abs(phi1 * 3600 - deg * 3600) - minn * 60 - #print " Flinders Peak = %3i\xF8%3i\' %6.3f\", " % ( deg, minn, sec ), - #deg = int(lembda1) - #minn = int(abs( ( lembda1 - deg) * 60.0 )) - #sec = abs(lembda1 * 3600 - deg * 3600) - minn * 60 - #print " %3i\xF8%3i\' %6.3f\" \n" % ( deg, minn, sec ) - - #phi2 = -(( 10.1561 / 60. + 39) / 60. + 37 ) - #lembda2 = ( 35.3839 / 60. + 55) / 60. + 143 - #print "\n Buninyon = %12.6f, %13.6f \n" % ( phi2, lembda2 ) - - #deg = int(phi2) - #minn = int(abs( ( phi2 - deg) * 60.0 )) - #sec = abs(phi2 * 3600 - deg * 3600) - minn * 60 - #print " Buninyon = %3i\xF8%3i\' %6.3f\", " % ( deg, minn, sec ), - #deg = int(lembda2) - #minn = int(abs( ( lembda2 - deg) * 60.0 )) - #sec = abs(lembda2 * 3600 - deg * 3600) - minn * 60 - #print " %3i\xF8%3i\' %6.3f\" \n" % ( deg, minn, sec ) - - #dist, alpha12, alpha21 = vinc_dist ( f, a, math.radians(phi1), math.radians(lembda1), math.radians(phi2), math.radians(lembda2) ) - - #alpha12 = math.degrees(alpha12) - #alpha21 = math.degrees(alpha21) - - #print "\n Ellipsoidal Distance = %15.3f metres\n should be 54972.271 m\n" % ( dist ) - #print "\n Forward and back azimuths = %15.6f, %15.6f \n" % ( alpha12, alpha21 ) - #deg = int(alpha12) - #minn =int( abs(( alpha12 - deg) * 60.0 ) ) - #sec = abs(alpha12 * 3600 - deg * 3600) - minn * 60 - #print " Forward azimuth = %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec ) - #deg = int(alpha21) - #minn =int(abs( ( alpha21 - deg) * 60.0 )) - #sec = abs(alpha21 * 3600 - deg * 3600) - minn * 60 - #print " Reverse azimuth = %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec ) + # a = 6378137.0 + # b = 6356752.3142 + # f = (a-b)/a + + # print("\n Ellipsoidal major axis = %12.3f metres\n" % ( a )) + # print("\n Inverse flattening = %15.9f\n" % ( 1.0/f )) + + # print("\n Test Flinders Peak to Buninyon") + # print("\n ****************************** \n") + # phi1 = -(( 3.7203 / 60. + 57) / 60. + 37 ) + # lembda1 = ( 29.5244 / 60. + 25) / 60. + 144 + # print("\n Flinders Peak = %12.6f, %13.6f \n" % ( phi1, lembda1 )) + # deg = int(phi1) + # minn = int(abs( ( phi1 - deg) * 60.0 )) + # sec = abs(phi1 * 3600 - deg * 3600) - minn * 60 + # print(" Flinders Peak = %3i\xF8%3i\' %6.3f\", % ( deg, minn, sec ),) + # deg = int(lembda1) + # minn = int(abs( ( lembda1 - deg) * 60.0 )) + # sec = abs(lembda1 * 3600 - deg * 3600) - minn * 60 + # print(" %3i\xF8%3i\' %6.3f\" \n" % ( deg, minn, sec )) + + # phi2 = -(( 10.1561 / 60. + 39) / 60. + 37 ) + # lembda2 = ( 35.3839 / 60. + 55) / 60. + 143 + # print("\n Buninyon = %12.6f, %13.6f \n" % ( phi2, lembda2 )) + + # deg = int(phi2) + # minn = int(abs( ( phi2 - deg) * 60.0 )) + # sec = abs(phi2 * 3600 - deg * 3600) - minn * 60 + # print(" Buninyon = %3i\xF8%3i\' %6.3f\", " % ( deg, minn, sec )) + # deg = int(lembda2) + # minn = int(abs( ( lembda2 - deg) * 60.0 )) + # sec = abs(lembda2 * 3600 - deg * 3600) - minn * 60 + # print(" %3i\xF8%3i\' %6.3f\" \n" % ( deg, minn, sec )) + + # dist, alpha12, alpha21 = vinc_dist ( f, a, math.radians(phi1), math.radians(lembda1), math.radians(phi2), math.radians(lembda2) ) + + # alpha12 = math.degrees(alpha12) + # alpha21 = math.degrees(alpha21) + + # print("\n Ellipsoidal Distance = %15.3f metres\n should be 54972.271 m\n" % ( dist )) + # print("\n Forward and back azimuths = %15.6f, %15.6f \n" % ( alpha12, alpha21 )) + # deg = int(alpha12) + # minn =int( abs(( alpha12 - deg) * 60.0 ) ) + # sec = abs(alpha12 * 3600 - deg * 3600) - minn * 60 + # print(" Forward azimuth = %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec )) + # deg = int(alpha21) + # minn =int(abs( ( alpha21 - deg) * 60.0 )) + # sec = abs(alpha21 * 3600 - deg * 3600) - minn * 60 + # print(" Reverse azimuth = %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec )) ## Test the direct function */ - #phi1 = -(( 3.7203 / 60. + 57) / 60. + 37 ) - #lembda1 = ( 29.5244 / 60. + 25) / 60. + 144 - #dist = 54972.271 - #alpha12 = ( 5.37 / 60. + 52) / 60. + 306 - #phi2 = lembda2 = 0.0 - #alpha21 = 0.0 - - ephi2, lembda2, alpha21 = vinc_pt ( f, a, math.radians(phi1), math.radians(lembda1), math.radians(alpha12), dist ) - - #phi2 = math.degrees(phi2) - #lembda2 = math.degrees(lembda2) - #alpha21 = math.degrees(alpha21) - - #print "\n Projected point =%11.6f, %13.6f \n" % ( phi2, lembda2 ) - #deg = int(phi2) - #minn =int(abs( ( phi2 - deg) * 60.0 )) - #sec = abs( phi2 * 3600 - deg * 3600) - minn * 60 - #print " Projected Point = %3i\xF8%3i\' %6.3f\", " % ( deg, minn, sec ), - #deg = int(lembda2) - #minn =int(abs( ( lembda2 - deg) * 60.0 )) - #sec = abs(lembda2 * 3600 - deg * 3600) - minn * 60 - #print " %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec ) - #print " Should be Buninyon \n" - #print "\n Reverse azimuth = %10.6f \n" % ( alpha21 ) - #deg = int(alpha21) - #minn =int(abs( ( alpha21 - deg) * 60.0 )) - #sec = abs(alpha21 * 3600 - deg * 3600) - minn * 60 - #print " Reverse azimuth = %3i\xF8%3i\' %6.3f\"\n\n" % ( deg, minn, sec ) + # phi1 = -(( 3.7203 / 60. + 57) / 60. + 37 ) + # lembda1 = ( 29.5244 / 60. + 25) / 60. + 144 + # dist = 54972.271 + # alpha12 = ( 5.37 / 60. + 52) / 60. + 306 + # phi2 = lembda2 = 0.0 + # alpha21 = 0.0 + + ephi2, lembda2, alpha21 = vinc_pt(f, a, math.radians(phi1), math.radians(lembda1), math.radians(alpha12), dist) + + # phi2 = math.degrees(phi2) + # lembda2 = math.degrees(lembda2) + # alpha21 = math.degrees(alpha21) + + # print("\n Projected point =%11.6f, %13.6f \n" % ( phi2, lembda2 )) + # deg = int(phi2) + # minn =int(abs( ( phi2 - deg) * 60.0 )) + # sec = abs( phi2 * 3600 - deg * 3600) - minn * 60 + # print(" Projected Point = %3i\xF8%3i\' %6.3f\", " % ( deg, minn, sec ),) + # deg = int(lembda2) + # minn =int(abs( ( lembda2 - deg) * 60.0 )) + # sec = abs(lembda2 * 3600 - deg * 3600) - minn * 60 + # print(" %3i\xF8%3i\' %6.3f\"\n" % ( deg, minn, sec )) + # print(" Should be Buninyon \n") + # print("\n Reverse azimuth = %10.6f \n" % ( alpha21 )) + # deg = int(alpha21) + # minn =int(abs( ( alpha21 - deg) * 60.0 )) + # sec = abs(alpha21 * 3600 - deg * 3600) - minn * 60 + # print(" Reverse azimuth = %3i\xF8%3i\' %6.3f\"\n\n" % ( deg, minn, sec )) ## lat/lon of New York - #lat1 = 40.78 - #lon1 = -73.98 + # lat1 = 40.78 + # lon1 = -73.98 ## lat/lon of London. - #lat2 = 51.53 - #lon2 = 0.08 - #print 'New York to London:' - #gc = GreatCircle((2*a+b)/3.,(2*a+b)/3.,lon1,lat1,lon2,lat2) - #print 'geodesic distance using a sphere with WGS84 mean radius = ',gc.distance - #print 'lon/lat for 10 equally spaced points along geodesic:' - #lons,lats = gc.points(10) - #for lon,lat in zip(lons,lats): - #print lon,lat - #gc = GreatCircle(a,b,lon1,lat1,lon2,lat2) - #print 'geodesic distance using WGS84 ellipsoid = ',gc.distance - #print 'lon/lat for 10 equally spaced points along geodesic:' - #lons,lats = gc.points(10) - #for lon,lat in zip(lons,lats): - #print lon,lat + # lat2 = 51.53 + # lon2 = 0.08 + # print('New York to London:') + # gc = GreatCircle((2*a+b)/3.,(2*a+b)/3.,lon1,lat1,lon2,lat2) + # print('geodesic distance using a sphere with WGS84 mean radius = ',gc.distance) + # print('lon/lat for 10 equally spaced points along geodesic:') + # lons,lats = gc.points(10) + # for lon,lat in zip(lons,lats): + # print(lon, lat) + # gc = GreatCircle(a,b,lon1,lat1,lon2,lat2) + # print('geodesic distance using WGS84 ellipsoid = ', gc.distance) + # print('lon/lat for 10 equally spaced points along geodesic:') + # lons,lats = gc.points(10) + # for lon,lat in zip(lons,lats): + # print(lon,lat) diff --git a/reflexible/plotting.py b/reflexible/plotting.py index 1cf3934..8f9473c 100644 --- a/reflexible/plotting.py +++ b/reflexible/plotting.py @@ -263,12 +263,12 @@ def plot_sensitivity(H, data, \ del FIGURE.ax.collections[FIGURE.indices.collections:] del FIGURE.ax.lines[FIGURE.indices.lines:] - if dropm != None: + if dropm is not None: try: del m plt.close('all') except: - print 'could not drop m' + print('could not drop m') # # make tick lables smaller mpl.rcParams['xtick.labelsize'] = 6 @@ -358,7 +358,7 @@ def plot_sensitivity(H, data, \ colmap = _gen_flexpart_colormap() colmap.set_over(color='k', alpha=0.8) # # Plotting METHODS (pcolormesh now default, imshow is smoother) - # print topodat.max(), topodat.min(), topodat.shape + # print(topodat.max(), topodat.min(), topodat.shape) if method == 'imshow': im = m.imshow(topodat, cmap=colmap, zorder=-1, norm=mpl.colors.LogNorm(vmin=clevs[0], @@ -690,7 +690,7 @@ def _gen_flexpart_colormap(ctbfile=None, colors=None): try: colors = np.loadtxt(ctbfile) except: - print "WARNING: cannot load ctbfile. using colors" + print("WARNING: cannot load ctbfile. using colors") if colors: name = 'user_colormap' if not colors: diff --git a/reflexible/tests/test_plotting.py b/reflexible/tests/test_plotting.py index 0b0c0f0..903658a 100644 --- a/reflexible/tests/test_plotting.py +++ b/reflexible/tests/test_plotting.py @@ -1,10 +1,8 @@ #!/usr/bin/env python from __future__ import absolute_import - -import os.path +from __future__ import print_function import pytest -import numpy as np import reflexible.plotting as pf import reflexible as rf @@ -12,9 +10,11 @@ # tuples locating test data, nested(True) and global(False) test_datasets = [('Fwd1_V10.0', False), ('Fwd1_V10.0', True)] + def monotonically_increasing(l): return all(x < y for x, y in zip(l, l[1:])) + class Dataset: def __init__(self, fp_dataset): self.fp_name = fp_dataset[0] @@ -35,12 +35,10 @@ class TestPlotting: @pytest.fixture(autouse=True, params=test_datasets) def setup(self, request): pass - #dataset = Dataset(request.param) - #self.H, self.fp_path, self.nc_path = dataset.setup() - #request.addfinalizer(dataset.cleanup) + # dataset = Dataset(request.param) + # self.H, self.fp_path, self.nc_path = dataset.setup() + # request.addfinalizer(dataset.cleanup) def test_flexpart_colormap(self): from matplotlib.colors import ListedColormap assert isinstance(pf._gen_flexpart_colormap(), ListedColormap) - - \ No newline at end of file From 6b4c8227cf8822e1938130437716936b86faab6f Mon Sep 17 00:00:00 2001 From: Francesc Alted Date: Tue, 13 Sep 2016 18:01:25 +0200 Subject: [PATCH 2/4] All tests passes now for Python 3.5 --- reflexible/base_read.py | 88 ++++++++++---------- reflexible/conv2netcdf4/flexpart_read.py | 21 ++--- reflexible/conv2netcdf4/tests/test_basics.py | 18 ++-- reflexible/data_structures.py | 4 +- 4 files changed, 66 insertions(+), 65 deletions(-) diff --git a/reflexible/base_read.py b/reflexible/base_read.py index 1edf10a..9c98b95 100755 --- a/reflexible/base_read.py +++ b/reflexible/base_read.py @@ -1,4 +1,3 @@ - """ functions for reading fp data in reflexible. @@ -36,20 +35,15 @@ import os import datetime as dt - # Dependencies: # Numpy import numpy as np - from .data_structures import Trajectory - - -def read_trajectories(H, trajfile='trajectories.txt', \ - ncluster=5, \ - ageclasses=20): +def read_trajectories(H, trajfile='trajectories.txt', ncluster=5, + ageclasses=20): """ Reads the trajectories.txt file in a FLEXPART run output directory. @@ -126,13 +120,12 @@ def read_trajectories(H, trajfile='trajectories.txt', \ if isinstance(H, str): try: - alltraj = file(H, 'r').readlines() + alltraj = open(H, 'r').readlines() except: raise IOError('Could not open file: %s' % H) else: path = H.fp_path - alltraj = file(os.path.join(path, trajfile), 'r').readlines() - + alltraj = open(os.path.join(path, trajfile), 'r').readlines() try: ibdate, ibtime, model, version = alltraj[0].strip().split()[:4] @@ -142,33 +135,37 @@ def read_trajectories(H, trajfile='trajectories.txt', \ version = 'V.x' tdelta = dt.datetime.strptime(ibdate + ibtime.zfill(6), '%Y%m%d%H%M%S') numpoint = int(alltraj[2].strip()) - + RelTraj = Trajectory() Trajectories = [] # format change?? Try block below because some trajectories.txt # files have linebreaks, but seems more recent FP v10. does not. try: for i in range(3, 3 + (numpoint * 2), 3): - i1, i2, xp1, yp1, xp2 , = tuple([float(j) for j in alltraj[i].strip().split()]) - yp2, zp1, zp2, k, npart , = tuple([float(j) for j in alltraj[i+1].strip().split()]) + i1, i2, xp1, yp1, xp2, = tuple((float(j) for j in + alltraj[i].strip().split())) + yp2, zp1, zp2, k, npart, = tuple((float(j) for j in + alltraj[i + 1].strip().split())) except: for i in range(3, 3 + (numpoint * 2), 2): - i1, i2, xp1, yp1, xp2, yp2, zp1, zp2, k, npart , = \ - tuple([float(j) for j in alltraj[i].strip().split()]) + i1, i2, xp1, yp1, xp2, yp2, zp1, zp2, k, npart, = \ + tuple([float(j) for j in alltraj[i].strip().split()]) itimerel1 = tdelta + dt.timedelta(seconds=i1) itimerel2 = tdelta + dt.timedelta(seconds=i2) Xp = (xp1 + xp2) / 2 Yp = (yp1 + yp2) / 2 Zp = (zp1 + zp2) / 2 - RelTraj[alltraj[i + 1].strip()] = np.array((itimerel1, itimerel2, Xp, Yp, Zp, k, npart)) + RelTraj[alltraj[i + 1].strip()] = np.array( + (itimerel1, itimerel2, Xp, Yp, Zp, k, npart)) for i in range(3 + (numpoint * 3), len(alltraj)): raw = alltraj[i] - FMT = [0, 5, 8, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 6, 6] + ncluster * [8, 8, 7, 6, 8] + FMT = [0, 5, 8, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 6, 6] + \ + ncluster * [8, 8, 7, 6, 8] - data = [raw[sum(FMT[:ii]):sum(FMT[:ii + 1])] for ii in range(1, len(FMT) - 1)] + \ - [raw[sum(FMT[:-1]):]] + data = [raw[sum(FMT[:ii]):sum(FMT[:ii + 1])] + for ii in range(1, len(FMT) - 1)] + [raw[sum(FMT[:-1]):]] ### FIX ### # ## To get rid of '******' that is now in trajectories.txt data = [float(r.replace('********', 'NaN')) for r in data] @@ -180,28 +177,31 @@ def read_trajectories(H, trajfile='trajectories.txt', \ RelTraj['date'] = tdelta RelTraj['Trajectories'] = data RelTraj['labels'] = \ - ['release number', 'seconds prior to release', 'lon', 'lat', 'height', 'mean topography', \ - 'mean mixing height', 'mean tropopause height', 'mean PV index', \ - 'rms distance', 'rms', 'zrms distance', 'zrms', \ - 'fraction mixing layer', 'fraction PV<2pvu', 'fraction in troposphere'] + \ - ncluster * ['xcluster', 'ycluster', 'zcluster', 'fcluster', 'rmscluster'] + ['release number', 'seconds prior to release', 'lon', 'lat', + 'height', 'mean topography', + 'mean mixing height', 'mean tropopause height', 'mean PV index', + 'rms distance', 'rms', 'zrms distance', 'zrms', + 'fraction mixing layer', 'fraction PV<2pvu', + 'fraction in troposphere'] + \ + ncluster * ['xcluster', 'ycluster', 'zcluster', 'fcluster', + 'rmscluster'] RelTraj['info'] = \ - """ - Returns a dictionary: - R['Trajectories'] = array_of_floats( - releasenum,it1,xi,yi,zi,topoi,hmixi,tropoi,pvi, - rmsdisti,rmsi,zrmsdisti,zrmsi,hfri,pvfri,trfri, - (xclusti(k),yclusti(k),zclusti(k),fclusti(k),rmsclusti(k),k=1,5)) - - R['RELEASE_ID'] = (dt_i1,dt_i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart) - R['info'] = this message - - To plot a trajectory timeseries: - RT = read_trajectories(H) - T = RT['Trajectories'] - rel = 1 - t = T[np.where(T[:,0]==rel),:][0] - plt.plot(t[:,1],t[:,14]) - plt.savefig('trajectories.png') - """ - return RelTraj \ No newline at end of file + """ + Returns a dictionary: + R['Trajectories'] = array_of_floats( + releasenum,it1,xi,yi,zi,topoi,hmixi,tropoi,pvi, + rmsdisti,rmsi,zrmsdisti,zrmsi,hfri,pvfri,trfri, + (xclusti(k),yclusti(k),zclusti(k),fclusti(k),rmsclusti(k),k=1,5)) + + R['RELEASE_ID'] = (dt_i1,dt_i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart) + R['info'] = this message + + To plot a trajectory timeseries: + RT = read_trajectories(H) + T = RT['Trajectories'] + rel = 1 + t = T[np.where(T[:,0]==rel),:][0] + plt.plot(t[:,1],t[:,14]) + plt.savefig('trajectories.png') + """ + return RelTraj diff --git a/reflexible/conv2netcdf4/flexpart_read.py b/reflexible/conv2netcdf4/flexpart_read.py index 0a27dc2..8750605 100644 --- a/reflexible/conv2netcdf4/flexpart_read.py +++ b/reflexible/conv2netcdf4/flexpart_read.py @@ -390,7 +390,8 @@ def read_trajectories(H, trajfile='trajectories.txt', ibdate, ibtime = alltraj[0].strip().split()[:2] model = 'Flexpart' version = 'V.x' - dt = datetime.datetime.strptime(ibdate + ibtime.zfill(6), '%Y%m%d%H%M%S') + dt = datetime.datetime.strptime(ibdate + ibtime.zfill(6), + '%Y%m%d%H%M%S') numpoint = int(alltraj[2].strip()) # Fill a dictionary with the Release points and information keyed by name # RelTraj['RelTraj_ID'] = (i1,i2,xp1,yp1,xp2,yp2,zp1,zp2,k,npart) @@ -702,12 +703,12 @@ def _get_header_version(bf): ret = bf.tell() bf.seek(12) # start of version string - version = bf.read('13S') + version = bf.read('13S').decode() # Somewhere in version 9.2 beta, the version length changed to 29 # However, one *must* check which is the final size for this if '9.2 b' in version: bf.seek(12) - version = bf.read('29S') + version = bf.read('29S').decode() bf.seek(ret) return version @@ -874,7 +875,7 @@ def read_header(pathname, **kwargs): h['jjjjmmdd'] = bf.read('i') h['hhmmss'] = bf.read('i') junk.append(bf.read('2i')) - h['nspec'] = bf.read('i') / 3 # why!? + h['nspec'] = bf.read('i') // 3 # why!? h['numpointspec'] = bf.read('i') junk.append(bf.read('2i')) @@ -896,15 +897,15 @@ def read_header(pathname, **kwargs): else: junk.append(bf.read('i')) h['wetdep'].append( - ''.join([bf.read('c') for i in range(10)]).strip()) + ''.join([bf.read('c').decode() for i in range(10)]).strip()) junk.append(bf.read('2i')) junk.append(bf.read('i')) h['drydep'].append( - ''.join([bf.read('c') for i in range(10)]).strip()) + ''.join([bf.read('c').decode() for i in range(10)]).strip()) junk.append(bf.read('2i')) h['nz_list'].append(bf.read('i')) h['species'].append( - ''.join([bf.read('c') for i in range(10)]).strip()) + ''.join([bf.read('c').decode() for i in range(10)]).strip()) junk.append(bf.read('2i')) if 'V6' in version: @@ -922,7 +923,7 @@ def read_header(pathname, **kwargs): I = {2: 'kindz', 3: 'xp1', 4: 'yp1', 5: 'xp2', 6: 'yp2', 7: 'zpoint1', 8: 'zpoint2', 9: 'npart', 10: 'mpart'} - for k, v in I.iteritems(): + for k, v in I.items(): # create zero-filled lists in H dict h[v] = np.zeros(h['numpoint']) @@ -960,9 +961,9 @@ def read_header(pathname, **kwargs): gt = bf.tell() + l # create 'goto' point sp = '' # collect the characters for the compoint - while re.search("\w", bf.read('c')): + while re.search(b"\w", bf.read('c')): bf.seek(-1, 1) - sp = sp + bf.read('c') + sp = sp + bf.read('c').decode() bf.seek(gt) # skip ahead to gt point diff --git a/reflexible/conv2netcdf4/tests/test_basics.py b/reflexible/conv2netcdf4/tests/test_basics.py index 2faea8a..b737364 100644 --- a/reflexible/conv2netcdf4/tests/test_basics.py +++ b/reflexible/conv2netcdf4/tests/test_basics.py @@ -16,7 +16,7 @@ def __init__(self, fp_name): self.fp_path = rf.datasets[fp_name] def setup(self, tmpdir): - self.tmpdir = tmpdir # bring the fixture to the Dataset instance + self.tmpdir = tmpdir # bring the fixture to the Dataset instance self.H = conv.Header(self.fp_path) self.nc_path = tmpdir.join("%s.nc" % self.fp_name).strpath return self.H, self.fp_path, self.nc_path @@ -38,15 +38,15 @@ def test_nc_create(self): def test_read_grid(self): FD = conv.read_grid(self.H, time_ret=0, nspec_ret=0) - fdkeys = sorted(FD.keys()) - assert fdkeys == ['grid_dates', 'options', (0, '20070121100000')] + fdkeys = sorted([k for k in FD.keys() if type(k) is str]) + assert fdkeys == ['grid_dates', 'options'] fdkeys_ = sorted(FD[(0, '20070121100000')].keys()) assert fdkeys_ == fd_keys def test_H_read_grid(self): self.H.read_grid(time_ret=0, nspec_ret=0) - fdkeys = sorted(self.H.FD.keys()) - assert fdkeys == ['grid_dates', 'options', (0, '20070121100000')] + fdkeys = sorted([k for k in self.H.FD.keys() if type(k) is str]) + assert fdkeys == ['grid_dates', 'options'] fdkeys_ = sorted(self.H.FD[(0, '20070121100000')].keys()) assert fdkeys_ == fd_keys @@ -74,14 +74,14 @@ def test_fill_backward(self): def test_read_grid(self): FD = conv.read_grid(self.H, time_ret=0, nspec_ret=0) - fdkeys = sorted(FD.keys()) - assert fdkeys == ['grid_dates', 'options', (0, self.H.available_dates[0])] + fdkeys = sorted([k for k in FD.keys() if type(k) is str]) + assert fdkeys == ['grid_dates', 'options'] fdkeys_ = sorted(FD[(0, self.H.available_dates[0])].keys()) assert fdkeys_ == fd_keys def test_H_read_grid(self): self.H.read_grid(time_ret=0, nspec_ret=0) - fdkeys = sorted(self.H.FD.keys()) - assert fdkeys == ['grid_dates', 'options', (0, self.H.available_dates[0])] + fdkeys = sorted([k for k in self.H.FD.keys() if type(k) is str]) + assert fdkeys == ['grid_dates', 'options'] fdkeys_ = sorted(self.H.FD[(0, self.H.available_dates[0])].keys()) assert fdkeys_ == fd_keys diff --git a/reflexible/data_structures.py b/reflexible/data_structures.py index b28b9de..37a47c4 100755 --- a/reflexible/data_structures.py +++ b/reflexible/data_structures.py @@ -735,11 +735,11 @@ def __init__(self, **options): self._overrides = options # set the default options as attributes - for key, value in self._OPTIONS.iteritems(): + for key, value in self._OPTIONS.items(): setattr(self, key.lower(), value[0]) # override the attributes with options - for key, value in options.iteritems(): + for key, value in options.items(): setattr(self, key.lower(), value) if self.ibdate is None: From 48e69c0697dccdaf93de19d37d83306751ffaf7e Mon Sep 17 00:00:00 2001 From: Francesc Alted Date: Tue, 13 Sep 2016 18:09:04 +0200 Subject: [PATCH 3/4] Fixed some other remaining prints --- reflexible/__init__.py | 16 +++++++--------- reflexible/conv2netcdf4/helpers.py | 4 +++- reflexible/scripts/readpartpositions.py | 11 ++++++----- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/reflexible/__init__.py b/reflexible/__init__.py index 642f66b..e84f695 100755 --- a/reflexible/__init__.py +++ b/reflexible/__init__.py @@ -38,14 +38,10 @@ This package follows creative commons usage. """ -#print WELCOME import os -import datetime as dt - from .version import __version__ - this_dir = __path__[0] @@ -70,6 +66,7 @@ def _get_hg_description(path_): except subprocess.CalledProcessError: # not in hg repo pass + _hg_description = _get_hg_description(this_dir) @@ -94,18 +91,19 @@ def print_versions(): datasets = { 'Fwd1_V10.0': os.path.join(this_dir, "uio_examples/Fwd1_V10.0"), 'Bwd1_V9.02': os.path.join(this_dir, "uio_examples/Bwd1_V9.02/outputs"), - 'Bwd2_V9.2beta': os.path.join(this_dir, "uio_examples/Bwd2_V9.2beta/outputs"), + 'Bwd2_V9.2beta': os.path.join(this_dir, + "uio_examples/Bwd2_V9.2beta/outputs"), 'Fwd1_V9.02': os.path.join(this_dir, "uio_examples/Fwd1_V9.02/outputs"), 'Fwd2_V9.02': os.path.join(this_dir, "uio_examples/Fwd2_V9.02/outputs"), - 'HelloWorld_V9.02': os.path.join(this_dir, "uio_examples/HelloWorld_V9.02/outputs"), - } - + 'HelloWorld_V9.02': os.path.join(this_dir, + "uio_examples/HelloWorld_V9.02/outputs"), +} # Import the public functions here from .tests.all import test from . import conv2netcdf4 from .scripts.create_ncfile import create_ncfile from .data_structures import (Header, Command, Release, - ReleasePoint, Ageclass) + ReleasePoint, Ageclass) from .base_read import read_trajectories diff --git a/reflexible/conv2netcdf4/helpers.py b/reflexible/conv2netcdf4/helpers.py index 9062880..ff58b20 100644 --- a/reflexible/conv2netcdf4/helpers.py +++ b/reflexible/conv2netcdf4/helpers.py @@ -1,5 +1,7 @@ ########### HELPER FUNCTIONS ########## +from __future__ import print_function + import datetime import numpy as np @@ -60,7 +62,7 @@ def _datarange(H, G, index=None): zpmax = np.max(G[:, :, 0, i] / Heightnn[:, :, 0]) if zpmax > fpmax: fpmax = zpmax - # print fpmax + # print(fpmax) tcmax = np.max(G) return ((0, fpmax), (0, tcmax)) diff --git a/reflexible/scripts/readpartpositions.py b/reflexible/scripts/readpartpositions.py index 23eca13..d080232 100644 --- a/reflexible/scripts/readpartpositions.py +++ b/reflexible/scripts/readpartpositions.py @@ -15,7 +15,8 @@ Date: 2015-01-30 """ -import sys +from __future__ import print_function + from time import time import numpy as np @@ -124,10 +125,10 @@ def main(): t0 = time() readout = readpartpositions(args.datafile, H.nspec, H.outlon0, H.outlat0, H.dxout, H.dyout, args.dataframe) if args.timming: - print "Time for reading the data: %.3fs" % (time()-t0,) - print "records read:", len(readout) - print "3 first records:\n", repr(readout[:3]) - print "3 last records:\n", repr(readout[-3:]) + print("Time for reading the data: %.3fs" % (time()-t0,)) + print("records read:", len(readout)) + print("3 first records:\n", repr(readout[:3])) + print("3 last records:\n", repr(readout[-3:])) if __name__ == "__main__": From b4134781c9dbe975276aeb5063edd3f1bbfddb14 Mon Sep 17 00:00:00 2001 From: Francesc Alted Date: Tue, 13 Sep 2016 18:18:39 +0200 Subject: [PATCH 4/4] Added python 3.5 to travis --- .travis.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.travis.yml b/.travis.yml index 2121326..c9bf976 100644 --- a/.travis.yml +++ b/.travis.yml @@ -43,6 +43,7 @@ os: python: - 2.7 + - 3.5 script: - py.test reflexible -v --cov reflexible --cov-report term-missing