Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 125

Last change on this file since 125 was 125, checked in by pinsard, 18 years ago
improvements of Interpolation/*.pro header
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 16.3 KB

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

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

Download in other formats:

Original Format