Changeset 766 in lmdz_wrf

Timestamp:

May 10, 2016, 1:55:00 PM (9 years ago)

Author:

lfita

Message:

Almost working version with `lonlatfrac'

File:

• trunk/tools/nc_var_tools.py (modified) (8 diffs)

trunk/tools/nc_var_tools.py

@@ -13740 +13740 @@
 
 #    mindiff = np.min([minlat, minlon])/1.
-    mindiff = np.min([minlat, minlon])*10.e6
+    mindiffin = np.sqrt(minlat*minlat + minlon*minlon)
 
     print '  ' + fname + ': Closest latitude to Equator:', minlat
@@ -13856 +13856 @@
         quit(-1)
 
+    sortedlonmap = np.sort(lonmap)
+    sortedlatmap = np.sort(latmap)
+    Nsortlonm = len(sortedlonmap)
+    Nsortlatm = len(sortedlatmap)
+
+    diffsortlon = sortedlonmap[1:Nsortlonm-1]-sortedlonmap[0:Nsortlonm-2]
+    diffsortlat = sortedlatmap[1:Nsortlatm-1]-sortedlatmap[0:Nsortlatm-2]
+
+    mindifflonmap = np.min(diffsortlon)
+    maxdifflonmap = np.max(diffsortlon)
+    mindifflatmap = np.min(diffsortlat)
+    maxdifflatmap = np.max(diffsortlat)
+
+    mindiffmap = np.sqrt(mindifflonmap*mindifflonmap+mindifflatmap*mindifflatmap)
+    maxdiffmap = np.sqrt(maxdifflonmap*maxdifflonmap+maxdifflatmap*maxdifflatmap)
+
+    print '  ' + fname + ': minimum distance in target projection:', mindiffmap,     \
+      'maximum:', maxdiffmap
 # Matrix map creation
 ##
@@ -13902 +13920 @@
         basicvardef(newvarindiff, 'locindiff', 'distance between input point and ' + \
           'its final location','degree')
+        set_attributek(newvar,'authorized_minimum_distance',mindiff,'R')
 
 # map variable
@@ -13923 +13942 @@
             newvarin = newnc.createVariable('inpts', 'i4', ('pts'))    
         else:
-            newvarin = newnc.createVariable('inpts', 'i4', ('lat','lon'))
-
-        basicvardef(newvarin, 'inpts', 'Equivalent point from the input source', '-')
+            if kind[Lkind-4:Lkind] == 'frac':
+                Nmaxcoinpt = np.int((maxdiffmap/mindiffin)*(maxdiffmap/mindiffin))
+                print '  ' + fname + ': maximum theoretical number of coincident '   + \
+                  'grid points in the same target gid point:', Nmaxcoinpt
+                newdim = newnc.createDimension('Ncoinpt', Nmaxcoinpt)
+                newvarin = newnc.createVariable('inpts', 'i4', ('Ncoinpt','lat',     \
+                  'lon'))
+                basicvardef(newvarin, 'inpts', 'Equivalent points from the input '   + 
+                  'source', '-')
+            else:
+                newvarin = newnc.createVariable('inpts', 'i4', ('lat','lon'))
+                basicvardef(newvarin, 'inpts', 'Equivalent point from the input ' +  \
+                  'source', '-')
+
         newvar.setncattr('coordinates', 'lon lat')
 
@@ -13950 +13980 @@
 
 # Reducing the searching points
-    Nsearchpt = 499999
+    Nsearchpt = 999999
     Ninpts=np.min([Nsearchpt, Ninpts-ipoint])
     newvarinvals = newvarinpt[ipoint:ipoint+Ninpts]
@@ -13957 +13987 @@
     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]]
-    print '  ' + fname + ': from:', Ninpts,'re-locating:',Nptsf,'points starting at',\
-      ipoint,'...'
 
     Nptsf = len(ptsf)
@@ -14002 +14027 @@
             idiff, ilonlatv = fin.module_forinterpolate.interpolate1dll(projlon,     \
               projlat, lonvss, latvss, np.float64(mindiff), inptss)
-            print 'Lluis shapes: newvar', newvar.shape,'ilonlatv:',ilonlatv.shape,'range:',np.min([len(idiff),fracs])
-            print 'Lluis shapes: ovar:',ovar.shape
 
             for i in range(np.min([len(idiff),fracs])):
@@ -14030 +14053 @@
 
             for i in range(np.min([len(idiff),fracs])):
-                print i,'Lluis:',ovar[iri+i], ilonlatv[i,1],',',ilonlatv[i,0],newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]], type(newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]])
-                if newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]].mask: 
-                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]].mask = False
-                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] = 1.
+                if type(newvar[ovar[iri+i]-1,ilonlatv[i,1],ilonlatv[i,0]]) ==        \
+                  type(ma.array(1)) and newvar[ovar[iri+i]-1,ilonlatv[i,1],          \
+                  ilonlatv[i,0]].mask: 
+                    newvar[ovar[iri+i]-1,ilonlatv[i,1],ilonlatv[i,0]].mask = False
+                    newvar[ovar[iri+i]-1,ilonlatv[i,1],ilonlatv[i,0]] = 1.
+                    newvarin[0,ilonlatv[i,1],ilonlatv[i,0]] = iri + i
                 else:
-                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] =                \
-                      newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] + 1
-                newvarin[ilonlatv[i,1],ilonlatv[i,0]] = iri + i
+                    Npt = np.sum(newvarin[:,ilonlatv[i,1],ilonlatv[i,0]] > 0)
+                    print ilonlatv[i,1],ilonlatv[i,0],'Npt:',Npt
+                    print newvarin[:,ilonlatv[i,1],ilonlatv[i,0]]
+                    newvar[ovar[iri+i]-1,ilonlatv[i,1],ilonlatv[i,0]] =                \
+                      newvar[ovar[iri+i]-1,ilonlatv[i,1],ilonlatv[i,0]] + 1.
+                    newvarin[Npt+1,ilonlatv[i,1],ilonlatv[i,0]] = iri + i
                 newvarinpt[iri+i] = 1
                 newvarindiff[iri+i] = idiff[i]