Changeset 202 for trunk/SRC/Interpolation

Timestamp:

01/12/07 12:00:06 (17 years ago)

Author:

smasson

Message:

add extrapsmooth + light bugfix in compute_fromirr_bilinear_weigaddr + improve some headers

Location:

trunk/SRC/Interpolation

Files:

• angle.pro (modified) (3 diffs)
• compute_fromirr_bilinear_weigaddr.pro (modified) (7 diffs)
• compute_fromreg_bilinear_weigaddr.pro (modified) (1 diff)
• compute_fromreg_imoms3_weigaddr.pro (modified) (1 diff)
• extrapolate.pro (modified) (2 diffs)
• extrapsmooth.pro (added)
• fromirr.pro (modified) (2 diffs)

trunk/SRC/Interpolation/angle.pro

@@ -8 +8 @@
 ;
 ; @param plam {in}{required}
+; longitude position
 ;
 ; @param pphi {in}{required}
+; latitude position
 ;
 ; @keyword DOUBLE {default=0}
@@ -15 +17 @@
 ;
 ; @returns
-;       gsinu,gcosu : sinus and cosinus of the angle
-;       gsinv,gcosv   between north-south direction
-;       gsint,gcost   and the j-direction of the mesh
-;
-; @restrictions
-; to compute the lateral boundary conditions, we assume that:
-;     (1) the first line is similar to the second line
-;       =>    gcosu[*, 0] = gcosu[*, 1]
-;       =>    gsinu[*, 0] = gsinu[*, 1]
-;     (2) the grid follows OPA x periodicity rule, first column is
-;     equal to the next to last column
-;       =>    gcosv[0, *] = gcosv[jpj-2, *]
-;       =>    gsinv[0, *] = gsinv[jpj-2, *]
-;
-;
-; @history
-;       Original :  96-07 (O. Marti)
-;                   98-06 (G. Madec)
-;       Feb 2005: IDL adaptation S. Masson
-;
-; @version $Id$
+; structure: {x:x, y:y} containing the point position in north stereographic polar projection
+;
+; @hidden
 ;
 ;-
@@ -56 +40 @@
 ;---------
 ;+
-; @file_comments Compute angles between grid lines and direction of the North
+; @file_comments Compute angles between grid lines and direction of the North pole
 ;(fom angle.F,v 2.2 in OPA8.2)
 ;
-; @param fileocemesh {in}{required} a netcdf file that contains (at least):
+; @categories 
+; Interpolation
+;
+; @param fileocemesh {in}{required}{type=scalar string}
+; a netcdf file that contains (at least) the following variables:
 ;        glamu, gphiu: longitudes and latitudes at U-points
 ;        glamv, gphiv: longitudes and latitudes at V-points
 ;        glamf, gphif: longitudes and latitudes at F-points
 ;
-; @param gcosu {in}{required}
-; @param gsinu {in}{required}
-; @param gcosv {in}{required}
-; @param gsinv {in}{required}
-; @param gcost {in}{required}
-; @param gsint {in}{required}
-; @keyword IODIRECTORY the directory path where is located fileocemesh
-; @keyword DOUBLE {default=0} use double precision (default is float)
+; @param gcosu {out}{type=2d array}
+; cosinus of the angle between grid lines at U points and direction of the North pole
+;
+; @param gsinu {out}{type=2d array}
+; sinus of the angle between grid lines at U points and direction of the North pole
+;
+; @param gcosv {out}{type=2d array}
+; cosinus of the angle between grid lines at V points and direction of the North pole
+;
+; @param gsinv {out}{type=2d array}
+; sinus of the angle between grid lines at V points and direction of the North pole
+;
+; @param gcost {out}{type=2d array}
+; cosinus of the angle between grid lines at T points and direction of the North pole
+;
+; @param gsint {out}{type=2d array}
+; sinus of the angle between grid lines at T points and direction of the North pole
+;
+; @keyword IODIRECTORY {type=scalar string}{default=''}
+; the directory path where is located fileocemesh
+;
+; @keyword DOUBLE {type=1 ou 2}{default=0}
+; put 1 to use double precision (default is float)
+;
+; @restrictions
+; to compute the lateral boundary conditions, we assume that:
+;     (1) the first line is similar to the second line
+;       =>    gcosu[*, 0] = gcosu[*, 1]
+;       =>    gsinu[*, 0] = gsinu[*, 1]
+;     (2) the grid follows OPA x periodicity rule, first column is
+;     equal to the next to last column
+;       =>    gcosv[0, *] = gcosv[jpj-2, *]
+;       =>    gsinv[0, *] = gsinv[jpj-2, *]
+;
+; @history
+;       Original :  96-07 (O. Marti)
+;                   98-06 (G. Madec)
+;       Feb 2005: IDL adaptation S. Masson
+;
+; @version $Id$
 ;-
 ;---------

trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

@@ -8 +8 @@
 ; Interpolation
 ;
-; @param olonin {in}{required}
+; @param olonin {in}{required}{type=2d array}
 ; longitude of the input data
 ;
-; @param olat {in}{required}
+; @param olat {in}{required}{type=2d array}
 ; latitude of the input data
 ;
-; @param omsk {in}{required}
+; @param omsk {in}{required}{type=2d array or -1}
 ; land/sea mask of the input data
-;
-; @param alonin {in}{required}
+; put -1 if input data are not masked
+;
+; @param alonin {in}{required}{type=2d array}
 ; longitude of the output data
 ;
-; @param alat {in}{required}
+; @param alat {in}{required}{type=2d array}
 ; latitude of the output data
 ;
-; @param amsk {in}{required}
+; @param amsk {in}{required}{type=2d array or -1}
 ; land/sea mask of the output data
-;
-; @param weig {out}
-; @param addr {out}
+; put -1 if output data are not masked
+;
+; @param weig {out}{type=2d array}
+; (see ADDR)
+;
+; @param addr {out}{type=2d array}
 ; 2D arrays, weig and addr are the weight and addresses used to
 ; perform the interpolation:
@@ -44 +48 @@
 ;  and the weight is redistributed on the remaining "water" corners
 ;  -  points located out of the southern and northern boundaries or in cells
-;  containing only land points are set the the same value as their closest neighbor l
+;  containing only land points are set the the same value as their closest neighbor
 ;
 ; @history
@@ -61 +65 @@
 ;
   jpio = (size(olonin, /dimensions))[0]
-  jpjo = (size(olonin, /dimensions))[1]
+  jpjo = (size(olonin, /dimensions))[1]  
+; mask check
+  IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN omsk = replicate(1b, jpio, jpjo)
+  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja)
+  IF n_elements(omsk) NE jpio*jpjo THEN BEGIN
+    print, 'input grid mask do not have the good size'
+    stop
+  ENDIF
+  IF n_elements(amsk) NE jpia*jpja THEN BEGIN 
+    print, 'output grid mask do not have the good size'
+    stop
+  ENDIF
 ;
 ; longitude, between 0 and 360
@@ -71 +86 @@
   IF out[0] NE -1 THEN olon[out] = olon[out]+360
 ;
-; we work only on the water points
-  owater = where(omsk EQ 1)
-  nowater = n_elements(owater)
+; we work only on the ouput grid water points
   awater = where(amsk EQ 1)
   nawater = n_elements(awater)
@@ -167 +180 @@
       IF ind NE -1 THEN BEGIN
         ind = good[ind]
-; we keep its address
-        foraddr[n] = ind
 ; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1)
 ; in addition we get the corrdinates of the atmospheric point in this new morphed square
@@ -189 +200 @@
         one = where(abs(1-xy) LT 1e-4)
         IF one[0] NE -1 THEN xy[one] = 1
