Changeset 2601 in lmdz_wrf for trunk/tools/generic_tools.py

Timestamp:

Jun 12, 2019, 8:40:32 PM (5 years ago)

Author:

lfita

Message:

Adding:

• `point_georef': Function to georeference the matrix coordinates of a given point using a couple of 2D longitude,latitude fields

File:

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

trunk/tools/generic_tools.py

@@ -191 +191 @@
 #  [YYYY][MM][DD][HH][MI][SS] format) in a 360 years calendar
 # points_2Dline: Function to define a line from a pair of points on a plane
+# point_georef: Function to georeference the matrix coordinates of a given point using a 
+#   couple of 2D longitude,latitude fields
 # PolyArea: Function to compute the area of the polygon following 'Shoelace formula'
 # pretty_int: Function to plot nice intervals
@@ -16768 +16770 @@
     return values
 
+def point_georef(lonv, latv, ivalue, jvalue, kind):
+    """ Function to georeference the matrix coordinates of a given point using a 
+        couple of 2D longitude,latitude fields
+      lonv: 2D longitude values
+      latv: 2D latitude values
+      ivalue: x-coordinate of the point to interpolate
+      jvalue: y-coordinate of the point to interpolate
+      kind: kind of georefernce methodology
+        - 'nearest': using the nearest values
+        - '4invdist': using the inverted distance of the 4 closest ones
+      >>> lon1D = [5., 10., 15.]
+      >>> lat1D = [-25., -30., -35.]
+      >>> lon, lat = np.meshgrid(lon1D, lat1D)
+      >>> ipt = 1.25
+      >>> jpt = 0.75
+      >>> point_georef(lon, lat, ipt, jpt, 'nearest')
+      (10.0, -30.0)
+      >>> point_georef(lon, lat, ipt, jpt, '4invdist')
+      (11.751864530113494, -28.248135469886506)
+    """ 
+    fname = 'point_georef'
+    availkind = ['4invdist', 'nearest']
+    
+    iix = int(ivalue)
+    iiy = int(jvalue)
+    dimx = lonv.shape[1]
+    dimy = lonv.shape[0]
+    iix1 = np.min([dimx,iix+1])
+    iiy1 = np.min([dimy,iiy+1])
+    if kind == '4invdist':
+        dists = np.zeros((4), dtype=np.float)
+        lons = np.zeros((4), dtype=np.float)
+        lats = np.zeros((4), dtype=np.float)
+
+        dists[0] = np.sqrt((iix*1.-ivalue*1.)**2+(iiy*1.-jvalue*1.)**2)
+        dists[1] = np.sqrt((iix1*1.-ivalue*1.)**2+(iiy*1.-jvalue*1.)**2)
+        dists[2] = np.sqrt((iix*1.-ivalue*1.)**2+(iiy1*1.-jvalue*1.)**2)
+        dists[3] = np.sqrt((iix1*1.-ivalue*1.)**2+(iiy1*1.-jvalue*1.)**2)
+
+        lons = [lonv[iiy,iix], lonv[iiy,iix1], lonv[iiy1,iix], lonv[iiy1,iix1]]
+        lats = [latv[iiy,iix], latv[iiy,iix1], latv[iiy1,iix], latv[iiy1,iix1]]
+
+        distsum = np.sum(dists)
+        wghts = dists/distsum
+
+        distsum = np.sum(1./dists)
+        wghts = 1./dists/distsum
+
+        plon = np.sum(lons*wghts)
+        plat = np.sum(lats*wghts)
+
+    elif kind == 'nearest':
+        if (ivalue - int(ivalue)) <= 0.5: iix = int(ivalue)
+        else: iix = iix1
+        if (jvalue - int(jvalue)) <= 0.5: iiy = int(jvalue)
+        else: iiy = iiy1
+
+        plon = lonv[iiy, iix]
+        plat = latv[iiy, iix]
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of georeference method '" + kind +              \
+          "' not ready !!"
+        print '   available ones:', availkind
+        quit(-1)
+
+    return plon, plat
+
 #quit()
 
 
+