Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 262

Last change on this file since 262 was 262, checked in by pinsard, 17 years ago
corrections of some headers and parameters and keywords case. change of pro2href to replace proidl
Property svn:keywords set to `Id`
File size: 16.5 KB

Line
1	;+
2	;
3	; @file_comments
4	; compute the weight and address need to interpolate data from a
5	; "regular grid" to any grid using the imoms3 method
6	;
7	; @categories
8	; Interpolation
9	;
10	; @param alonin {in}{required}{type=2d array}
11	; longitude of the input data
12	;
13	; @param alatin {in}{required}{type=2d array}
14	; latitude of the input data
15	;
16	; @param olonin {in}{required}{type=2d array}
17	; longitude of the output data
18	;
19	; @param olat {in}{required}{type=2d array}
20	; latitude of the output data
21	;
22	; @keyword NONORTHERNLINE {type=scalar 0 or 1}{default=0}
23	; put 1 if you don't want to take into
24	; account the northern line of the input data when performing the interpolation.
25	;
26	; @keyword NOSOUTHERNLINE {type=scalar 0 or 1}{default=0}
27	; put 1 if you don't want to take into
28	; account the southern line of the input data when performing the interpolation.
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, jpio, jpjo, /over)
38	;
39	; @restrictions
40	; - the input grid must be a "regular/rectangular grid", defined as a grid for
41	; which each longitude lines have the same latitude and each latitude columns
42	; have the same longitude.
43	; - We supposed the data are located on a sphere, with a periodicity along
44	; the longitude.
45	; - points located between the first/last 2 lines are interpolated
46	; using a imoms3 interpolation along the longitudinal direction and linear
47	; interpolation along the latitudinal direction
48	; - points located out of the southern and northern boundaries are interpolated
49	; using a imoms3 interpolation only along the longitudinal direction.
50	;
51	; @history
52	; November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)
53	; March 2006: works for rectangular grids
54	;
55	; @version
56	; $Id$
57	;
58	;-
59	;
60	PRO compute_fromreg_imoms3_weigaddr, alonin, alatin, olonin, olat, weig, addr $
61	, NONORTHERNLINE = nonorthernline, $
62	NOSOUTHERNLINE = nosouthernline
63	;
64	compile_opt idl2, strictarrsubs
65	;
66	alon = alonin
67	alat = alatin
68	olon = olonin
69	;
70	jpia = n_elements(alon)
71	jpja = n_elements(alat)
72	;
73	jpio = (size(olon, /dimensions))[0]
74	jpjo = (size(olon, /dimensions))[1]
75	;
76	; alon
77	minalon = min(alon, max = maxalon)
78	IF maxalon-minalon GE 360. THEN stop
79	; alon must be monotonically increasing
80	IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN
81	shiftx = -(where(alon EQ min(alon)))[0]
82	alon = shift(alon, shiftx)
83	IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop
84	ENDIF ELSE shiftx = 0
85	; alon is it regularly spaced?
86	step = alon-shift(alon, 1)
87	step[0] = step[0] + 360.
88	IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregx = 1
89	; we extend the longitude range of alon (-> easy interpolation even
90	; near minalon et maxalon)
91	toadd = 10*jpia/360+1
92	alon = [alon[jpia-toadd:jpia-1]-360., alon[*], alon[0:toadd-1]+360.]
93	jpia = jpia+2*toadd
94	; alat
95	revy = alat[0] GT alat[1]
96	IF revy THEN alat = reverse(alat)
97	; alat must be monotonically increasing
98	IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop
99	; alat is it regularly spaced?
100	step = alat-shift(alat, 1)
101	step = step[1:jpja - 1L]
102	IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregy = 1
103	;
104	if keyword_set(nonorthernline) then BEGIN
105	jpja = jpja - 1L
106	alat = alat[0: jpja-1L]
107	ENDIF
108	if keyword_set(nosouthernline) then BEGIN
109	alat = alat[1: jpja-1L]
110	jpja = jpja - 1L
111	ENDIF
112	; olon between minalon et minalon+360
113	out = where(olon LT minalon)
114	WHILE out[0] NE -1 DO BEGIN
115	olon[out] = olon[out]+360.
116	out = where(olon LT minalon)
117	ENDWHILE
118	out = where(olon GE minalon+360.)
119	WHILE out[0] NE -1 DO BEGIN
120	olon[out] = olon[out]- 360.
121	out = where(olon GE minalon+360.)
122	ENDWHILE
123	; make sure that all values of olon are located within values of alon
124	IF min(olon, max = ma) LT minalon THEN stop
125	IF ma GE minalon+360. THEN stop
126	;
127	xaddr = lonarr(16, jpio*jpjo)
128	yaddr = lonarr(16, jpio*jpjo)
129	weig = fltarr(16, jpio*jpjo)
130	;
131	indexlon = value_locate(alon, olon)
132	IF total(alon[indexlon] GT olon) NE 0 THEN stop
133	IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop
134	IF (where(indexlon LE 1L ))[0] NE -1 THEN stop
135	IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop
136	indexlat = value_locate(alat, olat)
137	;
138	; for the ocean points located below the atm line
139	; jpja-2 and above the line 1
140	; for those points we can always find 16 neighbors
141	; imoms interpolation along longitude and latitude
142	;
143	short = where(indexlat LT jpja-2L AND indexlat GE 1L)
144	ilon = indexlon[short]
145	ilat = indexlat[short]
146	;
147	IF NOT keyword_set(noregy) THEN BEGIN
148	delta = alat[ilat+1L]-alat[ilat]
149	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
150	delta = delta[0]
151	;
152	d0 = (alat[ilat-1L]-olat[short])/delta
153	IF min(d0, max = ma) LE -2 THEN stop
154	IF ma GT -1 THEN stop
155	wy0 = imoms3(temporary(d0))
156	d1 = (alat[ilat ]-olat[short])/delta
157	IF min(d1, max = ma) LE -1 THEN stop
158	IF ma GT 0 THEN stop
159	wy1 = imoms3(temporary(d1))
160	d2 = (alat[ilat+1L]-olat[short])/delta
161	IF min(d2, max = ma) LE 0 THEN stop
162	IF ma GT 1 THEN stop
163	wy2 = imoms3(temporary(d2))
164	d3 = (alat[ilat+2L]-olat[short])/delta
165	IF min(d3, max = ma) LE 1 THEN stop
166	IF ma GT 2 THEN stop
167	wy3 = imoms3(temporary(d3))
168	ENDIF ELSE BEGIN
169	nele = n_elements(short)
170	wy0 = fltarr(nele)
171	wy1 = fltarr(nele)
172	wy2 = fltarr(nele)
173	wy3 = fltarr(nele)
174	FOR i = 0L, nele-1 DO BEGIN
175	IF i MOD 10000 EQ 0 THEN print, i
176	newlat = spl_incr(alat[ilat[i]-1L:ilat[i]+2L], [-1., 0., 1., 2.], olat[short[i]])
177	IF newlat LE 0 THEN stop
178	IF newlat GT 1 THEN stop
179	wy0[i] = imoms3(newlat+1)
180	wy1[i] = imoms3(newlat)
181	wy2[i] = imoms3(1-newlat)
182	wy3[i] = imoms3(2-newlat)
183	ENDFOR
184	ENDELSE
185	;
186	mi = min(wy0+wy1+wy2+wy3, max = ma)
187	IF abs(mi-1) GE 1.e-6 THEN stop
188	IF abs(ma-1) GE 1.e-6 THEN stop
189	;
190	IF NOT keyword_set(noregx) THEN BEGIN
191	delta = alon[ilon]-alon[ilon-1L]
192	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
193	delta = delta[0]
194	;
195	d0 = (alon[ilon-1L]-olon[short])/delta
196	IF min(d0, max = ma) LE -2 THEN stop
197	IF ma GT -1 THEN stop
198	wx0 = imoms3(temporary(d0))
199	d1 = (alon[ilon ]-olon[short])/delta
200	IF min(d1, max = ma) LE -1 THEN stop
201	IF ma GT 0 THEN stop
202	wx1 = imoms3(temporary(d1))
203	d2 = (alon[ilon+1L]-olon[short])/delta
204	IF min(d2, max = ma) LE 0 THEN stop
205	IF ma GT 1 THEN stop
206	wx2 = imoms3(temporary(d2))
207	d3 = (alon[ilon+2L]-olon[short])/delta
208	IF min(d3, max = ma) LE 1 THEN stop
209	IF ma GT 2 THEN stop
210	wx3 = imoms3(temporary(d3))
211	ENDIF ELSE BEGIN
212	nele = n_elements(short)
213	wx0 = fltarr(nele)
214	wx1 = fltarr(nele)
215	wx2 = fltarr(nele)
216	wx3 = fltarr(nele)
217	FOR i = 0L, nele-1 DO BEGIN
218	IF i MOD 10000 EQ 0 THEN print, i
219	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
220	IF newlon LE 0 THEN stop
221	IF newlon GT 1 THEN stop
222	wx0[i] = imoms3(newlon+1)
223	wx1[i] = imoms3(newlon)
224	wx2[i] = imoms3(1-newlon)
225	wx3[i] = imoms3(2-newlon)
226	ENDFOR
227	ENDELSE
228	;
229	mi = min(wx0+wx1+wx2+wx3, max = ma)
230	IF abs(mi-1) GE 1.e-6 THEN stop
231	IF abs(ma-1) GE 1.e-6 THEN stop
232	;
233	; line 0
234	xaddr[0, short] = ilon - 1L
235	xaddr[1, short] = ilon
236	xaddr[2, short] = ilon + 1L
237	xaddr[3, short] = ilon + 2L
238	yaddr[0, short] = ilat - 1L
239	yaddr[1, short] = yaddr[0, short]
240	yaddr[2, short] = yaddr[0, short]
241	yaddr[3, short] = yaddr[0, short]
242	weig[0, short] = wx0 * wy0
243	weig[1, short] = wx1 * wy0
244	weig[2, short] = wx2 * wy0
245	weig[3, short] = wx3 * wy0
246	; line 1
247	xaddr[4, short] = ilon - 1L
248	xaddr[5, short] = ilon
249	xaddr[6, short] = ilon + 1L
250	xaddr[7, short] = ilon + 2L
251	yaddr[4, short] = ilat
252	yaddr[5, short] = ilat
253	yaddr[6, short] = ilat
254	yaddr[7, short] = ilat
255	weig[4, short] = wx0 * wy1
256	weig[5, short] = wx1 * wy1
257	weig[6, short] = wx2 * wy1
258	weig[7, short] = wx3 * wy1
259	; line 2
260	xaddr[8, short] = ilon - 1L
261	xaddr[9, short] = ilon
262	xaddr[10, short] = ilon + 1L
263	xaddr[11, short] = ilon + 2L
264	yaddr[8, short] = ilat + 1L
265	yaddr[9, short] = yaddr[8, short]
266	yaddr[10, short] = yaddr[8, short]
267	yaddr[11, short] = yaddr[8, short]
268	weig[8, short] = wx0 * wy2
269	weig[9, short] = wx1 * wy2
270	weig[10, short] = wx2 * wy2
271	weig[11, short] = wx3 * wy2
272	; line 3
273	xaddr[12, short] = ilon - 1L
274	xaddr[13, short] = ilon
275	xaddr[14, short] = ilon + 1L
276	xaddr[15, short] = ilon + 2L
277	yaddr[12, short] = ilat + 2L
278	yaddr[13, short] = yaddr[12, short]
279	yaddr[14, short] = yaddr[12, short]
280	yaddr[15, short] = yaddr[12, short]
281	weig[12, short] = wx0 * wy3
282	weig[13, short] = wx1 * wy3
283	weig[14, short] = wx2 * wy3
284	weig[15, short] = wx3 * wy3
285	;
286	mi = min(total(weig[*, short], 1), max = ma)
287	IF abs(mi-1) GE 1.e-6 THEN stop
288	IF abs(ma-1) GE 1.e-6 THEN stop
289	;
290	; for the ocean points located between the atm lines
291	; jpja-2 and jpja-1 or between the atm lines 0 and 1
292	; linear interpolation between line 1 and line 2
293	;
294	short = where(indexlat EQ jpja-2L OR indexlat EQ 0)
295	IF short[0] NE -1 THEN BEGIN
296	ilon = indexlon[short]
297	ilat = indexlat[short]
298	;
299	delta = alat[ilat+1L]-alat[ilat]
300	IF NOT keyword_set(noregy) THEN BEGIN
301	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
302	delta = delta[0]
303	ENDIF
304	;
305	d1 = (alat[ilat ]-olat[short])/delta
306	IF min(d1, max = ma) LE -1 THEN stop
307	IF ma GT 0 THEN stop
308	wy1 = 1.+ temporary(d1)
309	d2 = (alat[ilat+1L]-olat[short])/delta
310	IF min(d2, max = ma) LE 0 THEN stop
311	IF ma GT 1 THEN stop
312	wy2 = 1.- temporary(d2)
313	;
314	mi = min(wy1+wy2, max = ma)
315	IF abs(mi-1) GE 1.e-6 THEN stop
316	IF abs(ma-1) GE 1.e-6 THEN stop
317	; but imoms3 along the longitude
318	IF NOT keyword_set(noregx) THEN BEGIN
319	delta = alon[ilon]-alon[ilon-1L]
320	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
321	delta = delta[0]
322	;
323	d0 = (alon[ilon-1L]-olon[short])/delta
324	IF min(d0, max = ma) LE -2 THEN stop
325	IF ma GT -1 THEN stop
326	wx0 = imoms3(temporary(d0))
327	d1 = (alon[ilon ]-olon[short])/delta
328	IF min(d1, max = ma) LE -1 THEN stop
329	IF ma GT 0 THEN stop
330	wx1 = imoms3(temporary(d1))
331	d2 = (alon[ilon+1L]-olon[short])/delta
332	IF min(d2, max = ma) LE 0 THEN stop
333	IF ma GT 1 THEN stop
334	wx2 = imoms3(temporary(d2))
335	d3 = (alon[ilon+2L]-olon[short])/delta
336	IF min(d3, max = ma) LE 1 THEN stop
337	IF ma GT 2 THEN stop
338	wx3 = imoms3(temporary(d3))
339	ENDIF ELSE BEGIN
340	nele = n_elements(short)
341	wx0 = fltarr(nele)
342	wx1 = fltarr(nele)
343	wx2 = fltarr(nele)
344	wx3 = fltarr(nele)
345	FOR i = 0L, nele-1 DO BEGIN
346	IF i MOD 10000 EQ 0 THEN print, i
347	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
348	IF newlon LE 0 THEN stop
349	IF newlon GT 1 THEN stop
350	wx0[i] = imoms3(newlon+1)
351	wx1[i] = imoms3(newlon)
352	wx2[i] = imoms3(1-newlon)
353	wx3[i] = imoms3(2-newlon)
354	ENDFOR
355	ENDELSE
356	;
357	mi = min(wx0+wx1+wx2+wx3, max = ma)
358	IF abs(mi-1) GE 1.e-6 THEN stop
359	IF abs(ma-1) GE 1.e-6 THEN stop
360	; line 1
361	xaddr[0, short] = ilon - 1L
362	xaddr[1, short] = ilon
363	xaddr[2, short] = ilon + 1L
364	xaddr[3, short] = ilon + 2L
365	yaddr[0, short] = ilat
366	yaddr[1, short] = ilat
367	yaddr[2, short] = ilat
368	yaddr[3, short] = ilat
369	weig[0, short] = wx0 * wy1
370	weig[1, short] = wx1 * wy1
371	weig[2, short] = wx2 * wy1
372	weig[3, short] = wx3 * wy1
373	; line 2
374	xaddr[4, short] = ilon - 1L
375	xaddr[5, short] = ilon
376	xaddr[6, short] = ilon + 1L
377	xaddr[7, short] = ilon + 2L
378	yaddr[4, short] = ilat + 1L
379	yaddr[5, short] = yaddr[4, short]
380	yaddr[6, short] = yaddr[4, short]
381	yaddr[7, short] = yaddr[4, short]
382	weig[4, short] = wx0 * wy2
383	weig[5, short] = wx1 * wy2
384	weig[6, short] = wx2 * wy2
385	weig[7, short] = wx3 * wy2
386	;
387	mi = min(total(weig[*, short], 1), max = ma)
388	IF abs(mi-1) GE 1.e-6 THEN stop
389	IF abs(ma-1) GE 1.e-6 THEN stop
390	;
391	ENDIF
392	;
393	; for the ocean points located below the line 0
394	; Interpolation only along the longitude
395	;
396	short = where(indexlat EQ -1)
397	IF short[0] NE -1 THEN BEGIN
398	ilon = indexlon[short]
399	;
400	IF NOT keyword_set(noregx) THEN BEGIN
401	delta = alon[ilon]-alon[ilon-1L]
402	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
403	delta = delta[0]
404	;
405	d0 = (alon[ilon-1L]-olon[short])/delta
406	IF min(d0, max = ma) LE -2 THEN stop
407	IF ma GT -1 THEN stop
408	wx0 = imoms3(temporary(d0))
409	d1 = (alon[ilon ]-olon[short])/delta
410	IF min(d1, max = ma) LE -1 THEN stop
411	IF ma GT 0 THEN stop
412	wx1 = imoms3(temporary(d1))
413	d2 = (alon[ilon+1L]-olon[short])/delta
414	IF min(d2, max = ma) LE 0 THEN stop
415	IF ma GT 1 THEN stop
416	wx2 = imoms3(temporary(d2))
417	d3 = (alon[ilon+2L]-olon[short])/delta
418	IF min(d3, max = ma) LE 1 THEN stop
419	IF ma GT 2 THEN stop
420	wx3 = imoms3(temporary(d3))
421	ENDIF ELSE BEGIN
422	nele = n_elements(short)
423	wx0 = fltarr(nele)
424	wx1 = fltarr(nele)
425	wx2 = fltarr(nele)
426	wx3 = fltarr(nele)
427	FOR i = 0L, nele-1 DO BEGIN
428	IF i MOD 10000 EQ 0 THEN print, i
429	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
430	IF newlon LE 0 THEN stop
431	IF newlon GT 1 THEN stop
432	wx0[i] = imoms3(newlon+1)
433	wx1[i] = imoms3(newlon)
434	wx2[i] = imoms3(1-newlon)
435	wx3[i] = imoms3(2-newlon)
436	ENDFOR
437	ENDELSE
438	;
439	mi = min(wx0+wx1+wx2+wx3, max = ma)
440	IF abs(mi-1) GE 1.e-6 THEN stop
441	IF abs(ma-1) GE 1.e-6 THEN stop
442	; line 1
443	xaddr[0, short] = ilon - 1L
444	xaddr[1, short] = ilon
445	xaddr[2, short] = ilon + 1L
446	xaddr[3, short] = ilon + 2L
447	yaddr[0:3, short] = 0
448	weig[0, short] = wx0
449	weig[1, short] = wx1
450	weig[2, short] = wx2
451	weig[3, short] = wx3
452	;
453	mi = min(total(weig[*, short], 1), max = ma)
454	IF abs(mi-1) GE 1.e-6 THEN stop
455	IF abs(ma-1) GE 1.e-6 THEN stop
456	;
457	ENDIF
458	;
459	; for the ocean points located above jpia-1
460	; Interpolation only along the longitude
461	;
462	short = where(indexlat EQ jpja-1L)
463	IF short[0] NE -1 THEN BEGIN
464	ilon = indexlon[short]
465	;
466	IF NOT keyword_set(noregx) THEN BEGIN
467	delta = alon[ilon]-alon[ilon-1L]
468	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
469	delta = delta[0]
470	;
471	d0 = (alon[ilon-1L]-olon[short])/delta
472	IF min(d0, max = ma) LE -2 THEN stop
473	IF ma GT -1 THEN stop
474	wx0 = imoms3(temporary(d0))
475	d1 = (alon[ilon ]-olon[short])/delta
476	IF min(d1, max = ma) LE -1 THEN stop
477	IF ma GT 0 THEN stop
478	wx1 = imoms3(temporary(d1))
479	d2 = (alon[ilon+1L]-olon[short])/delta
480	IF min(d2, max = ma) LE 0 THEN stop
481	IF ma GT 1 THEN stop
482	wx2 = imoms3(temporary(d2))
483	d3 = (alon[ilon+2L]-olon[short])/delta
484	IF min(d3, max = ma) LE 1 THEN stop
485	IF ma GT 2 THEN stop
486	wx3 = imoms3(temporary(d3))
487	ENDIF ELSE BEGIN
488	nele = n_elements(short)
489	wx0 = fltarr(nele)
490	wx1 = fltarr(nele)
491	wx2 = fltarr(nele)
492	wx3 = fltarr(nele)
493	FOR i = 0L, nele-1 DO BEGIN
494	IF i MOD 10000 EQ 0 THEN print, i
495	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
496	IF newlon LE 0 THEN stop
497	IF newlon GT 1 THEN stop
498	wx0[i] = imoms3(newlon+1)
499	wx1[i] = imoms3(newlon)
500	wx2[i] = imoms3(1-newlon)
501	wx3[i] = imoms3(2-newlon)
502	ENDFOR
503	ENDELSE
504	;
505	mi = min(wx0+wx1+wx2+wx3, max = ma)
506	IF abs(mi-1) GE 1.e-6 THEN stop
507	IF abs(ma-1) GE 1.e-6 THEN stop
508	; line 1
509	xaddr[0, short] = ilon-1L
510	xaddr[1, short] = ilon
511	xaddr[2, short] = ilon+1L
512	xaddr[3, short] = ilon+2L
513	yaddr[0:3, short] = jpja-1L
514	weig[0, short] = wx0
515	weig[1, short] = wx1
516	weig[2, short] = wx2
517	weig[3, short] = wx3
518	;
519	mi = min(total(weig[*, short], 1), max = ma)
520	IF abs(mi-1) GE 1.e-6 THEN stop
521	IF abs(ma-1) GE 1.e-6 THEN stop
522	;
523	ENDIF
524	;
525	; Come back to the original index of atm grid without longitudinal overlap.
526	;
527	;
528	xaddr = temporary(xaddr) - toadd
529	jpia = jpia - 2*toadd
530	; make sure all values are ge 0
531	xaddr = temporary(xaddr) + jpia
532	; range the values between 0 and jpia-1
533	xaddr = temporary(xaddr) mod jpia
534	;
535	; take into account shiftx if needed
536	IF shiftx NE 0 THEN xaddr = (temporary(xaddr) - shiftx) MOD jpia
537	; take into account nosouthernline and nonorthernline
538	if keyword_set(nosouthernline) then BEGIN
539	yaddr = temporary(yaddr) + 1L
540	jpja = jpja + 1L
541	ENDIF
542	if keyword_set(nonorthernline) then jpja = jpja + 1L
543	; take into account revy if needed
544	IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr)
545	; ;
546	addr = temporary(yaddr)*jpia+temporary(xaddr)
547	;
548	RETURN
549	END

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

Download in other formats:

Original Format