Changeset 714 in lmdz_wrf for trunk

Timestamp:

Apr 25, 2016, 6:56:05 PM (9 years ago)

Author:

lfita

Message:

Adding `Partialmap_EntiremapFor': Version using Fortran subroutines

File:

• trunk/tools/nc_var_tools.py (modified) (1 diff)

trunk/tools/nc_var_tools.py

@@ -19203 +19203 @@
     return 
 
+def Partialmap_EntiremapFor(values, filen, varn):
+    """ Function to transform from a partial global map (e.g.: only land points) to an entire one usinf Fortran code
+      Coincidence of points is done throughout a first guess from fractions of the total domain of search
+      Using fortran codes: module_ForInterpolate.F90, module_generic.F90
+      f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
+      values= [lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile],[fracd]
+        [lonname]: name of the longitude variable
+        [latname]: name of the latitude variable
+        [fillVal]: value for '_FillValue', 'std' for the standard value
+           Float = 1.e20
+           Character = '-'
+           Integer = -99999
+           Float64 = 1.e20
+           Integer32 = -99999
+        [resolution]: resolution of the map
+        [kind]: kind of target projection
+         'lonlat': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator
+         'lonlat_dxdyFix': projection with a fixed grid distance
+         'Goode': Goode projection
+        [lonlatProjfile]: file with the lon,lat of the desired projection. 'None' to be computed and written on fly
+        [fracd]: Percentage of the fractions within perform the first guess search
+      filen= name of the netCDF file
+      varn= name of the variable
+    """
+    import numpy.ma as ma
+    import subprocess as sub
+
+    fname = 'Partialmap_Entiremap'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print Partialmap_Entiremap.__doc__
+        quit()
+
+    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind],' +               \
+      '[lonlatProjfile],[fracd]'
+    check_arguments(fname, values, arguments, ',')
+
+    ofile = 'EntireGlobalMap.nc'
+
+    onc = NetCDFFile(filen, 'r')
+    if not onc.variables.has_key(varn):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "variable '" + varn + "' !!"
+        quit(-1)
+
+    lonname = values.split(',')[0]
+    latname = values.split(',')[1]
+    fval0 = values.split(',')[2]
+    resolution = np.float(values.split(',')[3])
+    kind = values.split(',')[4]
+    Projfile = values.split(',')[5]
+    fracd = np.float(values.split(',')[6])
+
+    if Projfile == 'None':
+        lonlatProjfile = None
+    else:
+        lonlatProjfile = Projfile
+
+    if not onc.variables.has_key(lonname):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "longitude '" + lonname + "' !!"
+        quit(-1)
+    if not onc.variables.has_key(latname):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have " +       \
+          "latitude '" + latname + "' !!"
+        quit(-1)
+
+    olon = onc.variables[lonname]
+    olat = onc.variables[latname]
+
+    lonvs = olon[:]
+    latvs = olat[:]
+
+    Ninpts = len(lonvs)
+
+    nlat = np.min(latvs)
+    nlon = np.min(lonvs)
+
+    minlat = np.min(np.abs(latvs))
+    minlon = np.min(np.abs(lonvs))
+
+#    mindiff = np.min([minlat, minlon])/1.
+    mindiff = np.min([minlat, minlon])*10.e6
+
+    print '  ' + fname + ': Closest latitude to Equator:', minlat
+    print '  ' + fname + ': Closest longitude to Greenwich Meridian:', minlon
+    print '  ' + fname + ': Minium distance for point coincidence:', mindiff
+
+# Existing longitudes along/closest to the Equator
+    eqlons = []
+    for i in range(len(latvs)):
+        if latvs[i] == minlat: eqlons.append(lonvs[i])
+
+# Existing latitudes along/closest to the Greenwich Meridian
+    grlats = []
+    for i in range(len(lonvs)):
+        if lonvs[i] == minlon: grlats.append(latvs[i])
+
+    sorteqlons = np.sort(eqlons)
+    sortgrlats = np.sort(grlats)
+ 
+    Neqlons = len(sorteqlons)
+    Ngrlats = len(sortgrlats)
+
+    if Neqlons > 1:
+        diffeqlons = sorteqlons[1:Neqlons-1] - sorteqlons[0:Neqlons-2]     
+        mindeqlon = np.min(diffeqlons)
+        print 'N sorteqlons:',Neqlons,'min deqlon:', mindeqlon
+    else:
+        mindeqlon = None
+
+    if Ngrlats > 1:
+        mindgrlat = np.min(diffgrlats)
+        diffgrlats = sortgrlats[1:Ngrlats-1] - sortgrlats[0:Ngrlats-2]
+        print 'N sortgrlats:',Ngrlats,'min dgrlat:', mindgrlat
+        latmap = np.range(0,360.,360./mindeqlon)
+    else:
+        mindgrlat = None
+
+# Fixing in case it has not worked
+    if mindeqlon is not None and mindgrlat is None: mindgrlat = mindeqlon
+    if mindeqlon is None and mindgrlat is not None: mindeqlon = mindgrlat
+    if mindeqlon is None and mindgrlat is None:
+        print errormsg
+        print fname + ': Not enough values along the Equator and Greenwich!!'
+        quit(-1)
+
+    if lonlatProjfile is None:
+        lonlatProjfile = kind + '.nc'
+        print warnmsg
+        print '  ' + fname + ": no reference map !!"
+        print "    creation of '" + lonlatProjfile + "'"
+        print '    creating it via: lonlatProj(', resolution, ',', resolution,       \
+          ', ' + kind + ', True)'
+        lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
+        sub.call(['mv','lonlatProj.nc',lonlatProjfile])
+
+    oproj = NetCDFFile(lonlatProjfile, 'r')
+    if kind == 'Goode':
+        olon = oproj.variables['lon']
+        olat = oproj.variables['lat']
+        lonmap = olon[:]
+        latmap = olat[:]
+        projlon = olon[:]
+        projlat = olat[:]
+        omapvals = oproj.variables['ptid']
+        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
+        Nmaplon = olon.shape[1]
+        Nmaplat = olat.shape[0]
+
+    elif kind == 'lonlat':
+        olon = oproj.variables['lon']
+        olat = oproj.variables['lat']
+        lonmap = olon[:]
+        latmap = olat[:]
+        projlon = olon[:]
+        projlat = olat[:]
+        omapvals = oproj.variables['points']
+        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
+        Nmaplon = olon.shape[0]
+        Nmaplat = olat.shape[0]
+
+    elif kind == 'lonlat_dxdyFix':
+        oprojlon = oproj.variables['loncirc']
+        oprojlat = oproj.variables['latcirc']
+        oNprojlon = oproj.variables['Nloncirc']
+        projlat = oprojlat[:]
+        Nprojlon = oNprojlon[:]
+        Ntotvals = len(oproj.dimensions['Npts'])
+        olon = oproj.variables['lon']
+        olat = oproj.variables['lat']
+        lonmap = olon[:]
+        latmap = olat[:]
+        omapvals = oproj.variables['points']
+
+#    Ntotvals = len(lonmap)
+        Nmaplon = Ntotvals
+        Nmaplat = Ntotvals
+    else:
+        print errormsg
+        print '  ' + fname + ": projection kind '" + kind + "' not ready!!"
+        quit(-1)
+
+# Matrix map creation
+##
+    ovar = onc.variables[varn]
+    vartype = type(ovar[0])
+    varlongname = ovar.long_name
+    varunits = ovar.units
+
+    fval = fillvalue_kind(type(ovar[0]), fval0)
+
+# Final map creation
+##
+    if not os.path.isfile(ofile):
+        print '  ' + fname + "File '" + ofile + "' does not exist !!"
+        print '    cration of output file'
+        newnc = NetCDFFile(ofile, 'w')
+# dimensions
+        newdim = newnc.createDimension('lon',Nmaplon)
+        newdim = newnc.createDimension('lat',Nmaplat)
+        newdim = newnc.createDimension('inpts', Ninpts)
+        newdim = newnc.createDimension('pts', Ntotvals)
+
+# Variables
+        if kind == 'Goode':
+            newvar = newnc.createVariable('lon','f8',('lat', 'lon'))
+        else:
+            newvar = newnc.createVariable('lon','f8',('lon'))
+        basicvardef(newvar, 'lon', 'Longitudes','degrees_East')
+        newvar[:] = lonmap
+        newvar.setncattr('axis', 'X')
+        newvar.setncattr('_CoordinateAxisType', 'Lon')   
+
+        if kind == 'Goode':
+            newvar = newnc.createVariable('lat','f8',('lat','lon'))
+        else:
+            newvar = newnc.createVariable('lat','f8',('lat'))
+        basicvardef(newvar, 'lat', 'Latitudes','degrees_North')
+        newvar[:] = latmap
+        newvar.setncattr('axis', 'Y')
+        newvar.setncattr('_CoordinateAxisType', 'Lat')   
+
+        newvarinpt = newnc.createVariable('locinpt','i',('inpts'))
+        basicvardef(newvarinpt, 'locinpt', 'input point located: 0: no, 1: yes','-')
+
+        newvarindiff = newnc.createVariable('locindiff','f4',('inpts'))
+        basicvardef(newvarindiff, 'locindiff', 'distance between input point and its final location','degree')
+
+# map variable
+        Lmapvalsshape = len(omapvals.shape)
+        if Lmapvalsshape == 1:
+            newvar = newnc.createVariable(varn, nctype(vartype), ('pts'), fill_value=fval)
+        else:
+            newvar = newnc.createVariable(varn, nctype(vartype), ('lat','lon'), fill_value=fval)
+
+        basicvardef(newvar, varn, varlongname, varunits)
+        newvar.setncattr('coordinates', 'lon lat')
+
+        if Lmapvalsshape == 1:
+            newvarin = newnc.createVariable('inpts', 'i4', ('pts'))    
+        else:
+            newvarin = newnc.createVariable('inpts', 'i4', ('lat','lon'))
+
+        basicvardef(newvarin, 'inpts', 'Equivalent point from the input source', '-')
+        newvar.setncattr('coordinates', 'lon lat')
+
+    else:
+        print '  ' + fname + "File '" + ofile + "' exists !!"
+        print '    reading values from file'
+        newnc = NetCDFFile(ofile, 'a')
+        newvar = newnc.variables[varn]
+        newvarin = newnc.variables['inpts']
+        newvarinpt = newnc.variables['locinpt']
+        newvarindiff = newnc.variables['locindiff']
+
+    amsk = np.arange(3)
+    amsk = ma.masked_equal(amsk, 0)
+
+#    fraclon = projlon[::Nmaplon*0.1]
+#    fraclat = projlat[::Nmaplat*0.1]
+
+#    print 'fraclon________________', fraclon.shape
+#    print fraclon
+
+#    print 'fraclat________________', fraclat.shape
+#    print fraclat
+
+# Reducing the searching points
+    newvarinvals = newvarinpt[:]
+    maskpt = np.where(newvarinvals.mask == True, False, True)
+    points = np.arange(Ninpts)
+    mapoints = ma.array(points, mask=maskpt)
+    ptsf = mapoints.compressed()
+
+    Nptsf = len(ptsf)
+    print Ninpts,'Npoints to find:', len(ptsf), ptsf[0:10], newvarindiff[ptsf[0:10]]
+# Error at 150024, 150025, 151709, 153421
+    print '  ' + fname + ': from:', Ninpts,'re-locating:',Nptsf,'points...'
+    fracs = 1000
+    if kind == 'Goode':
+#        newvar,newvarin,newvarinpt,newvarindiff =                                    \
+#          fin.module_forinterpolate.coarseinterpolate(projlon, projlat, lonvs,       \
+#          latvs, percen, mindiff, ivar, dimx, dimy, ninpts)
+        for ir in range(0,Ninpts,fracs):
+            iri = ir
+            ire = ir + fracs + 1
+            print iri,',',ire
+            newvar,newvarin,newvarinpt,newvarindiff =                                \
+              fin.module_forinterpolate.coarseinterpolate(projlon, projlat,          \
+              lonvs[iri:ire], latvs[iri:ire], fracd, mindiff, ivar, Nmaplon,         \
+              Nmaplat, fracs)
+            newnc.sync()
+
+    elif kind == 'lonlat':
+        for iv in range(Ntotvals):
+            difflat = np.abs(projlat - latvs[iv])
+            mindiffL = np.min(difflat)
+            ilat = index_mat(difflat, mindiffL)
+            difflon = np.abs(projlon - lonvs[iv])
+            mindiffl = difflon.min()
+            ilon = index_mat(difflon, mindiffl)
+            percendone(iv,Ntotvals,0.5,'done:')
+            if mindiffl > mindiff or mindiffL > mindiff:
+                print errormsg
+                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
+                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
+                  'completed map closer than: ', mindiff, '!!'
+                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
+                quit()
+
+            if ilon >= 0 and ilat >= 0:
+                newvar[ilat,ilon] = ovar[iv]         
+                newnc.sync()
+            else:
+                print errormsg
+                print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
+                print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
+                quit(-1)
+
+        newnc.sync()
+
+    elif kind == 'lonlat_dxdyFix':
+
+        for iv in range(Ntotvals):
+#        for iv in range(15):
+#            print 'Lluis:',iv,'lon lat:',lonvs[iv],',',latvs[iv]
+            difflat = np.abs(projlat - latvs[iv])
+            mindiffL = np.min(difflat)
+            ilat = index_mat(difflat, mindiffL)
+#            print '  Lluis mindiffL:',mindiffL,'<> ilat:',ilat,'lat_ilat:',projlat[ilat],'Nprojlon:',Nprojlon[ilat],'shape oporjlon:',oprojlon[ilat,:].shape, type(oprojlon[ilat,:])
+            loncirc = np.zeros((Nprojlon[ilat]), dtype=np.float)
+            loncirc[:] = np.asarray(oprojlon[ilat,0:int(Nprojlon[ilat])].flatten())
+            difflon = np.abs(loncirc - lonvs[iv])
+            mindiffl = difflon.min()
+            ilon = index_mat(difflon, mindiffl)
+#            print '  difflon:',type(difflon),'shape difflon',difflon.shape,'shape loncirc:', loncirc.shape
+#            print '  Lluis mindiffl:',mindiffl,'<> ilon:',ilon,'oprojlon:', loncirc[ilon]
+
+            percendone(iv,Ntotvals,0.5,'done:')
+            if mindiffl > mindiff or mindiffL > mindiff:
+                print errormsg
+                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
+                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
+                  'completed map closer than: ', mindiff, '!!'
+                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
+                quit()
+
+            if ilon >= 0 and ilat >= 0:
+                if Lmapvalsshape ==1:
+                    newvar[iv] = ovar[iv]            
+                newnc.sync()
+            else:
+                print errormsg
+                print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
+                print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
+                quit(-1)
+
+        newnc.sync()
+# Global attributes
+##
+    newnc.setncattr('script',  fname)
+    newnc.setncattr('version',  '1.0')
+    newnc.setncattr('author',  'L. Fita')
+    newnc.setncattr('institution',  'Laboratoire de Meteorology Dynamique')
+    newnc.setncattr('university',  'Pierre et Marie Curie')
+    newnc.setncattr('country',  'France')
+    for attrs in onc.ncattrs():
+        attrv = onc.getncattr(attrs)
+        attr = set_attribute(newnc, attrs, attrv)
+
+    newnc.sync()
+    onc.close()
+    newnc.close()
+
+    print fname + ": Successfull written of file: '" + ofile + "' !!"
+
+    return 
+
 #Partialmap_Entiremap('longitude,latitude,std,5000.,Goode,Goode_5km.nc', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/#carteveg5km.nc', 'vegetation_map')
 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat_dxdyFix', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')