Changeset 693 in lmdz_wrf


Ignore:
Timestamp:
Feb 11, 2016, 3:29:13 PM (9 years ago)
Author:
lfita
Message:

Adding `Partialmap_Entiremap': Function to transform from a partial global map (e.g.: only land points) to an entire one
Adding `lonlatProj': Function to compute longitudes and latitudes for a given projection

Location:
trunk/tools
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var.py

    r674 r693  
    3434  'maskvar',                                                                         \
    3535  'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
    36   'remapnn',                                                                         \
     36  'Partialmap_Entiremap', 'remapnn',                                                 \
    3737  'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'selvar', 'setvar_asciivalues',     \
    3838  'sorttimesmat', 'spacemean', 'SpatialWeightedMean', 'statcompare_files', 'submns', \
    … …  
    207207elif oper == 'remapnn':
    208208    ncvar.remapnn(opts.values, opts.ncfile, opts.varname)
     209elif oper == 'Partialmap_Entiremap':
     210    ncvar.Partialmap_Entiremap(opts.values, opts.ncfile, opts.varname)
    209211elif oper == 'seasmean':
    210212    ncvar.seasmean(timename, opts.ncfile, opts.varname)

  • trunk/tools/nc_var_tools.py

    r691 r693  
    1774717747        Nlat = int(np.pi * EarthR / dyres) + 1
    1774817748
     17749        if np.mod(Nlon,2) == 0: Nlon = Nlon + 1
     17750        if np.mod(Nlat,2) == 0: Nlat = Nlat + 1
     17751
    1774917752        lonProj = np.zeros((Nlon), dtype=np.float)
    1775017753        latProj = np.zeros((Nlat), dtype=np.float)
    … …  
    1782117824        newdim = onewnc.createDimension('y', dimy)
    1782217825        newdim = onewnc.createDimension('Npts', Npts)
    17823         newdim = onewnc.createDimension('Ncoord', 2)
    17824         newdim = onewnc.createDimension('xcirc', Nlonmaxcirc)
    17825         newdim = onewnc.createDimension('ycirc', Nlat)
     17826        if project == 'lonlat_dxdyFix':
     17827            newdim = onewnc.createDimension('xcirc', Nlonmaxcirc)
     17828            newdim = onewnc.createDimension('ycirc', Nlat)
    1782617829
    1782717830# Variables
    … …  
    1783217835        newvar.setncattr('_CoordinateAxisType', 'Lon')   
    1783317836
    17834 # Longitudes sorted by latitudes
    17835         if len(Npoints.shape) == 1:
    17836             newvar = onewnc.createVariable('Nloncirc', 'i', ('ycirc'))
    17837             newvar[:] = Nloncirc
    17838             basicvardef(newvar, 'Nloncirc', 'Number of longitudes by latitude circs', '-')
    17839 
    17840             newvar = onewnc.createVariable('loncirc', 'f8', ('ycirc','xcirc'))
    17841             basicvardef(newvar, 'longitude', 'Longitude by latitude circs', 'degrees_East')
    17842             for iL in range(Nlat):
    17843                 newvar[iL,0:Nloncirc[iL]] = np.asarray(lonProjcirc[iL])
    17844 
    1784517837        newvar = onewnc.createVariable('lat', 'f8', ('y'))
    1784617838        basicvardef(newvar, 'latitude', 'Latitude', 'degrees_North')
    … …  
    1784917841        newvar.setncattr('_CoordinateAxisType', 'Lat')   
    1785017842
    17851         newvar = onewnc.createVariable('latcirc', 'f8', ('ycirc'))
    17852         basicvardef(newvar, 'latitude', 'Latitude', 'degrees_North')
    17853         newvar[:] = latProjcirc
     17843# Specific variables
     17844        if project == 'lonlat_dxdyFix':
     17845            newvar = onewnc.createVariable('Nloncirc', 'i', ('ycirc'))
     17846            newvar[:] = Nloncirc
     17847            basicvardef(newvar, 'Nloncirc', 'Number of longitudes by latitude circs', '-')
     17848
     17849            newvar = onewnc.createVariable('loncirc', 'f8', ('ycirc','xcirc'))
     17850            basicvardef(newvar, 'longitude', 'Longitude by latitude circs', 'degrees_East')
     17851            for iL in range(Nlat):
     17852                newvar[iL,0:Nloncirc[iL]] = np.asarray(lonProjcirc[iL])
     17853
     17854            newvar = onewnc.createVariable('latcirc', 'f8', ('ycirc'))
     17855            basicvardef(newvar, 'latitude', 'Latitude', 'degrees_North')
     17856            newvar[:] = latProjcirc
    1785417857
    1785517858        if len(Npoints.shape) == 2:
    … …  
    1787617879        print fname + "Successfull written of '" + ofile + "' !!'"
    1787717880
    17878     if len(Npoints.shape) == 2:
     17881    if project == 'lonlat':
    1787917882        return lonProj, latProj
    17880     else:
     17883    elif project == 'lonlat_dxdyFix':
    1788117884        return lonProjcirc, latProj
    1788217885
    17883 #lonlatProj(100000.,100000.,'lonlat_dxdyFix',True)
     17886#lonlatProj(100000.,100000.,'lonlat',True)
    1788417887#quit()
    1788517888
    1788617889def Partialmap_Entiremap(values, filen, varn):
    1788717890    """ Function to transform from a partial global map (e.g.: only land points) to an entire one
    17888       values= [lonmame],[latname],[fillVal],[resolution]
     17891      values= [lonmame],[latname],[fillVal],[resolution],[kind]
    1788917892        [lonname]: name of the longitude variable
    1789017893        [latname]: name of the latitude variable
    … …  
    1789617899           Integer32 = -99999
    1789717900        [resolution]: resolution of the map
     17901        [kind]: kind of target projection
     17902         'lonlat': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator
     17903         'lonlat_dxdyFix': projection with a fixed grid distance
    1789817904      filen= name of the netCDF file
    1789917905      varn= name of the variable
    … …  
    1790417910    fname = 'Partialmap_Entiremap'
    1790517911
     17912    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind]'
     17913    check_arguments(fname, values, arguments, ',')
     17914
    1790617915    if values == 'h':
    1790717916        print fname + '_____________________________________________________________'
    … …  
    1791017919
    1791117920    ofile = 'EntireGlobalMap.nc'
    17912 
    17913     lonlatProjfile = filen.split('.')[0] + '_lonlatProj.nc'
    1791417921
    1791517922    onc = NetCDFFile(filen, 'r')
    … …  
    1792417931    fval0 = values.split(',')[2]
    1792517932    resolution = np.float(values.split(',')[3])
     17933    kind = values.split(',')[4]
     17934
     17935# Input reference lonlat map
     17936    lonlatProjfile = filen.split('.')[0] + '_' + kind + '.nc'
    1792617937
    1792717938    if not onc.variables.has_key(lonname):
    … …  
    1799518006
    1799618007    if not os.path.isfile(lonlatProjfile):
    17997         print lonlatProjfile
    17998         quit(-1)
    17999         lonmap, latmap = lonlatProj(resolution, resolution, 'lonlat_dxdyFix', True)
    18000         sub.call(['mv','lonlatProj.nc',lonlatProjfile], shell=True)
    18001         Ntotvals = len(lonmap)
    18002     else:
    18003         oproj = NetCDFFile(lonlatProjfile, 'r')
     18008        print warnmsg
     18009        print '  ' + fname + ": no reference map '" + lonlatProjfile + "' found !!"
     18010        print '  creating it via: lonlatProj(', resolution, ',', resolution,         \
     18011          ', ' + kind + ', True)'
     18012        lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
     18013        sub.call(['mv','lonlatProj.nc',lonlatProjfile])
     18014
     18015    oproj = NetCDFFile(lonlatProjfile, 'r')
     18016    if kind == 'lonlat':
     18017        olon = oproj.variables['lon']
     18018        olat = oproj.variables['lat']
     18019        lonmap = olon[:]
     18020        latmap = olat[:]
     18021        projlon = olon[:]
     18022        projlat = olat[:]
     18023        omapvals = oproj.variables['points']
     18024        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
     18025        Nmaplon = olon.shape[0]
     18026        Nmaplat = olat.shape[0]
     18027
     18028    elif kind == 'lonlat_dxdyFix':
    1800418029        oprojlon = oproj.variables['loncirc']
    1800518030        oprojlat = oproj.variables['latcirc']
    … …  
    1801518040
    1801618041#    Ntotvals = len(lonmap)
    18017     Nmaplon = Ntotvals
    18018     Nmaplat = Ntotvals
     18042        Nmaplon = Ntotvals
     18043        Nmaplat = Ntotvals
     18044    else:
     18045        print errormsg
     18046        print '  ' + fname + ": projection kind '" + kind + "' not ready!!"
     18047        quit(-1)
    1801918048
    1802018049# Matrix map creation
    … …  
    1806018089
    1806118090    print '  ' + fname + ': re-locating:',Ntotvals,'points...'
    18062     for iv in range(Ntotvals):
    18063 #    for iv in range(15):
    18064 #        print 'Lluis:',iv,'lon lat:',lonvs[iv],',',latvs[iv]
    18065         difflat = np.abs(projlat - latvs[iv])
    18066         mindiffL = np.min(difflat)
    18067         ilat = index_mat(difflat, mindiffL)
    18068 #        print '  Lluis mindiffL:',mindiffL,'<> ilat:',ilat,'lat_ilat:',projlat[ilat],'Nprojlon:',Nprojlon[ilat],'shape oporjlon:',oprojlon[ilat,:].shape, type(oprojlon[ilat,:])
    18069         loncirc = np.zeros((Nprojlon[ilat]), dtype=np.float)
    18070         loncirc[:] = np.asarray(oprojlon[ilat,0:int(Nprojlon[ilat])].flatten())
    18071         difflon = np.abs(loncirc - lonvs[iv])
    18072         mindiffl = difflon.min()
    18073         ilon = index_mat(difflon, mindiffl)
    18074 #        print '  difflon:',type(difflon),'shape difflon',difflon.shape,'shape loncirc:', loncirc.shape
    18075 #        print '  Lluis mindiffl:',mindiffl,'<> ilon:',ilon,'oprojlon:', loncirc[ilon]
    18076 
    18077         percendone(iv,Ntotvals,0.5,'done:')
    18078         if mindiffl > mindiff or mindiffL > mindiff:
    18079             print errormsg
    18080             print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
    18081               lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
    18082               'completed map closer than: ', mindiff, '!!'
    18083             print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
    18084             quit()
    18085 
    18086         if ilon > 0 and ilat > 0:
    18087 #            print 'Lluis shapes:',newvar.shape,'i lon/lat:',ilon,ilat
    18088             if Lmapvalsshape ==1:
    18089                 newvar[iv] = ovar[iv]
     18091    if kind == 'lonlat':
     18092        for iv in range(Ntotvals):
     18093            difflat = np.abs(projlat - latvs[iv])
     18094            mindiffL = np.min(difflat)
     18095            ilat = index_mat(difflat, mindiffL)
     18096            difflon = np.abs(projlon - lonvs[iv])
     18097            mindiffl = difflon.min()
     18098            ilon = index_mat(difflon, mindiffl)
     18099            percendone(iv,Ntotvals,0.5,'done:')
     18100            if mindiffl > mindiff or mindiffL > mindiff:
     18101                print errormsg
     18102                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
     18103                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
     18104                  'completed map closer than: ', mindiff, '!!'
     18105                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
     18106                quit()
     18107
     18108            if ilon > 0 and ilat > 0:
     18109                newvar[ilat,ilon] = ovar[iv]         
     18110                newnc.sync()
    1809018111            else:
    18091                 newvar[ilat,ilon] = ovar[iv]
     18112                print errormsg
     18113                print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
     18114                print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
     18115                quit(-1)
     18116
     18117        newnc.sync()
     18118                       
     18119
     18120    elif kind == 'lonlat_dxdyFix':
     18121
     18122        for iv in range(Ntotvals):
     18123#        for iv in range(15):
     18124#            print 'Lluis:',iv,'lon lat:',lonvs[iv],',',latvs[iv]
     18125            difflat = np.abs(projlat - latvs[iv])
     18126            mindiffL = np.min(difflat)
     18127            ilat = index_mat(difflat, mindiffL)
     18128#            print '  Lluis mindiffL:',mindiffL,'<> ilat:',ilat,'lat_ilat:',projlat[ilat],'Nprojlon:',Nprojlon[ilat],'shape oporjlon:',oprojlon[ilat,:].shape, type(oprojlon[ilat,:])
     18129            loncirc = np.zeros((Nprojlon[ilat]), dtype=np.float)
     18130            loncirc[:] = np.asarray(oprojlon[ilat,0:int(Nprojlon[ilat])].flatten())
     18131            difflon = np.abs(loncirc - lonvs[iv])
     18132            mindiffl = difflon.min()
     18133            ilon = index_mat(difflon, mindiffl)
     18134#            print '  difflon:',type(difflon),'shape difflon',difflon.shape,'shape loncirc:', loncirc.shape
     18135#            print '  Lluis mindiffl:',mindiffl,'<> ilon:',ilon,'oprojlon:', loncirc[ilon]
     18136
     18137            percendone(iv,Ntotvals,0.5,'done:')
     18138            if mindiffl > mindiff or mindiffL > mindiff:
     18139                print errormsg
     18140                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
     18141                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
     18142                  'completed map closer than: ', mindiff, '!!'
     18143                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
     18144                quit()
     18145
     18146            if ilon > 0 and ilat > 0:
     18147#                print 'Lluis shapes:',newvar.shape,'i lon/lat:',ilon,ilat
     18148                if Lmapvalsshape ==1:
     18149                    newvar[iv] = ovar[iv]
    1809218150           
    18093             newnc.sync()
    18094         else:
    18095             print errormsg
    18096             print '  ' + fname + ': point:',lonvs[iv],',',latvs[iv],' not relocated !!'
    18097             quit(-1)
    18098 
    18099     newnc.sync()
     18151                newnc.sync()
     18152            else:
     18153                print errormsg
     18154                print '  ' + fname + ': point:',lonvs[iv],',',latvs[iv],' not relocated !!'
     18155                quit(-1)
     18156
     18157        newnc.sync()
    1810018158# Global attributes
    1810118159##
    … …  
    1811818176    return
    1811918177
    18120 #Partialmap_Entiremap('longitude,latitude,std,5000.', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
     18178#Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat_dxdyFix', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
     18179#Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
    1812118180#quit()
    1812218181
Note: See TracChangeset for help on using the changeset viewer.