1 | ;+ |
---|
2 | ; @file_comments |
---|
3 | ; |
---|
4 | ; |
---|
5 | ; @categories |
---|
6 | ; |
---|
7 | ; |
---|
8 | ; @param FIELD |
---|
9 | ; |
---|
10 | ; |
---|
11 | ; @param RES |
---|
12 | ; |
---|
13 | ; |
---|
14 | ; @param GLAMAXE |
---|
15 | ; |
---|
16 | ; |
---|
17 | ; @param GPHIAXE |
---|
18 | ; |
---|
19 | ; |
---|
20 | ; @keyword ENDPOINTS |
---|
21 | ; |
---|
22 | ; |
---|
23 | ; @keyword BOXZOOM |
---|
24 | ; |
---|
25 | ; |
---|
26 | ; @keyword TYPE |
---|
27 | ; |
---|
28 | ; |
---|
29 | ; @keyword WDEPTH |
---|
30 | ; |
---|
31 | ; |
---|
32 | ; @keyword DIREC |
---|
33 | ; |
---|
34 | ; |
---|
35 | ; @keyword SHOWBUILD |
---|
36 | ; |
---|
37 | ; |
---|
38 | ; @keyword ONLYBOX |
---|
39 | ; |
---|
40 | ; |
---|
41 | ; @keyword _EXTRA |
---|
42 | ; Used to pass your keywords |
---|
43 | ; |
---|
44 | ; @returns |
---|
45 | ; |
---|
46 | ; |
---|
47 | ; @uses |
---|
48 | ; common.pro |
---|
49 | ; |
---|
50 | ; @restrictions |
---|
51 | ; |
---|
52 | ; |
---|
53 | ; @examples |
---|
54 | ; |
---|
55 | ; |
---|
56 | ; @history |
---|
57 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
58 | ; |
---|
59 | ; @version |
---|
60 | ; $Id$ |
---|
61 | ; |
---|
62 | ;- |
---|
63 | |
---|
64 | PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $ |
---|
65 | , BOXZOOM = boxzoom, TYPE = type, WDEPTH = wdepth $ |
---|
66 | , DIREC = direc, SHOWBUILD = showbuild, ONLYBOX = onlybox $ |
---|
67 | , _extra = ex |
---|
68 | ; |
---|
69 | ;--------------------------------------------------------- |
---|
70 | ; include common |
---|
71 | ; |
---|
72 | compile_opt idl2, strictarrsubs |
---|
73 | ; |
---|
74 | @cm_4mesh |
---|
75 | @cm_4data |
---|
76 | @cm_4cal |
---|
77 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
78 | @updatenew |
---|
79 | @updatekwd |
---|
80 | ENDIF |
---|
81 | ;-------------------------------------------------------------- |
---|
82 | ;--------------------- |
---|
83 | ;--------------------- |
---|
84 | ; definition of boxzoom in function of endpoints, then redefinition of the domain |
---|
85 | boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $ |
---|
86 | , min([endpoints[1], endpoints[3]], max = ma13), ma13] |
---|
87 | ; |
---|
88 | minprof = 0. |
---|
89 | profdefault = 200. |
---|
90 | if n_elements(type) EQ 0 then type = 'nothing' |
---|
91 | Case N_Elements(Boxzoom) OF |
---|
92 | 0:localbox = [boxzoom2d, minprof, profdefault] |
---|
93 | 1:localbox = [boxzoom2d, minprof, boxzoom[0]] |
---|
94 | 2:localbox = [boxzoom2d, boxzoom[0]] |
---|
95 | 4:if strpos(type, 'z') NE -1 THEN $ |
---|
96 | localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d |
---|
97 | 5:localbox = [boxzoom2d, minprof, boxzoom[4]] |
---|
98 | 6:localbox = [boxzoom2d, boxzoom[4:5]] |
---|
99 | Else:BEGIN |
---|
100 | print, report('Bad definition of the box') |
---|
101 | stop |
---|
102 | END |
---|
103 | ENDCASE |
---|
104 | nelbox = n_elements(localbox) |
---|
105 | ; |
---|
106 | if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $ |
---|
107 | ELSE grillechoice = vargrid |
---|
108 | domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex |
---|
109 | grille, -1, -1, -1, -1, nx, ny |
---|
110 | ; if less than 10 points where found, we apply domdef over the whole domain |
---|
111 | ; -> problem... why 10 points as a test value??? |
---|
112 | ; how can we find a good test value? |
---|
113 | IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex |
---|
114 | ; We redefine lon1, ... in case findalways has been used in domdef |
---|
115 | lon1 = min([endpoints[0], endpoints[2]], max = lon2) |
---|
116 | lat1 = min([endpoints[1], endpoints[3]], max = lat2) |
---|
117 | ; we extend the box along the z axis -> i that way the plot will be drawn |
---|
118 | ; until its bottom part. |
---|
119 | if strpos(type, 'z') NE -1 THEN BEGIN |
---|
120 | ; We keep yranges (axis z) before changing the boxzoom. |
---|
121 | !y.range = [localbox[nelbox-1], localbox[nelbox-2]] |
---|
122 | if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN |
---|
123 | firstzw = 0 > (firstzw-1) |
---|
124 | lastzw = (lastzw+1) < (jpk-1) |
---|
125 | nzw = lastzw - firstzw + 1 |
---|
126 | ENDIF ELSE BEGIN |
---|
127 | firstzt = 0 > (firstzt-1) |
---|
128 | lastzt = (lastzt+1) < (jpk-1) |
---|
129 | nzt = lastzt - firstzt + 1 |
---|
130 | ENDELSE |
---|
131 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
132 | @updateold |
---|
133 | ENDIF |
---|
134 | ENDIF |
---|
135 | !y.range = [localbox[nelbox-1], localbox[nelbox-2]] |
---|
136 | ; We recuperate the grid on the boxzoom. |
---|
137 | ; |
---|
138 | IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat' |
---|
139 | grille, -1, -1, -1, -1, nx, ny, nz $ |
---|
140 | , firstx, firsty, firstz, lastx, lasty, lastz |
---|
141 | ; the extend the box in the x and y direction in order to maximise |
---|
142 | ; the number of triangles used to build the section |
---|
143 | firstx = 0 > (firstx - 1) |
---|
144 | lastx = (lastx + 1) < (jpi -1) |
---|
145 | firsty = 0 > (firsty - 1) |
---|
146 | lasty = (lasty + 1) < (jpj -1) |
---|
147 | |
---|
148 | domdef, firstx, lastx, firsty, lasty, firstz, lastz $ |
---|
149 | , /index, gridtype = vargrid |
---|
150 | |
---|
151 | IF keyword_set(onlybox) THEN return |
---|
152 | ; |
---|
153 | grille, mask, glam, gphi, gdep, nx, ny, nz $ |
---|
154 | , firstx, firsty, firstz, lastx, lasty, lastz |
---|
155 | |
---|
156 | ;--------------------- |
---|
157 | ; We define the triangulation which will allows us to determinate the section. |
---|
158 | ; We recalculate it because it must be defined on the Earth and on oceans. |
---|
159 | ; Following the direction of the section -rather longitude or rather latitude- |
---|
160 | ; we define the way to triangulate. |
---|
161 | if strpos(type, 'x') NE -1 then BEGIN |
---|
162 | downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] |
---|
163 | tri = definetri(nx, ny, (downward)[*]) |
---|
164 | ENDIF ELSE tri = definetri(nx, ny) |
---|
165 | ; If we have an irregular grid that is periodic, then it is possible that |
---|
166 | ; some of the triangle have a very large size (neighborg points on the |
---|
167 | ; sphere but far away when doing the projection) and should not be |
---|
168 | ; taken into account.. |
---|
169 | IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN |
---|
170 | glamtri = glam[tri] |
---|
171 | glamtri = abs(glamtri - shift(glamtri, 1, 0)) |
---|
172 | good = temporary(glamtri) LT (10.*max(glam)/nx) |
---|
173 | good = where(total(temporary(good), 1) EQ 3) |
---|
174 | tri = (temporary(tri))[*, temporary(good)] |
---|
175 | ENDIF |
---|
176 | ;--------------------- |
---|
177 | ; Equation of the line on which we do the section. |
---|
178 | ;--------------------- |
---|
179 | abc = linearequation(endpoints[0:1], endpoints[2:3]) |
---|
180 | glamtri = glam[tri] |
---|
181 | gphitri = gphi[tri] |
---|
182 | ; Which points of the triangulation are above and below the line? |
---|
183 | if abc[1] NE 0 THEN $ |
---|
184 | test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ |
---|
185 | ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) |
---|
186 | |
---|
187 | zero123 = total(test, 1) |
---|
188 | ; to keep: triangles of the triangulation which are over the line. |
---|
189 | tokeep1 = where(zero123 EQ 1) |
---|
190 | tokeep2 = where(temporary(zero123) EQ 2) |
---|
191 | tokeep = [tokeep1, tokeep2] |
---|
192 | |
---|
193 | test = test[*, tokeep] |
---|
194 | tri = tri[*, tokeep] |
---|
195 | ; Which summit of the triangle is alone in a side of the line? |
---|
196 | single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) |
---|
197 | single1 = single1-(single1/3)*3 |
---|
198 | single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) |
---|
199 | single2 = single2-(single2/3)*3 |
---|
200 | |
---|
201 | undefine, tokeep |
---|
202 | undefine, tokeep1 |
---|
203 | undefine, tokeep2 |
---|
204 | undefine, test |
---|
205 | |
---|
206 | single = [temporary(single1), temporary(single2)] |
---|
207 | ; points1 the point (of the triangle) alone in a side of the line. |
---|
208 | ; point2 the other point of the triangle in the other side of the line. |
---|
209 | point1 = [single, single] |
---|
210 | point2 = [single EQ 0, 1 + (single LE 1)] |
---|
211 | |
---|
212 | undefine, single |
---|
213 | |
---|
214 | ntri = (size(tri))[2] |
---|
215 | index = [lindgen(ntri), lindgen(ntri)] |
---|
216 | |
---|
217 | points1 = tri[point1, index] |
---|
218 | points2 = tri[point2, temporary(index)] |
---|
219 | ; points : complex containing couples of points in a side and the other |
---|
220 | ; side of the line. We have to delete duplicates. |
---|
221 | points = dcomplex(points1, points2) |
---|
222 | points = points[uniq(points, sort(points))] |
---|
223 | symetrique = dcomplex(imaginary(points), double(points)) |
---|
224 | points = points[where(points-shift(temporary(symetrique), 1) NE 0)] |
---|
225 | ; points1 coordinates of the point of the triangle which is alone in a side of the line. |
---|
226 | ; point2 coordinates of the other point of the triangle in the other side of the line. |
---|
227 | points1 = complex(glam[ double(points)], gphi[ double(points)]) |
---|
228 | points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) |
---|
229 | ; droites equations of line whose we look for the intersection wit the section. |
---|
230 | droites = linearequation(points1, points2) |
---|
231 | inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) |
---|
232 | |
---|
233 | ; Geographic coordinates of points we look for on the section. |
---|
234 | glamaxe = float(inter) |
---|
235 | gphiaxe = imaginary(inter) |
---|
236 | ; We arrange them in the growing order between boundaries of the section. |
---|
237 | if strpos(type, 'x') NE -1 then BEGIN |
---|
238 | sort = sort(glamaxe) |
---|
239 | glamaxe = glamaxe[sort] |
---|
240 | inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) |
---|
241 | glamaxe = glamaxe[inbox] |
---|
242 | sort = sort[inbox] |
---|
243 | gphiaxe = gphiaxe[sort] |
---|
244 | ENDIF ELSE BEGIN |
---|
245 | sort = sort(gphiaxe) |
---|
246 | gphiaxe = gphiaxe[sort] |
---|
247 | inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) |
---|
248 | gphiaxe = gphiaxe[inbox] |
---|
249 | sort = sort[inbox] |
---|
250 | glamaxe = glamaxe[sort] |
---|
251 | ENDELSE |
---|
252 | points = points[sort] |
---|
253 | points1 = points1[sort] |
---|
254 | points2 = points2[sort] |
---|
255 | inter = inter[sort] |
---|
256 | poids = abs(points2-inter)/abs(points2-points1) |
---|
257 | |
---|
258 | ;--------------------- |
---|
259 | array = litchamp(field) |
---|
260 | array = fitintobox(array) |
---|
261 | if array[0] EQ -1 THEN BEGIN |
---|
262 | res = -1 |
---|
263 | return |
---|
264 | ENDIF |
---|
265 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
---|
266 | taille = size(array) |
---|
267 | if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN |
---|
268 | direc = 't' |
---|
269 | array = grossemoyenne(array, 't') |
---|
270 | taille = size(array) |
---|
271 | jpt = 1 |
---|
272 | ENDIF |
---|
273 | case 1 of |
---|
274 | ;---------------------------------------------------------------------------- |
---|
275 | ;xy |
---|
276 | ;---------------------------------------------------------------------------- |
---|
277 | taille[0] EQ 2:BEGIN |
---|
278 | value1 = array[double(points)] |
---|
279 | terre = where(value1 GT valmask/10) |
---|
280 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
281 | value2 = array[imaginary(points)] |
---|
282 | terre = where(value2 GT valmask/10) |
---|
283 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
284 | res = poids*value1+(1-poids)*value2 |
---|
285 | END |
---|
286 | ;---------------------------------------------------------------------------- |
---|
287 | ;xyz |
---|
288 | ;---------------------------------------------------------------------------- |
---|
289 | taille[0] EQ 3 AND jpt EQ 1:BEGIN |
---|
290 | npoints = n_elements(points) |
---|
291 | index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) |
---|
292 | value1 = array[index] |
---|
293 | terre = where(value1 GT valmask/10) |
---|
294 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
295 | index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) |
---|
296 | value2 = array[index] |
---|
297 | terre = where(value2 GT valmask/10) |
---|
298 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
299 | poids = poids#replicate(1, nz) |
---|
300 | res = poids*value1+(1-poids)*value2 |
---|
301 | ; average following z ? |
---|
302 | if strpos(type, 'z') EQ -1 then begin |
---|
303 | nan = where(finite(res) EQ 0) |
---|
304 | if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] |
---|
305 | weight = replicate(1, npoints)#e3 |
---|
306 | if nan[0] NE -1 then weight[nan] = !values.f_nan |
---|
307 | totalweight = total(weight, 2, /nan) |
---|
308 | zero = where(totalweight EQ 0) |
---|
309 | if zero[0] NE -1 then totalweight[zero] = !values.f_nan |
---|
310 | res = total(res*weight, 2, /nan)/totalweight |
---|
311 | direc = 'z'+string(byte(testvar(var = toto))) |
---|
312 | endif |
---|
313 | END |
---|
314 | ;---------------------------------------------------------------------------- |
---|
315 | ;xyt |
---|
316 | ;---------------------------------------------------------------------------- |
---|
317 | taille[0] EQ 3 AND jpt NE 1:BEGIN |
---|
318 | npoints = n_elements(points) |
---|
319 | index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) |
---|
320 | value1 = array[index] |
---|
321 | terre = where(value1 GT valmask/10) |
---|
322 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
323 | index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) |
---|
324 | value2 = array[index] |
---|
325 | terre = where(value2 GT valmask/10) |
---|
326 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
327 | poids = poids#replicate(1, jpt) |
---|
328 | res = poids*value1+(1-poids)*value2 |
---|
329 | END |
---|
330 | ;---------------------------------------------------------------------------- |
---|
331 | ;xyzt |
---|
332 | ;---------------------------------------------------------------------------- |
---|
333 | taille[0] EQ 4:BEGIN |
---|
334 | npoints = n_elements(points) |
---|
335 | index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) |
---|
336 | index = reform(index, npoints, nz, jpt, /over) |
---|
337 | value1 = array[index] |
---|
338 | terre = where(value1 GT valmask/10) |
---|
339 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
340 | index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) |
---|
341 | index = reform(index, npoints, nz, jpt, /over) |
---|
342 | value2 = array[index] |
---|
343 | terre = where(value2 GT valmask/10) |
---|
344 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
345 | poids = poids#replicate(1, nz*jpt) |
---|
346 | poids = reform(poids, npoints, nz, jpt, /over) |
---|
347 | res = poids*value1+(1-poids)*value2 |
---|
348 | ; average following z ? |
---|
349 | if strpos(type, 'z') EQ -1 then begin |
---|
350 | nan = where(finite(res) EQ 0) |
---|
351 | if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] |
---|
352 | weight = replicate(1, npoints)#e3 |
---|
353 | weight = weight[*]#replicate(1, jpt) |
---|
354 | weight = reform(weight, npoints, nz, jpt, /over) |
---|
355 | if nan[0] NE -1 then weight[nan] = !values.f_nan |
---|
356 | totalweight = total(weight, 2, /nan) |
---|
357 | zero = where(totalweight EQ 0) |
---|
358 | if zero[0] NE -1 then totalweight[zero] = !values.f_nan |
---|
359 | res = total(res*weight, 2, /nan)/totalweight |
---|
360 | direc = 'z'+string(byte(testvar(var = toto))) |
---|
361 | endif |
---|
362 | END |
---|
363 | endcase |
---|
364 | ;--------------------- |
---|
365 | |
---|
366 | terre = where(finite(res) EQ 0) |
---|
367 | if terre[0] NE -1 then res[terre] = valmask |
---|
368 | |
---|
369 | if n_elements(showbuild) then BEGIN |
---|
370 | winsave = !window |
---|
371 | psave = !p |
---|
372 | xsave = !x |
---|
373 | ysave = !y |
---|
374 | plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '' $ |
---|
375 | , coast_thick = 2, window = showbuild |
---|
376 | !p.title = '' |
---|
377 | !p.subtitle = '' |
---|
378 | |
---|
379 | plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50 |
---|
380 | plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2 |
---|
381 | |
---|
382 | FOR i = 0, n_elements(points1)-1 DO $ |
---|
383 | plots, [float(points1[i]), float(points2[i])] $ |
---|
384 | , [imaginary(points1[i]), imaginary(points2[i])], color = 150 |
---|
385 | |
---|
386 | plots, float(points1), imaginary(points1), color = 150, psym = 1 |
---|
387 | plots, float(points2), imaginary(points2), color = 150, psym = 1 |
---|
388 | plots, float(inter), imaginary(inter), color = 250, psym = 1 |
---|
389 | |
---|
390 | ; ?? bug ?? IF terre[0] NE -1 THEN plots, float(terre[inter]), imaginary(terre[inter]), color = 0, psym = 1 |
---|
391 | |
---|
392 | ; dummy = '' |
---|
393 | ; read, dummy, prompt = 'press return to continue' |
---|
394 | IF !d.name EQ 'PS' THEN erase ELSE wset, winsave |
---|
395 | !p = psave |
---|
396 | !x = xsave |
---|
397 | !y = ysave |
---|
398 | ENDIF |
---|
399 | |
---|
400 | restoreboxparam, 'boxparam4section.dat' |
---|
401 | |
---|
402 | ;--------------------- |
---|
403 | return |
---|
404 | end |
---|