Changeset 409 in lmdz_wrf for trunk

Timestamp:

May 4, 2015, 6:58:59 PM (10 years ago)

Author:

lfita

Message:

Adding 'remapnn' Function to remap to the nearest neightbor a variable using projection from another file

Location:

trunk/tools

Files:

• nc_var.py (modified) (2 diffs)
• nc_var_tools.py (modified) (3 diffs)

trunk/tools/nc_var.py

@@ -29 +29 @@
   'maskvar',                                                                         \
   'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
-  'compute_opersvarsfiles',                                                          \
+  'remapnn',                                                                         \
   'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'setvar_asciivalues',               \
   'sorttimesmat', 'spacemean', 'statcompare_files', 'submns', 'subyrs', 'TimeInf',   \
@@ -189 +189 @@
 elif oper == 'opersvarsfiles':
     ncvar.compute_opersvarsfiles(opts.values, opts.varname)
+elif oper == 'remapnn':
+    ncvar.remapnn(opts.values, opts.ncfile, opts.varname)
 elif oper == 'seasmean':
     ncvar.seasmean(timename, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -7334 +7334 @@
             self.quantiles = mask_quantiles(finalvalues, Nq) 
 
-def remapnn(values, filename, varn, ddegL, ddegl):
+def remapnn(values, filename, varname):
+    """ Function to remap to the nearest neightbor a variable using projection from another file
+    values=[newprojectionfile]|[newlonname,newlatname]|[newlondim,newlatdim]|[oldlonname,oldlatname]|
+     [oldlondim,oldlatdim]
+      [newprojectionfile]: name of the file with the new projection
+      [newlonname,newlatname]: name of the longitude and the latitude variables in the new file
+      [newlondim,newlatdim]: name of the dimensions for the longitude and the latitude in the new file
+      [oldlonname,oldlatname]: name of the longitude and the latitude variables in the old file
+      [oldlondim,oldlatdim]: name of the dimensions for the longitude and the latitude in the old file
+    filename= netCDF file name
+    varn= variable name ('all', for all variables)
+    """ 
+    fname = 'remapnn'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print remapnn.__doc__
+        quit()
+
+    expectargs = '[newprojectionfile]|[newlonname,newlatname]|[newlondim,' +         \
+      'newlatdim]|[oldlonname,oldlatname]|[oldlondim,oldlatdim]'
+ 
+    check_arguments(fname,len(expectargs.split('|')),values,'|',expectargs)
+
+    newprojectionfile = values.split('|')[0]
+    newlonlatname = values.split('|')[1]
+    newlonlatdim = values.split('|')[2].split(',')
+    oldlonlatname = values.split('|')[3]
+    oldlonlatdim = values.split('|')[4].split(',')
+
+    if varname != 'all':
+        ofile = 'remapnn_' + varname + '.nc'
+    else:
+        ofile = 'remapnn_' + filename
+
+    filexist(filename, errormsg, 'old file')
+    filexist(newprojectionfile, errormsg, 'new file')
+
+    oldlon=oldlonlatname.split(',')[0]
+    oldlat=oldlonlatname.split(',')[1]
+
+    oldncfile = NetCDFFile(filename,'r')
+    varinfile(oldncfile, filename, errormsg, 'old', oldlon)
+    varinfile(oldncfile, filename, errormsg, 'old', oldlat)
+
+    newlon=newlonlatname.split(',')[0]
+    newlat=newlonlatname.split(',')[1]
+
+    newncfile = NetCDFFile(newprojectionfile,'r')
+    varinfile(newncfile, newprojectionfile, errormsg, 'old', newlon)
+    varinfile(newncfile, newprojectionfile, errormsg, 'old', newlat)
+
+    oldlon2d, oldlat2d = lonlat2D(oldncfile, oldlon, oldlat)
+    newlon2d, newlat2d = lonlat2D(newncfile, newlon, newlat)
+
+    newdimx = newlon2d.shape[1]
+    newdimy = newlon2d.shape[0]
+
+    olddimx = oldlon2d.shape[1]
+    olddimy = oldlon2d.shape[0]
+
+# We are going to do the search in two phases
+#   1.- Searching on grid points every 10% of the size of the origina domain
+#   2.- Within the correspondant box of 20%x20% percent
+
+# 1.- first pairs old-->new
+##
+    dx10 = int(olddimx/10)
+    dy10 = int(olddimy/10)
+
+    oldlon2d10 = oldlon2d[slice(dy10,olddimy,dy10), slice(dx10,olddimx,dx10)]
+    oldlat2d10 = oldlat2d[slice(dy10,olddimy,dy10), slice(dx10,olddimx,dx10)]
+
+    firstmappairs = np.ones((2, newdimy, newdimx), dtype=int)
+    print '     ' + fname + ' first searching nearest neighbor by',dx10,',',dy10,'...'
+    for jnew in range(newdimy):
+        for inew in range(newdimx):
+            dist=np.sqrt((newlon2d[jnew, inew] - oldlon2d10)**2. +                   \
+              (newlat2d[jnew, inew] - oldlat2d10)**2.)
+            mindist = np.min(dist)
+            minydist10, minxdist10 = index_mat(dist,mindist)
+            firstmappairs[:,jnew,inew] = [int(dy10*(1+minydist10)),                  \
+              int(dx10*(1+minxdist10))]
+
+# 2.- Looking around the closest coarse point
+##
+    mappairs = np.ones((2, newdimy, newdimx), dtype=int)
+    dx2 = dx10/2
+    dy2 = dy10/2
+
+    print '     ' + fname + ' searching nearest neighbor...'
+    for jnew in range(newdimy):
+        for inew in range(newdimx):
+            i10 = firstmappairs[1,jnew,inew]
+            j10 = firstmappairs[0,jnew,inew]
+
+            oldlon2d2 = oldlon2d[slice(j10-dy2,j10+dy2), slice(i10-dx2,i10+dx2)]
+            oldlat2d2 = oldlat2d[slice(j10-dy2,j10+dy2), slice(i10-dx2,i10+dx2)]
+
+            dist=np.sqrt((newlon2d[jnew, inew] - oldlon2d2)**2. +                    \
+              (newlat2d[jnew, inew] - oldlat2d2)**2.)
+#            print 'oldlon2d2________________',i10,j10
+#            print oldlon2d2
+#            print 'oldlat2d2________________'
+#            print oldlat2d2
+#            print 'dist________________'
+#            print dist
+
+            mindist = np.min(dist)
+            minydist, minxdist = index_mat(dist,mindist)
+            mappairs[:,jnew,inew] = [int(j10+minydist), int(i10+minxdist)]
+#            print jnew,inew,'|',firstmappairs[:,jnew,inew],':',minxdist,minydist,    \
+#              'final:',mappairs[:,jnew,inew],'lon,lat old:',oldlon2d[minydist,minxdist],\
+#              oldlat2d[minydist,minxdist], 'new:', newlon2d[jnew, inew],             \
+#              newlat2d[jnew, inew],'mindist:',mindist
+
+
+#            quit()
+
+# Creation of the new file
+##
+    ncfnew = NetCDFFile(ofile,'w')
+   
+    for dn in oldncfile.dimensions:
+        if dn == oldlon:
+            ncfnew.createDimension(newlonlatdim[0], newdimx)
+        elif dn == oldlat:
+            ncfnew.createDimension(newlonlatdim[1], newdimy)
+        else:
+            if oldncfile.dimensions[dn].isunlimited():
+                ncfnew.createDimension(dn, None)
+            else:
+                ncfnew.createDimension(dn, len(oldncfile.dimensions[dn]))
+
+    ncfnew.sync()
+
+# looping along all the variables
+##
+    if varname == 'all':
+        oldvarnames = oldncfile.variables
+    else:
+        oldvarnames = [varname]
+
+    for varn in oldvarnames:
+        print '  adding:',varn,'...'
+        ovar = oldncfile.variables[varn]
+        vardims = ovar.dimensions
+        vartype = ovar.dtype
+        varattrs = ovar.ncattrs()
+
+#       Looking if variable has lon and/or lat:
+        varlonlat = False
+        if searchInlist(vardims, oldlonlatdim[0]): varlonlat = True
+        if searchInlist(vardims, oldlonlatdim[1]): varlonlat = True
+
+        if varn != oldlon and varn != oldlat and varlonlat:
+            newvarshape = []
+            newvardimn = []
+            newslice = []
+            oldslice = []
+            for dn in vardims:
+                if dn == oldlat:
+                    newvarshape.append(newdimy)
+                    newvardimn.append(newlonlatdim[1])
+                elif dn == oldlon:
+                    newvarshape.append(newdimx)
+                    newvardimn.append(newlonlatdim[0])
+                else:
+                    newvarshape.append(len(oldncfile.dimensions[dn]))
+                    newvardimn.append(dn)
+                    newslice.append(slice(0,len(oldncfile.dimensions[dn])))
+                    oldslice.append(slice(0,len(oldncfile.dimensions[dn])))
+
+            newvarval = np.zeros(tuple(newvarshape), dtype=vartype)
+
+            for jnew in range(newdimy):
+                for inew in range(newdimx):
+                    oldslice0 = oldslice[:]
+                    newslice0 = newslice[:]
+
+                    if searchInlist(vardims,oldlat):
+                        oldslice0.append(mappairs[0,jnew,inew])
+                        newslice0.append(jnew)
+
+                    if searchInlist(vardims,oldlon):
+                        oldslice0.append(mappairs[1,jnew,inew])
+                        newslice0.append(inew)
+
+#                    if ((jnew*newdimx + inew)%100 == 0): print '      ' + fname + ': ',  \
+#                      '{0:5.2f}'.format(float((jnew*newdimx + inew)*100./(newdimx*newdimy))),\
+#                       '% done'
+
+                    newvarval[tuple(newslice0)] = ovar[tuple(oldslice0)]
+
+        elif varn == oldlon:
+            ovar = newncfile.variables[newlon]
+            vardims = ovar.dimensions
+            vartype = ovar.dtype
+            varattrs = ovar.ncattrs()
+
+            newvardimn = (newlonlatdim[1], newlonlatdim[0])
+            newvarval = newlon2d
+
+        elif varn == oldlat:
+            ovar = newncfile.variables[newlat]
+            vardims = ovar.dimensions
+            vartype = ovar.dtype
+            varattrs = ovar.ncattrs()
+
+            newvardimn = (newlonlatdim[1], newlonlatdim[0])
+            newvarval = newlat2d
+
+        else:
+            newvardimn = vardims
+            newvarval = ovar[:]
+
+        if searchInlist(varattrs, '_FillValue'):
+            fvalv = ovar.getncattr('_FillValue')
+            newvar = ncfnew.createVariable(varn, vartype, tuple(newvardimn),         \
+              fill_value=fvalv)
+        else:
+            newvar = ncfnew.createVariable(varn, vartype, tuple(newvardimn))
+
+        for ncattr in varattrs:
+            if ncattr != '_FillValue':
+                attrv = ovar.getncattr(ncattr)
+                newattr = set_attribute(newvar,ncattr,attrv)
+                
+        attr = newvar.setncattr('mask', 'variable remapped at nearest neighbors ' +  \
+          'using ' + values.split(':')[0] + ' using as longitude and latitude values ')
+
+        newvar[:] = newvarval
+
+
+    ncfnew.sync()
+    ncfnew.close()
+
+    newncfile.close()
+    oldncfile.close()
+
+# Adding attributes
+    fgaddattr(filename, ofile)
+
+    print 'remapnn: new file "' + ofile + '" with remapped variable written !!'
+
+#remapnn('met_em.d01.1979-01-01_00:00:00.nc|XLONG_M,XLAT_M|east_weast,south_north|longitude,latitude|longitude,latitude', \
+#  'Albedo.nc', 'all')
+#quit()
+
+def remapnn_old(values, filename, varn, ddegL, ddegl):
     """ Function to remap to the nearest neightbor a variable using projection from another file
     values=[newprojectionfile]:[newlonname,newlatname]:[oldlonname,oldlatname]
@@ -7381 +7630 @@
 # pairs old-->new
 ##
-    mappairs = np.ones((newdimy, newdimx, 2), dtype=int)
+    mappairs = np.ones((2, newdimy, newdimx), dtype=int)
     mappairs = mappairs*fillValue
     print '     ' + fname + ' searching nearest neighbor...'
@@ -7420 +7669 @@
                                mindist = dist
 ##                               print jold, iold, 'dist: ', dist, 'L', maskindlon[iold,:], 'l', maskindlat[jold,:]
-                               mappairs[jnew, inew, 0] = maskindlon[iold,1]
-                               mappairs[jnew, inew, 1] = maskindlat[jold,0]
+                               mappairs[0, jnew, inew] = maskindlon[iold,1]
+                               mappairs[1, jnew, inew] = maskindlat[jold,0]
 
 ##                quit()