-; some (useless) checks...
-        IF  xy[0] LT 0 OR xy[0] GT 1 THEN stop
-        IF  xy[0] LT 0 OR xy[0] GT 1 THEN stop
+; some checks...
+        tmpmsk = omsk[celladdr[*, ind]]
+        CASE 1 OF
+          xy[0] LT 0 OR xy[0] GT 1 : stop
+          xy[1] LT 0 OR xy[1] GT 1 : stop
+          xy[0] EQ 0 AND xy[1] EQ 0 AND tmpmsk[0] EQ 0 : foraddr[n] = -1
+          xy[0] EQ 1 AND xy[1] EQ 0 AND tmpmsk[1] EQ 0 : foraddr[n] = -1
+          xy[0] EQ 1 AND xy[1] EQ 1 AND tmpmsk[2] EQ 0 : foraddr[n] = -1
+          xy[0] EQ 0 AND xy[1] EQ 1 AND tmpmsk[3] EQ 0 : foraddr[n] = -1
+          xy[0] EQ 0 AND (tmpmsk[0]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
+          xy[0] EQ 1 AND (tmpmsk[1]+tmpmsk[2]) EQ 0    : foraddr[n] = -1
+          xy[1] EQ 0 AND (tmpmsk[0]+tmpmsk[1]) EQ 0    : foraddr[n] = -1
+          xy[1] EQ 1 AND (tmpmsk[2]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
+          ELSE: BEGIN 
+; we keep its address
+        foraddr[n] = ind
 ; keep the new coordinates
-        forweight[n, 0] = xy[0]
-        forweight[n, 1] = xy[1]
+            forweight[n, 0] = xy[0]
+            forweight[n, 1] = xy[1]
+          END
+        ENDCASE
+
+
+
       ENDIF ELSE foraddr[n] = -1
     ENDIF ELSE foraddr[n] = -1
@@ -237 +266 @@
 ; weights normalization
   totalweig = total(newaweig, 1)
+  IF min(totalweig, max = ma) EQ 0 then stop
+  IF ma GT 1 then stop
   newaweig = newaweig/(replicate(1., 4)#totalweig)
   totalweig = total(newaweig, 1)

trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro

@@ -7 +7 @@
 ; Interpolation
 ;
-; @param alonin{in}{required}
+; @param alonin{in}{required}{type=2d array}
 ; longitude of the input data
 ;
-; @param alatin {in}{required}
+; @param alatin {in}{required}{type=2d array}
 ; latitude of the input data
 ;
-; @param olonin {in}{required}
+; @param olonin {in}{required}{type=2d array}
 ; longitude of the output data
 ;
-; @param olat {in}{required}
+; @param olat {in}{required}{type=2d array}
 ; latitude of the output data
 ;
-; @keyword NONORTHERNLINE
-; activate if you don't want to take into
+; @keyword NONORTHERNLINE {type=scalar 0 or 1}{default=0}
+; put 1 if you don't want to take into
 ; account the northen line of the input data when performing the interpolation.
 ;
-; @keyword NOSOUTHERNLINE
-; activate if you don't want to take into
+; @keyword NOSOUTHERNLINE {type=scalar 0 or 1}{default=0}
+; put 1 if you don't want to take into
 ; account the southern line of the input data when performing the interpolation.
 ;
-; @param weig {out}
-; @param addr {out}
+; @param weig {out}{type=2d array}
+; (see ADDR)
+;
+; @param addr {out}{type=2d array}
 ; 2D arrays, weig and addr are the weight and addresses used to
 ; perform the interpolation:

trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro

@@ -8 +8 @@
 ; Interpolation
 ;
-; @param alonin {in}{required}
+; @param alonin {in}{required}{type=2d array}
 ; longitude of the input data
 ;
-; @param alatin {in}{required}
+; @param alatin {in}{required}{type=2d array}
 ; latitude of the input data
 ;
-; @param olonin {in}{required}
+; @param olonin {in}{required}{type=2d array}
 ; longitude of the output data
-; @param olat {in}{required}
+;
+; @param olat {in}{required}{type=2d array}
 ; latitude of the output data
 ;
-; @keyword NONORTHERNLINE
-; @keyword NOSOUTHERNLINE
-; activate if you don't want to take into account the northen/southern line
-; of the input data when performing the interpolation.
-;
-; @param weig {out}
-; @param addr {out}
+; @keyword NONORTHERNLINE {type=scalar 0 or 1}{default=0}
+; put 1 if you don't want to take into
+; account the northen line of the input data when performing the interpolation.
+;
+; @keyword NOSOUTHERNLINE {type=scalar 0 or 1}{default=0}
+; put 1 if you don't want to take into
+; account the southern line of the input data when performing the interpolation.
+;
+; @param weig {out}{type=2d array}
+; (see ADDR)
+;
+; @param addr {out}{type=2d array}
 ; 2D arrays, weig and addr are the weight and addresses used to
 ; perform the interpolation:

trunk/SRC/Interpolation/extrapolate.pro

@@ -1 +1 @@
 ;+
 ; @file_comments
-; extrapolate data (zinput) where maskinput eq 0 by filling
-; step by step the coastline points with the mean value of the 8 neighbourgs.
+; extrapolate data (zinput) where maskinput eq 0 by filling step by 
+; step the coastline points with the mean value of the 8 neighbourgs
+; (weighted by their mask value).
 ;
 ; @categories 
 ; Interpolation
 ;
-; @param zinput {in}{required}
+; @param zinput {in}{required}{type=2d array}
 ; data to be extrapolate
 ;
-; @param maskinput {in}{required}
+; @param maskinput {in}{required}{type=2d array or -1}
+; a 2D array, the land-sea mask of the output data (1 on ocean, 0 on land)
+; put -1 if input data are not masked
 ;
-; @param nb_iteration {in}{optional}
-; number of iteration
+; @param nb_iteration {in}{optional}{type=integer scalar}{default=10.E20}
+; Maximum number if iterations done in the extrapolation process. If there
+; is no more masked values we exit extrapolate before reaching nb_iteration
+; (to be sure to fill everything, you can use a very large value)
 ;
-; @keyword x_periodic
-; @keyword MINVAL
-; @keyword MAXVAL
-; @keyword GE0
+; @keyword x_periodic {type=scalar, 0 or 1}{default=0}
+; put 1 to specify that the data are periodic along x axis
+;
+; @keyword MINVAL {type=scalar}{default=not used}
+; to specify a minimum value to the extrapolated values
+;
+; @keyword MAXVAL {type=scalar}{default=not used}
+; to specify a maximum value to the extrapolated values
+;
+; @keyword GE0 {type=scalar 0 or 1}{default=0}
+; put 1 to force the extrapolated values to be larger than 0, same as using minval=0.
+;
+; @returns {type=2d array}
+; the extrapolated array
+;
+; @examples
+; IDL> a=extrapolate(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
+; IDL> tvplus, a
+; IDL> tvplus, a*(1-tmask[*,*,0])
+; get the coastline:
+; IDL> a=extrapolate(tmask[*,*,0],tmask[*,*,0],1,/x_periodic)
+; IDL> tvplus, a-tmask[*,*,0]
+;
+; @history
+;  Originaly written by G. Roulet
+;  Sebastien Masson (smasson\@lodyc.jussieu.fr)
 ;
 ; @version $Id$
 ;
 ;-
-FUNCTION extrapolate, zinput, maskinput, nb_iteration, x_periodic = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0
+FUNCTION extrapolate, zinput, maskinput, nb_iteration, x_periodic = x_periodic $
+                      , MINVAL = minval, MAXVAL = maxval, GE0 = ge0
 ;
   compile_opt idl2, strictarrsubs
@@ -42 +70 @@
 ;
   ztmp = bytarr(nx+2, ny+2)
+  IF n_elements(maskinput) EQ 1 AND maskinput[0] EQ -1 THEN maskinput = replicate(1b, nx, ny)
+  IF n_elements(maskinput) NE nx*ny THEN BEGIN 
+    print, 'input grid mask do not have the good size'
+    return, -1
+  ENDIF
   ztmp[1:nx, 1:ny] = byte(maskinput)
   msk = temporary(ztmp)

trunk/SRC/Interpolation/fromirr.pro

@@ -20 +20 @@
 ; a 2D array defining the latitude of the input data.
 ;
-; @param mskin {in}{optional}{type=2d array}
+; @param mskin {in}{optional}{type=2d array or -1}
 ; a 2D array, the land-sea mask of the input data (1 on ocean, 0 on land)
+; put -1 if input data are not masked
 ;
 ; @param lonout {in}{optional}{type=1d or 2d array}
@@ -29 +30 @@
 ; 1D or 2D array defining the latitude of the output data.
 ;
-; @param mskout {in}{required}{type=2d array}
+; @param mskout {in}{required}{type=2d array or -1}
 ; a 2D array, the land-sea mask of the output data (1 on ocean, 0 on land)
+; put -1 if output data are not masked
 ;
 ; @keyword WEIG {type=2d array}