Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 296

Last change on this file since 296 was 296, checked in by pinsard, 17 years ago
typo
Property svn:keywords set to `Id`
File size: 16.6 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 and 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 and 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, /double), 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, /double), 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, /double), 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, /double), 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