Changeset 719 in lmdz_wrf for trunk

Timestamp:

Apr 26, 2016, 4:52:53 PM (9 years ago)

Author:

lfita

Message:

Adding `CoarselonlatFindAll?': Function to search all values from a coarser version of the data

File:

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

trunk/tools/nc_var_tools.py

@@ -18721 +18721 @@
     return ilatlon, mindiffLl
 
+def CoarselonlatFindAll(onc, ilon, ilat, longvs, latvs, per, frac, out):
+    """ Function to search all values from a coarser version of the data
+      onc= netCDF object from an existing file
+      ilon, ilat= original 2D matrices with the longitudes and the latitudes
+      lonv, latv= longitude and latitude to find
+      per= period (as fraction over 1) of the fractions of the original grid to use to explore
+      frac= frequency of points at which syncronize file with calculations
+      out= True/False if outcome is required
+      >>> npt=100
+      >>> lonval0 = np.arange(0.,npt*1.,1.)
+      >>> latval0 = np.arange(0.,npt*1.,1.)
+      >>> lonval = np.zeros((npt,npt), dtype=np.float)
+      >>> latval = np.zeros((npt,npt), dtype=np.float)
+      >>> lonvalues = np.arange(0.25,10.,0.25)
+      >>> latvalues = np.arange(0.25,10.,0.25)
+
+      >>> for i in range(npt):
+      >>>     lonval[:,i] = lonval0[i]
+      >>>     latval[i,:] = latval0[i]
+      >>> CoarselonlatFindAll(onewnc, lonval, latval, lovalues, latvalues, .1, 100)
+      (array([5, 5]), 0.35355339059327379)
+    """
+    fname = 'CoarselonlatFindAll'
+
+    dx = ilon.shape[1]
+    dy = ilon.shape[0]
+
+    nlon = np.min(ilon)
+    xlon = np.max(ilon)
+    nlat = np.min(ilat)
+    xlat = np.max(ilat)
+
+    Nallpts = len(longvs)
+
+    fracx = int(dx*per)
+    fracy = int(dy*per)
+
+    fraclon = ilon[::fracy,::fracx]
+    fraclat = ilat[::fracy,::fracx]
+    fracdx = fraclon.shape[1]
+    fracdy = fraclon.shape[0]
+
+#    print 'fraclon _______'
+#    print fraclon
+
+#    print 'fraclat _______'
+#    print fraclat
+
+    if out:
+        ilatlon = np.zeros((2,Nallpts), dtype=int)
+        mindiffLl = np.zeros((Nallpts), dtype=np.float)
+
+    for iv in range(Nallpts):
+        lonv = longvs[iv]
+        latv = latvs[iv]
+
+        if lonv < nlon or lonv > xlon:
+            print errormsg
+            print '  ' + fname + ': longitude outside data range!!'
+            print '    given value:', lonv,'outside [',nlon,',',xlon,']'
+            quit(-1) 
+        if latv < nlat or latv > xlat:
+            print errormsg
+            print '  ' + fname + ': latitude outside data range!!'
+            print '    given value:', latv,'outside [',nlat,',',xlat,']'
+            quit(-1) 
+
+# Fraction point
+        difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
+        mindifffracLl = np.min(difffraclonlat)
+        ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)
+
+#    print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
+#    print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',             \
+#      fraclat[ilatlonfrac[0],ilatlonfrac[1]]
+#    print 'values lon, lat:', lonv, latv
+      
+# Providing fraction range
+        if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                        \
+          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
+#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
+#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
+            if ilatlonfrac[0] > 0: 
+                iybeg = (ilatlonfrac[0]-1)*fracy
+                iyend = ilatlonfrac[0]*fracy+1
+            else:
+                iybeg = 0
+                iyend = fracy+1
+            if ilatlonfrac[1] > 0: 
+                ixbeg = (ilatlonfrac[1]-1)*fracx
+                ixend = ilatlonfrac[1]*fracx+1
+            else:
+                ixbeg = 0
+                ixend = fracx+1
+        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
+          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
+#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
+#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
+            if ilatlonfrac[0] > 0: 
+                iybeg = (ilatlonfrac[0]-1)*fracy
+                iyend = ilatlonfrac[0]*fracy+1
+            else:
+                iybeg = 0
+                iyend = fracy+1
+            if ilatlonfrac[1] < fracdx:
+                if ilatlonfrac[1] != 0:
+                    ixbeg = (ilatlonfrac[1]-1)*fracx
+                    ixend = ilatlonfrac[1]*fracx+1
+                else:
+                    ixbeg = 0
+                    ixend = fracx+1
+            else:
+                ixbeg = fracdx*fracx
+                ixend = dx+1
+        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
+          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
+#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
+#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
+            if ilatlonfrac[0] < fracdy: 
+                if ilatlonfrac[0] != 0:
+                    iybeg = (ilatlonfrac[0]-1)*fracy
+                    iyend = ilatlonfrac[0]*fracy+1
+                else:
+                    iybeg = 0
+                    iyend = fracy+1
+            else:
+                iybeg = fracdy*fracy
+                iyend = dy+1
+            if ilatlonfrac[1] < fracdx: 
+                if ilatlonfrac[1] != 0:
+                    ixbeg = (ilatlonfrac[1]-1)*fracx
+                    ixend = ilatlonfrac[1]*fracx+1
+                else:
+                    ixbeg = 0
+                    ixend = fracx+1
+            else:
+                ixbeg = fracdx*fracx
+                ixend = dx+1
+        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                      \
+          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
+#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
+#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
+            if ilatlonfrac[0] < fracdy: 
+                if ilatlonfrac[0] != 0:
+                    iybeg = (ilatlonfrac[0]-1)*fracy
+                    iyend = ilatlonfrac[0]*fracy+1
+                else:
+                    iybeg = 0
+                    iyend = fracy+1
+            else:
+                iybeg = fracdy*fracy
+                iyend = dy+1
+            if ilatlonfrac[1] > 0: 
+                ixbeg = (ilatlonfrac[1]-1)*fracx
+                ixend = ilatlonfrac[1]*fracx+1
+            else:
+                ixbeg = 0
+                ixend = fracx+1
+    
+#        ibeg = [iybeg, ixbeg]
+#        iend = [iyend, ixend]
+#        print 'find window; beginning', ibeg, 'end:', iend
+        lon = ilon[iybeg:iyend,ixbeg:ixend]
+        lat = ilat[iybeg:iyend,ixbeg:ixend]
+    
+#        print 'lon _______'
+#        print lon
+#        print 'lat _______'
+#        print lat
+    
+# Find point
+        difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
+        if out:
+            mindiffLl[iv] = np.min(difflonlat)
+            ilatlon[:,iv] = index_mat(difflonlat, mindiffLl)
+            ilatlon[0,iv] = ilatlon[0,iv] + iybeg
+            ilatlon[1,iv] = ilatlon[1,iv] + ixbeg
+        else:
+            mindiffLl = np.min(difflonlat)
+            ilatlon = index_mat(difflonlat, mindiffLl)
+            ilatlon[0] = ilatlon[0] + iybeg
+            ilatlon[1] = ilatlon[1] + ixbeg
+
+#        print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
+#        print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
+#        quit()
+
+        if mindiffLl == 0. and type(newvar[ilatlon[0],ilatlon[1]]) == type(amsk): 
+            percendone(iv,Ninpts,0.5,'done:')
+            if mindiffLl > 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 difference:', mindiffLl
+                quit(-1)
+
+            if ilatlon[0] >= 0 and ilatlon[1] >= 0:
+                newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
+                newvarin[ilatlon[0],ilatlon[1]] = iv
+                newvarinpt[iv] = 1
+                newvarindiff[iv] = mindiffLl
+#                print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],                  \
+#                  'localized:', newvarinpt[iv], 'values:',                           \
+#                  newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],              \
+#                  'mindist:', newvarindiff[iv], 'point:',ilatlon    
+            else:
+                print errormsg
+                print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',',        \
+                  latvs[iv],' not relocated !!'
+                print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],               \
+                  'ilatlon:', ilatlon[1]
+                quit(-1)
+
+            if np.mod(iv,fracs) == 0: 
+                onc.sync()
+    if out:
+        return ilatlon, mindiffLl
+    else:
+        return
+
 def lonlatFind(ilon, ilat, lonv, latv):
     """ Function to search a given value from a lon, lat 2D data
@@ -19050 +19272 @@
     print '  ' + fname + ': from:', Ninpts,'re-locating:',Nptsf,'points...'
     if kind == 'Goode':
-        for iv0 in range(Nptsf):
-            iv = mapoints[ptsf[iv0]]
-            if newvarinpt[iv] == 0:
-
-#                difflonlat = np.sqrt((projlon-lonvs[iv])**2. + (projlat-latvs[iv])**2.)
-#                mindiffLl = np.min(difflonlat)
-#                ilatlon = index_mat(difflonlat, mindiffLl)
-                ilatlon, mindiffLl = CoarselonlatFind(projlon,projlat,lonvs[iv],latvs[iv],fracd)
-
-#                if ilatlon[0] != ilatlon2[0] or ilatlon[1] != ilatlon2[1]:
-#                    print ilatlon, '<', ilatlon2
-#                    print 'diff:', mindiffLl, '<', mindiffLl2
-#                    quit(-1)
-
-#                if mindiffLl != 0. or type(newvar[ilatlon[0],ilatlon[1]]) != type(amsk): 
-#                    print errormsg
-#                    if mindiffLl !=0.:
-#                        print '  ' + main + ': no zero', newvarindiff[iv], 'distance!!'
+        empty = CoarselonlatFind(newnc,projlon,projlat,lonvs[ptfs],latvs[ptfs],fracd,fracs,False)
+
+#        for iv0 in range(Nptsf):
+#            iv = mapoints[ptsf[iv0]]
+#            if newvarinpt[iv] == 0:
+
+##                difflonlat = np.sqrt((projlon-lonvs[iv])**2. + (projlat-latvs[iv])**2.)
+##                mindiffLl = np.min(difflonlat)
+##                ilatlon = index_mat(difflonlat, mindiffLl)
+#                ilatlon, mindiffLl = CoarselonlatFind(projlon,projlat,lonvs[iv],latvs[iv],fracd)
+
+##                if ilatlon[0] != ilatlon2[0] or ilatlon[1] != ilatlon2[1]:
+##                    print ilatlon, '<', ilatlon2
+##                    print 'diff:', mindiffLl, '<', mindiffLl2
+##                    quit(-1)
+
+##                if mindiffLl != 0. or type(newvar[ilatlon[0],ilatlon[1]]) != type(amsk): 
+##                    print errormsg
+##                    if mindiffLl !=0.:
+##                        print '  ' + main + ': no zero', newvarindiff[iv], 'distance!!'
+##                    else:
+##                        print '  ' + main+': point',projlon[ilatlon[0],ilatlon[1]],  \
+##                          ',', projlat[ilatlon[0],ilatlon[1]], 'already filled!!'
+##                        print '    value:',newvar[ilatlon[0],ilatlon[1]],'distance:',\
+##                          newvarindiff[iv]
+##
+###                    print '   iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
+###                    print '   mindiffLl:',mindiffLl,'ilatlon:',ilatlon
+###                    quit(-1)
+###            print '  Lluis; ' + fname + ' iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
+###            print '  Lluis; ' + fname + ' mindiffLl:',mindiffLl,'ilatlon:',ilatlon
+###            print '  Lluis; ' + fname + ' A found lon:',projlon[ilatlon[0], ilatlon[1]], 'lat:', projlat[ilatlon[0], ilatlon[1]]
+###            quit()
+##                else:
+#                if mindiffLl == 0. and type(newvar[ilatlon[0],ilatlon[1]]) == type(amsk): 
+#                    percendone(iv,Ninpts,0.5,'done:')
+#                    if mindiffLl > 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 difference:', mindiffLl
+#                        quit(-1)
+
+#                    if ilatlon[0] >= 0 and ilatlon[1] >= 0:
+#                        newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
+#                        newvarin[ilatlon[0],ilatlon[1]] = iv
+#                        newvarinpt[iv] = 1
+#                        newvarindiff[iv] = mindiffLl
+##                        print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],          \
+##                          'localized:', newvarinpt[iv], 'values:',                   \
+##                          newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],      \
+##                          'mindist:', newvarindiff[iv], 'point:',ilatlon    
 #                    else:
-#                        print '  ' + main+': point',projlon[ilatlon[0],ilatlon[1]],  \
-#                          ',', projlat[ilatlon[0],ilatlon[1]], 'already filled!!'
-#                        print '    value:',newvar[ilatlon[0],ilatlon[1]],'distance:',\
-#                          newvarindiff[iv]
-#
-##                    print '   iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
-##                    print '   mindiffLl:',mindiffLl,'ilatlon:',ilatlon
+#                        print errormsg
+#                        print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',',\
+#                          latvs[iv],' not relocated !!'
+#                        print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],       \
+#                          'ilatlon:', ilatlon[1]
+#                        quit(-1)
+
+#                if np.mod(iv,1000) == 0: 
+#                    newnc.sync()
+##                    print '  ' + fname + 'values localized:', newvar[:].compressed()
+##                if iv > 10:
+##                    newnc.sync()
+##                    newnc.close()
 ##                    quit(-1)
-##            print '  Lluis; ' + fname + ' iv:', iv, 'lon:', lonvs[iv],'lat:',latvs[iv]
-##            print '  Lluis; ' + fname + ' mindiffLl:',mindiffLl,'ilatlon:',ilatlon
-##            print '  Lluis; ' + fname + ' A found lon:',projlon[ilatlon[0], ilatlon[1]], 'lat:', projlat[ilatlon[0], ilatlon[1]]
-##            quit()
-#                else:
-                if mindiffLl == 0. and type(newvar[ilatlon[0],ilatlon[1]]) == type(amsk): 
-                    percendone(iv,Ninpts,0.5,'done:')
-                    if mindiffLl > 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 difference:', mindiffLl
-                        quit(-1)
-
-                    if ilatlon[0] >= 0 and ilatlon[1] >= 0:
-                        newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
-                        newvarin[ilatlon[0],ilatlon[1]] = iv
-                        newvarinpt[iv] = 1
-                        newvarindiff[iv] = mindiffLl
-#                        print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],          \
-#                          'localized:', newvarinpt[iv], 'values:',                   \
-#                          newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],      \
-#                          'mindist:', newvarindiff[iv], 'point:',ilatlon    
-                    else:
-                        print errormsg
-                        print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',',\
-                          latvs[iv],' not relocated !!'
-                        print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],       \
-                          'ilatlon:', ilatlon[1]
-                        quit(-1)
-
-                if np.mod(iv,1000) == 0: 
-                    newnc.sync()
-#                    print '  ' + fname + 'values localized:', newvar[:].compressed()
-#                if iv > 10:
-#                    newnc.sync()
-#                    newnc.close()
-#                    quit(-1)
 
     elif kind == 'lonlat':