source: trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro @ 257

Last change on this file since 257 was 257, checked in by smasson, 17 years ago

allow interpolation from ORCA2 whithout mask

  • Property svn:keywords set to Id
File size: 11.0 KB
Line 
1;+
2;
3; @file_comments
4; compute the weight and address needed to interpolate data from
5; an "irregular 2D grid" (defined as a grid made of quadrilateral cells)
6; to any grid using the bilinear method
7;
8; @categories
9; Interpolation
10;
11; @param olonin {in}{required}{type=2d array}
12; longitude of the input data
13;
14; @param olat {in}{required}{type=2d array}
15; latitude of the input data
16;
17; @param omsk {in}{required}{type=2d array or -1}
18; land/sea mask of the input data
19; put -1 if input data are not masked
20;
21; @param alonin {in}{required}{type=2d array}
22; longitude of the output data
23;
24; @param alat {in}{required}{type=2d array}
25; latitude of the output data
26;
27; @param amsk {in}{required}{type=2d array or -1}
28; land/sea mask of the output data
29; put -1 if output data are not masked
30;
31; @param weig {out}{type=2d array}
32; (see ADDR)
33;
34; @param addr {out}{type=2d array}
35; 2D arrays, weig and addr are the weight and addresses used to
36; perform the interpolation:
37;  dataout = total(weig*datain[addr], 1)
38;  dataout = reform(dataout, jpia, jpja, /over)
39;
40; @restrictions
41;  -  the input grid must be an "irregular 2D grid", defined as a grid made
42;  of quadrilateral cells which corners positions are defined with olonin and olat
43;  -  We supposed the data are located on a sphere, with a periodicity along
44;  the longitude
45;  -  to perform the bilinear interpolation within quadrilateral cells, we
46;  first morph the cell into a square cell and then compute the bilinear
47;  interpolation.
48;  -  if some corners of the cell are land points, their weight is set to 0
49;  and the weight is redistributed on the remaining "water" corners
50;  -  points located out of the southern and northern boundaries or in cells
51;  containing only land points are set the same value as their closest neighbor
52;
53; @history
54;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr)
55;
56; @version
57; $Id$
58;
59;-
60;
61PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk, weig, addr
62;
63  compile_opt idl2, strictarrsubs
64;
65  jpia = (size(alonin, /dimensions))[0]
66  jpja = (size(alonin, /dimensions))[1]
67;
68  jpio = (size(olonin, /dimensions))[0]
69  jpjo = (size(olonin, /dimensions))[1]
70; mask check
71  IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN
72    omsk = replicate(1b, jpio, jpjo)
73   IF jpio EQ 182 AND jpjo EQ 149 THEN BEGIN
74      lontmp = (olonin[1:180, *] + 3600) MOD 360
75      lattmp = olat[1:180, *]
76      bad = abs(abs(lontmp - shift(lontmp, 1, 0)) - 180) LT 176 $
77            OR  abs(abs(lontmp - shift(lontmp, 0, 1)) - 180) LT 176 $
78            OR  abs(abs(lontmp - shift(lontmp, 1, 1)) - 180) LT 176 $
79            OR  abs(abs(lattmp - shift(lattmp, 1, 0)) - 180) LT 176 $
80            OR  abs(abs(lattmp - shift(lattmp, 0, 1)) - 180) LT 176 $
81            OR  abs(abs(lattmp - shift(lattmp, 1, 1)) - 180) LT 176
82      bad = bad AND lattmp LT 50
83      omsk[1:180, 1:148] = 1b - bad[*, 1:148]
84      omsk[0, *] = omsk[180, *]
85      omsk[181, *] = omsk[1, *]
86    ENDIF ELSE omsk = replicate(1b, jpio, jpjo)
87  ENDIF
88  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja)
89  IF n_elements(omsk) NE jpio*jpjo THEN BEGIN
90    ras = report('input grid mask do not have the good size')
91    stop
92  ENDIF
93  IF n_elements(amsk) NE jpia*jpja THEN BEGIN
94    ras = report('output grid mask do not have the good size')
95    stop
96  ENDIF
97;
98; longitude, between 0 and 360
99  alon = alonin MOD 360
100  out = where(alon LT 0)
101  IF out[0] NE -1 THEN alon[out] = alon[out]+360
102  olon = olonin MOD 360
103  out = where(olon LT 0)
104  IF out[0] NE -1 THEN olon[out] = olon[out]+360
105;
106; we work only on the ouput grid water points
107  awater = where(amsk EQ 1)
108  nawater = n_elements(awater)
109;
110; define all cells that have corners located at olon, olat
111; we define the cell with the address of all corners
112;
113;            3        2
114;             +------+
115;             |      |
116;             |      |
117;             |      |
118;             +------+
119;            0        1
120;
121  alladdr = lindgen(jpio, jpjo-1)
122  celladdr = lonarr(4, jpio*(jpjo-1))
123  celladdr[0, *] = alladdr
124  celladdr[1, *] = shift(alladdr, -1)
125  celladdr[2, *] = shift(alladdr + jpio, -1)
126  celladdr[3, *] = alladdr + jpio
127;
128; list the cells that have at least 1 ocean point as corner
129  good = where(total(omsk[celladdr], 1) GT 0)
130; keep only those cells
131  celladdr = celladdr[*, temporary(good)]
132;
133  xcell = olon[celladdr]
134  minxcell = min(xcell, dimension = 1, max = maxxcell)
135  ycell = olat[celladdr]
136  minycell = min(ycell, dimension = 1, max = maxycell)
137; foraddr: address of the ocean water cell associated to each atmosphere water point
138  foraddr = lonarr(nawater)
139; forweight: x/y position of the atmosphere water point in the ocean water cell
140  forweight = dblarr(nawater, 2)
141;
142; Loop on all the water point of the atmosphere
143; We look for which ocean water cell contains the atmosphere water point
144;
145  delta = max([(360./jpio), (180./jpjo)])* 4.
146  FOR n = 0L, nawater-1 DO BEGIN
147; control print
148    IF (n MOD 5000) EQ 0 THEN print, n
149; longitude and latitude of the atmosphere water point
150    xx = alon[awater[n]]
151    yy = alat[awater[n]]
152; 1) we reduce the number of ocean cells for which we will try to know if
153; xx,yy is inside.
154    CASE 1 OF
155; if we are near the norh pole
156      yy GE (90-delta):BEGIN
157        lat1 = 90-2*delta
158        good = where(maxycell GE lat1)
159        onsphere = 1
160      END
161; if we are near the longitude periodicity area
162      xx LE delta OR xx GE (360-delta):BEGIN
163        lat1 = yy-delta
164        lat2 = yy+delta
165        good = where((minxcell LE 2*delta OR maxxcell GE (360-2*delta)) AND maxycell GE lat1 AND minycell LE lat2)
166        onsphere = 1
167      END
168; other cases
169      ELSE:BEGIN
170        lon1 = xx-delta
171        lon2 = xx+delta
172        lat1 = yy-delta
173        lat2 = yy+delta
174        good = where(maxxcell GE lon1 AND minxcell LE lon2 AND maxycell GE lat1 AND minycell le lat2)
175; ORCA cases : orca grid is irregular only northward of 40N
176        CASE 1 OF
177          jpio EQ 92   AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40
178          jpio EQ 180  AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40
179          jpio EQ 720  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40
180          jpio EQ 1440 AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40
181          ELSE:onsphere = 1
182        ENDCASE
183      END
184    ENDCASE
185; we found a short list of possible ocean water cells containing the atmosphere water point
186    IF good[0] NE -1 THEN BEGIN
187; in which cell is located the atmosphere water point?
188; Warning! inquad use clockwise quadrilateral definition
189      ind = inquad(xx, yy $
190                   , xcell[0, good], ycell[0, good] $
191                   , xcell[3, good], ycell[3, good] $
192                   , xcell[2, good], ycell[2, good] $
193                   , xcell[1, good], ycell[1, good] $
194                   , onsphere = onsphere, newcoord = newcoord, /noprint)
195; keep only the first cell (if the atmospheric point was located in several ocean cells)
196      ind = ind[0]
197; we found one ocean water cell containing the atmosphere water point
198      IF ind NE -1 THEN BEGIN
199        ind = good[ind]
200; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1)
201; in addition we get the coordinates of the atmospheric point in this new morphed square
202        IF onsphere THEN BEGIN
203; Warning! quadrilateral2square use anticlockwise quadrilateral definition
204          xy = quadrilateral2square(newcoord[0, 0], newcoord[1, 0] $
205                                    , newcoord[0, 3], newcoord[1, 3] $
206                                    , newcoord[0, 2], newcoord[1, 2] $
207                                    , newcoord[0, 1], newcoord[1, 1] $
208                                    , newcoord[0, 4], newcoord[1, 4])
209        ENDIF ELSE BEGIN
210          xy = quadrilateral2square(xcell[0, ind], ycell[0, ind] $
211                                    , xcell[1, ind], ycell[1, ind] $
212                                    , xcell[2, ind], ycell[2, ind] $
213                                    , xcell[3, ind], ycell[3, ind], xx, yy)
214        ENDELSE
215; take care of rounding errors...
216        zero = where(abs(xy) LT 1e-4)
217        IF zero[0] NE -1 THEN xy[zero] = 0
218        one = where(abs(1-xy) LT 1e-4)
219        IF one[0] NE -1 THEN xy[one] = 1
220; some checks...
221        tmpmsk = omsk[celladdr[*, ind]]
222        CASE 1 OF
223          xy[0] LT 0 OR xy[0] GT 1 : stop
224          xy[1] LT 0 OR xy[1] GT 1 : stop
225          xy[0] EQ 0 AND xy[1] EQ 0 AND tmpmsk[0] EQ 0 : foraddr[n] = -1
226          xy[0] EQ 1 AND xy[1] EQ 0 AND tmpmsk[1] EQ 0 : foraddr[n] = -1
227          xy[0] EQ 1 AND xy[1] EQ 1 AND tmpmsk[2] EQ 0 : foraddr[n] = -1
228          xy[0] EQ 0 AND xy[1] EQ 1 AND tmpmsk[3] EQ 0 : foraddr[n] = -1
229          xy[0] EQ 0 AND (tmpmsk[0]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
230          xy[0] EQ 1 AND (tmpmsk[1]+tmpmsk[2]) EQ 0    : foraddr[n] = -1
231          xy[1] EQ 0 AND (tmpmsk[0]+tmpmsk[1]) EQ 0    : foraddr[n] = -1
232          xy[1] EQ 1 AND (tmpmsk[2]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
233          ELSE: BEGIN
234; we keep its address
235        foraddr[n] = ind
236; keep the new coordinates
237            forweight[n, 0] = xy[0]
238            forweight[n, 1] = xy[1]
239          END
240        ENDCASE
241
242
243
244      ENDIF ELSE foraddr[n] = -1
245    ENDIF ELSE foraddr[n] = -1
246  ENDFOR
247; do we have some water atmospheric points that are not located in an water oceanic cell?
248  bad = where(foraddr EQ -1)
249  IF bad[0] NE -1 THEN BEGIN
250; yes?
251; we look for neighbor water atmospheric point located in water oceanic cell
252    badaddr = awater[bad]
253    good = where(foraddr NE -1)
254; list the atmospheric points located in water oceanic cell
255    goodaddr = awater[good]
256; there longitude and latitude
257    goodlon = alon[goodaddr]
258    goodlat = alat[goodaddr]
259; for all the bad points, look for a neighbor
260    neig = lonarr(n_elements(bad))
261    FOR i = 0, n_elements(bad)-1 DO BEGIN
262      neig[i] = (neighbor(alon[badaddr[i]], alat[badaddr[i]], goodlon, goodlat, /sphere))[0]
263    ENDFOR
264; get the address regarding foraddr
265    neig = good[neig]
266; associate each bad point with its neighbor (get its address and weight)
267    foraddr[bad] = foraddr[neig]
268    forweight[bad, *] = forweight[neig, *]
269  endif
270; transform the address of the ocean cell into the address of its 4 corners
271  newaaddr = celladdr[*, temporary(foraddr)]
272; now we compute the weight to give at each corner
273  newaweig = dblarr(4, nawater)
274  a = reform(forweight[*, 0], 1, nawater)
275  b = reform(forweight[*, 1], 1, nawater)
276  forweight =  -1 ; free memory
277  newaweig = [(1-a)*(1-b), (1-b)*a, a*b, (1-a)*b]
278  a = -1 &  b = -1 ; free memory
279; mask the weight to suppress the corner located on land
280  newaweig = newaweig*((omsk)[newaaddr])
281  totalweig = total(newaweig, 1)
282; for cell with some land corner,
283; we have to redistribute the weight on the reaining water corners
284; weights normalization
285  totalweig = total(newaweig, 1)
286  IF min(totalweig, max = ma) EQ 0 then stop
287  IF ma GT 1 then stop
288  newaweig = newaweig/(replicate(1., 4)#totalweig)
289  totalweig = total(newaweig, 1)
290
291; weights
292  weig = dblarr(4, jpia*jpja)
293  weig[*, awater] = temporary(newaweig)
294; address
295  addr = dblarr(4, jpia*jpja)
296  addr[*, awater] = temporary(newaaddr)
297;
298  RETURN
299END
Note: See TracBrowser for help on using the repository browser.