Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 72

Last change on this file since 72 was 69, checked in by smasson, 18 years ago
debug + new xxx
Property svn:executable set to ``*
File size: 16.6 KB

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

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

Download in other formats:

Original Format