Context Navigation

section.pro @ 134

Last change on this file since 134 was 134, checked in by navarro, 18 years ago
change *.pro file properties (del eof-style, del executable, set keywords Id
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 13.7 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; NAME:
6	;
7	; PURPOSE:
8	;
9	; CATEGORY:
10	;
11	; CALLING SEQUENCE:
12	;
13	; INPUTS:
14	;
15	; KEYWORD PARAMETERS:
16	;
17	; OUTPUTS:
18	;
19	; COMMON BLOCKS:common.pro
20	;
21	; SIDE EFFECTS:
22	;
23	; RESTRICTIONS:
24	;
25	; EXAMPLE:
26	;
27	; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
28	;
29	;-
30	;------------------------------------------------------------
31	;------------------------------------------------------------
32	;------------------------------------------------------------
33
34	PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $
35	, BOXZOOM = boxzoom, TYPE = type, WDEPTH = wdepth $
36	, DIREC = direc, SHOWBUILD = showbuild, ONLYBOX = onlybox $
37	, _extra = ex
38	;
39	;---------------------------------------------------------
40	; include common
41	;
42	compile_opt idl2, strictarrsubs
43	;
44	@cm_4mesh
45	@cm_4data
46	@cm_4cal
47	IF NOT keyword_set(key_forgetold) THEN BEGIN
48	@updatenew
49	@updatekwd
50	ENDIF
51	;--------------------------------------------------------------
52	;---------------------
53	;---------------------
54	; definition de boxzoom en fonction de endpoints
55	; puis redefinition du domaine
56	boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $
57	, min([endpoints[1], endpoints[3]], max = ma13), ma13]
58	;
59	minprof = 0.
60	profdefault = 200.
61	if n_elements(type) EQ 0 then type = 'nothing'
62	Case N_Elements(Boxzoom) OF
63	0:localbox = [boxzoom2d, minprof, profdefault]
64	1:localbox = [boxzoom2d, minprof, boxzoom[0]]
65	2:localbox = [boxzoom2d, boxzoom[0]]
66	4:if strpos(type, 'z') NE -1 THEN $
67	localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d
68	5:localbox = [boxzoom2d, minprof, boxzoom[4]]
69	6:localbox = [boxzoom2d, boxzoom[4:5]]
70	Else:BEGIN
71	print, report('Bad definition of the box')
72	stop
73	END
74	ENDCASE
75	nelbox = n_elements(localbox)
76	;
77	if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $
78	ELSE grillechoice = vargrid
79	domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex
80	grille, -1, -1, -1, -1, nx, ny
81	; if less than 10 points where found, we apply domdef over the whole domain
82	; -> problem... why 10 points as a test value???
83	; how can we find a good test value?
84	IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex
85	; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef
86	lon1 = min([endpoints[0], endpoints[2]], max = lon2)
87	lat1 = min([endpoints[1], endpoints[3]], max = lat2)
88	; we extend the box along the z axis -> i that way the plot will be drawn
89	; until its bottom part.
90	if strpos(type, 'z') NE -1 THEN BEGIN
91	;on garde les yranges (axe z) avant de changer la boxzoom.
92	!y.range = [localbox[nelbox-1], localbox[nelbox-2]]
93	if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN
94	firstzw = 0 > (firstzw-1)
95	lastzw = (lastzw+1) < (jpk-1)
96	nzw = lastzw - firstzw + 1
97	ENDIF ELSE BEGIN
98	firstzt = 0 > (firstzt-1)
99	lastzt = (lastzt+1) < (jpk-1)
100	nzt = lastzt - firstzt + 1
101	ENDELSE
102	IF NOT keyword_set(key_forgetold) THEN BEGIN
103	@updateold
104	ENDIF
105	ENDIF
106	!y.range = [localbox[nelbox-1], localbox[nelbox-2]]
107	; on recupere la grille sur la boxzoom
108	;
109	IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat'
110	grille, -1, -1, -1, -1, nx, ny, nz $
111	, firstx, firsty, firstz, lastx, lasty, lastz
112	; the extend the box in the x and y direction in order to maximise
113	; the number of triangles used to build the section
114	firstx = 0 > (firstx - 1)
115	lastx = (lastx + 1) < (jpi -1)
116	firsty = 0 > (firsty - 1)
117	lasty = (lasty + 1) < (jpj -1)
118
119	domdef, firstx, lastx, firsty, lasty, firstz, lastz $
120	, /index, gridtype = vargrid
121
122	IF keyword_set(onlybox) THEN return
123	;
124	grille, mask, glam, gphi, gdep, nx, ny, nz $
125	, firstx, firsty, firstz, lastx, lasty, lastz
126
127	;---------------------
128	; on definit la triangulation qui va nous permetre de determiner la
129	; section. on la recalcule car elle doit etre definie sur la terre
130	; aussi bien que sur la mer. suivant le sens de la section -plutot
131	; longitude ou plutot latitude- on definit la facon de trianguler
132	if strpos(type, 'x') NE -1 then BEGIN
133	downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2]
134	tri = definetri(nx, ny, (downward)[*])
135	ENDIF ELSE tri = definetri(nx, ny)
136	; If we have an irregular grid that is periodic, then it is possible that
137	; some of the triangle have a very large size (neighborg points on the
138	; sphere but far away when doing the projection) and should not be
139	; taken into account..
140	IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN
141	glamtri = glam[tri]
142	glamtri = abs(glamtri - shift(glamtri, 1, 0))
143	good = temporary(glamtri) LT (10.*max(glam)/nx)
144	good = where(total(temporary(good), 1) EQ 3)
145	tri = (temporary(tri))[*, temporary(good)]
146	ENDIF
147	;---------------------
148	; equation de la droite suivant laquelle on fait la section
149	;---------------------
150	abc = linearequation(endpoints[0:1], endpoints[2:3])
151	glamtri = glam[tri]
152	gphitri = gphi[tri]
153	; quels sont les points de la triangulation qui sont au dessus et au
154	; dessous de la droite ?
155	if abc[1] NE 0 THEN $
156	test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $
157	ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0])
158
159	zero123 = total(test, 1)
160	; to keep: triangles de la triangulation qui sont a cheval sur la droite.
161	tokeep1 = where(zero123 EQ 1)
162	tokeep2 = where(temporary(zero123) EQ 2)
163	tokeep = [tokeep1, tokeep2]
164
165	test = test[*, tokeep]
166	tri = tri[*, tokeep]
167	; quel est le sommet du triangle qui est seul d''un cote de la droite?
168	single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1)
169	single1 = single1-(single1/3)*3
170	single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0)
171	single2 = single2-(single2/3)*3
172
173	undefine, tokeep
174	undefine, tokeep1
175	undefine, tokeep2
176	undefine, test
177
178	single = [temporary(single1), temporary(single2)]
179	; points1 le point du triangles qui est seul d''un cote de la droite.
180	; point2 l''autre point du triangle de l''autre cote de la droite
181	point1 = [single, single]
182	point2 = [single EQ 0, 1 + (single LE 1)]
183
184	undefine, single
185
186	ntri = (size(tri))[2]
187	index = [lindgen(ntri), lindgen(ntri)]
188
189	points1 = tri[point1, index]
190	points2 = tri[point2, temporary(index)]
191	; points : complexe contenant les couples de points de part et
192	; d''autre de la droite. Ils faut supprimer les doublons.
193	points = dcomplex(points1, points2)
194	points = points[uniq(points, sort(points))]
195	symetrique = dcomplex(imaginary(points), double(points))
196	points = points[where(points-shift(temporary(symetrique), 1) NE 0)]
197	; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite.
198	; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite
199	points1 = complex(glam[ double(points)], gphi[ double(points)])
200	points2 = complex(glam[imaginary(points)], gphi[imaginary(points)])
201	; droites les equations des droites dont on cherche l''intersection
202	; avec la section.
203	droites = linearequation(points1, points2)
204	inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) )
205
206	; les ccordonnes geographiques des points que l''on cherche sur la section.
207	glamaxe = float(inter)
208	gphiaxe = imaginary(inter)
209	; on les range ds l''ordre croissant entre les bornes de la section
210	if strpos(type, 'x') NE -1 then BEGIN
211	sort = sort(glamaxe)
212	glamaxe = glamaxe[sort]
213	inbox = where(glamaxe GE lon1 AND glamaxe LE lon2)
214	glamaxe = glamaxe[inbox]
215	sort = sort[inbox]
216	gphiaxe = gphiaxe[sort]
217	ENDIF ELSE BEGIN
218	sort = sort(gphiaxe)
219	gphiaxe = gphiaxe[sort]
220	inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2)
221	gphiaxe = gphiaxe[inbox]
222	sort = sort[inbox]
223	glamaxe = glamaxe[sort]
224	ENDELSE
225	points = points[sort]
226	points1 = points1[sort]
227	points2 = points2[sort]
228	inter = inter[sort]
229	poids = abs(points2-inter)/abs(points2-points1)
230
231	;---------------------
232	array = litchamp(field)
233	array = fitintobox(array)
234	if array[0] EQ -1 THEN BEGIN
235	res = -1
236	return
237	ENDIF
238	if n_elements(valmask) EQ 0 THEN valmask = 1e20
239	taille = size(array)
240	if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN
241	direc = 't'
242	array = grossemoyenne(array, 't')
243	taille = size(array)
244	jpt = 1
245	ENDIF
246	case 1 of
247	;----------------------------------------------------------------------------
248	;xy
249	;----------------------------------------------------------------------------
250	taille[0] EQ 2:BEGIN
251	value1 = array[double(points)]
252	terre = where(value1 GT valmask/10)
253	if terre[0] NE -1 then value1[terre] = !values.f_nan
254	value2 = array[imaginary(points)]
255	terre = where(value2 GT valmask/10)
256	if terre[0] NE -1 then value2[terre] = !values.f_nan
257	res = poidsvalue1+(1-poids)value2
258	END
259	;----------------------------------------------------------------------------
260	;xyz
261	;----------------------------------------------------------------------------
262	taille[0] EQ 3 AND jpt EQ 1:BEGIN
263	npoints = n_elements(points)
264	index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
265	value1 = array[index]
266	terre = where(value1 GT valmask/10)
267	if terre[0] NE -1 then value1[terre] = !values.f_nan
268	index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
269	value2 = array[index]
270	terre = where(value2 GT valmask/10)
271	if terre[0] NE -1 then value2[terre] = !values.f_nan
272	poids = poids#replicate(1, nz)
273	res = poidsvalue1+(1-poids)value2
274	; moyenne suivant z ?
275	if strpos(type, 'z') EQ -1 then begin
276	nan = where(finite(res) EQ 0)
277	if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
278	weight = replicate(1, npoints)#e3
279	if nan[0] NE -1 then weight[nan] = !values.f_nan
280	totalweight = total(weight, 2, /nan)
281	zero = where(totalweight EQ 0)
282	if zero[0] NE -1 then totalweight[zero] = !values.f_nan
283	res = total(res*weight, 2, /nan)/totalweight
284	direc = 'z'+string(byte(testvar(var = toto)))
285	endif
286	END
287	;----------------------------------------------------------------------------
288	;xyt
289	;----------------------------------------------------------------------------
290	taille[0] EQ 3 AND jpt NE 1:BEGIN
291	npoints = n_elements(points)
292	index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
293	value1 = array[index]
294	terre = where(value1 GT valmask/10)
295	if terre[0] NE -1 then value1[terre] = !values.f_nan
296	index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
297	value2 = array[index]
298	terre = where(value2 GT valmask/10)
299	if terre[0] NE -1 then value2[terre] = !values.f_nan
300	poids = poids#replicate(1, jpt)
301	res = poidsvalue1+(1-poids)value2
302	END
303	;----------------------------------------------------------------------------
304	;xyzt
305	;----------------------------------------------------------------------------
306	taille[0] EQ 4:BEGIN
307	npoints = n_elements(points)
308	index = double(points)#replicate(1, nzjpt)+replicate(nxny, npoints)#lindgen(nz*jpt)
309	index = reform(index, npoints, nz, jpt, /over)
310	value1 = array[index]
311	terre = where(value1 GT valmask/10)
312	if terre[0] NE -1 then value1[terre] = !values.f_nan
313	index = imaginary(points)#replicate(1, nzjpt)+replicate(nxny, npoints)#lindgen(nz*jpt)
314	index = reform(index, npoints, nz, jpt, /over)
315	value2 = array[index]
316	terre = where(value2 GT valmask/10)
317	if terre[0] NE -1 then value2[terre] = !values.f_nan
318	poids = poids#replicate(1, nz*jpt)
319	poids = reform(poids, npoints, nz, jpt, /over)
320	res = poidsvalue1+(1-poids)value2
321	; moyenne suivant z ?
322	if strpos(type, 'z') EQ -1 then begin
323	nan = where(finite(res) EQ 0)
324	if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
325	weight = replicate(1, npoints)#e3
326	weight = weight[*]#replicate(1, jpt)
327	weight = reform(weight, npoints, nz, jpt, /over)
328	if nan[0] NE -1 then weight[nan] = !values.f_nan
329	totalweight = total(weight, 2, /nan)
330	zero = where(totalweight EQ 0)
331	if zero[0] NE -1 then totalweight[zero] = !values.f_nan
332	res = total(res*weight, 2, /nan)/totalweight
333	direc = 'z'+string(byte(testvar(var = toto)))
334	endif
335	END
336	endcase
337	;---------------------
338
339	terre = where(finite(res) EQ 0)
340	if terre[0] NE -1 then res[terre] = valmask
341
342	if n_elements(showbuild) then BEGIN
343	winsave = !window
344	psave = !p
345	xsave = !x
346	ysave = !y
347	plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '' $
348	, coast_thick = 2, window = showbuild
349	!p.title = ''
350	!p.subtitle = ''
351
352	plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50
353	plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2
354
355	FOR i = 0, n_elements(points1)-1 DO $
356	plots, [float(points1[i]), float(points2[i])] $
357	, [imaginary(points1[i]), imaginary(points2[i])], color = 150
358
359	plots, float(points1), imaginary(points1), color = 150, psym = 1
360	plots, float(points2), imaginary(points2), color = 150, psym = 1
361	plots, float(inter), imaginary(inter), color = 250, psym = 1
362
363	; ?? bug ?? IF terre[0] NE -1 THEN plots, float(terre[inter]), imaginary(terre[inter]), color = 0, psym = 1
364
365	; dummy = ''
366	; read, dummy, prompt = 'press return to continue'
367	IF !d.name EQ 'PS' THEN erase ELSE wset, winsave
368	!p = psave
369	!x = xsave
370	!y = ysave
371	ENDIF
372
373	restoreboxparam, 'boxparam4section.dat'
374
375	;---------------------
376	return
377	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format