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