Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 105

Last change on this file since 105 was 105, checked in by pinsard, 18 years ago
bugfix compute_fromreg_imoms3_weigaddr and square2quadrilateral
Property svn:executable set to ``*
File size: 16.5 KB

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

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

Download in other formats:

Original Format