Context Navigation

extract_amsu.pro @ 545

Last change on this file since 545 was 545, checked in by pinsard, 12 years ago
fix for doc
File size: 17.9 KB

Line
1	; docformat = 'rst'
2	;+
3	;
4	; .. _extract_amsu.pro:
5	;
6	; =================
7	; extract_amsu.pro
8	; =================
9	;
10	; DESCRIPTION
11	; ===========
12	;
13	; prgm de lecture des fichiers AMSU (A et B)
14	;
15	; decode les noms des fichiers donnés en argument dans files_list
16	; dans la date choisie, puis appelle le prgm de lecture
17	;
18	; appelle :ref:`interpol_correc.pro`, qui fournit les fonctions de correction au
19	; nadir et les applique, en fonction du masque terre-mer issu de la
20	; carte de bathymetrie de S. Masson
21	;
22	; applique une interpolation tenant compte de la dimension de la tache
23	; au sol selon la position dans la fauchee
24	;
25	; ecrit les donnees pour la zone (lat/long) choisie
26	;
27	; ++ le fichier de sortie contient dans l'ordre mm, jj (jour de
28	; l'annee), le temps (identique pour les points d'une meme fauchee), le
29	; no satellite, le no de fov (1 a 40 ou 1 a 90), long, lat, tb
30	;
31	; Les valeurs manquantes sont codées par NaN.
32	;
33	; .. only:: man
34	;
35	; Figure is visible on PDF and HTML documents only.
36	;
37	; .. only:: html or latex
38	;
39	; .. graphviz::
40	;
41	; digraph extract_amsu {
42	;
43	; amsu_t1 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/AMSU/AMSU/L1C/yyyy/yyyy_mm_dd/.L1C"];
44	; amsu_t2 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/yyyy/mm/numch_yyyymmdd_geomin_geomax.dat"];
45	;
46	; extract_amsua [shape=box,
47	; fontname=Courier,
48	; color=blue,
49	; URL="http://forge.ipsl.jussieu.fr/varamma/browser/trunk/src/extract_amsua.pro",
50	; label="${PROJECT}/src/extract_amsua.pro"];
51	;
52	; {amsu_t1} -> {extract_amsua} -> {amsu_t2};
53	; }
54	;
55	; :param numch:
56	; :param files_list:
57	; :param yyyy:
58	; :param mm:
59	; :param resol:
60	; :param lon_min:
61	; :param lon_max:
62	; :param lat_min:
63	; :param lat_max:
64	;
65	; EXAMPLES
66	; ========
67	;
68	; Creating an articial list of AMSU-A files::
69	;
70	; files_list = strarr(2)
71	; files_list[0] = project_id_env + '/AMSU/AMSUAN15/L1C/2006/2006_08_01/NSS.AMAX.NK.D06213.S0110.E0255.B4271112.GC.L1C'
72	; files_list[1] = project_id_env + 'AMSU/AMSUAN15/L1C/2006/2006_08_01/NSS.AMAX.NK.D06213.S0250.E0437.B4271213.GC.L1C'
73	;
74	; This above tip can replace the call to :func:`search_amsufiles`.
75	;
76	; Pour ne garder que le dernier fichier::
77	;
78	; files_list = search_amsufiles(numch, yyyy, mm, dd)
79	; help, files_list
80	; files_list_last_array = strarr(1)
81	; files_list_last_array[0] = files_list[N_ELEMENTS(files_list) - 1]
82	; help, files_list_last_array
83	; files_list = files_list_last_array
84	; help, files_list
85	;
86	; Using AMSU-A channel 5::
87	;
88	; numch = 'a5'
89	; yyyy=2006L
90	; mm=8
91	; dd=13
92	; resol=1
93	; lon_min=-25.
94	; lon_max=25.
95	; lat_min=-5.
96	; lat_max=20.
97	; files_list = search_amsufiles(numch, yyyy, mm, dd)
98	; extract_amsua, numch, files_list, yyyy, mm, dd, resol, lon_min, lon_max, lat_min, lat_max
99	;
100	; :file:`${PROJECT_ID}/AMSU/2006/08/a5_20060813_025w05s_025e20n.dat`
101	; must have been created.
102	;
103	; SEE ALSO
104	; ========
105	;
106	; :ref:`data_amsu`
107	;
108	; called by :ref:`traite_amsuab.sh`
109	;
110	; previous step : :func:`search_amsufiles`
111	;
112	; use :ref:`read_amsua1c.pro`, :ref:`read_amsub1c.pro`,
113	; :func:`geolocation_to_string_idl`, :func:`mem_to_file_amsu_t2`
114	;
115	; next step : ++
116	;
117	; TODO
118	; ====
119	;
120	; track why end with::
121	;
122	; % Program caused arithmetic error: Floating illegal operand
123	;
124	; traiter les longitudes autour de 180 degres (passage de 180 a -180)
125	; interpolswath ne gere pas cette situation
126	;
127	; lever le doute sur le contenu du fichier écrit par le printf (je (fplod)
128	; crains des zeros louches!!).louche
129	;
130	;
131	; ajouter MHS
132	;
133	; eviter la creation d'un fichier vide (si pas de points dans la zone)
134	;
135	; tester la cohérence entre les fichiers présents dans files_list
136	; en terme de date avec yyyy, mm et dd donnés en paramètre
137	;
138	; improve log
139	;
140	; get rid of uppercase
141	;
142	; check args
143	;
144	; decrire la limitation avec les seuils
145	;
146	; EVOLUTIONS
147	; ==========
148	;
149	;
150	; $Id: extract_amsua.pro 513 2012-04-10 14:50:06Z pinsard $
151	;
152	; $URL: svn+ssh://lelod@forge.ipsl.jussieu.fr/ipsl/forge/projets/varamma/svn/trunk/src/extract_amsua.pro $
153	;
154	; - lelod 20120420
155	;
156	; * extension de extract_amsua à amsub
157	;
158	; - lelod 20120418
159	;
160	; * neutralisation de la variable desc: remplacement par numero de l'orbite lue
161	; - fplod 20120323
162	;
163	; * homogenize examples to traite_amsuab.sh
164	; * fix output name in example
165	; * add a paramater files_list to remplace usage of :file:`list_file`
166	; * replace usage of :file:`list_file` by the usage of the string array
167	; files_list
168	; * add a test if files_list no empty (see nofile label)
169	; * reorder SEE ALSO and TODO
170	;
171	; - lelod 20120223
172	;
173	; * modif lecture du fichier bathy (ETOPO) pour limiter le temps
174	; necessaire. Introduction de coordonnees limite dans initncdf, puis
175	; domdef
176	;
177	; Attention, si on change la zone, il faut changer les coordonnees!!!
178	;
179	; - lelod 20111215
180	;
181	; * tests de fonctionnement effectues: la correction nadir est
182	; correctement implementee
183	;
184	; l'interpolation de la fauchee fonctionne
185	; mais attention test sur la longitude code en dur:
186	; pbs avec les fauchees qui traversent le meridien 180deg!!!!
187	; l'interpolation ajoute des points parasites dans la zone...
188	;
189	;
190	; - fplod 20111214T082300Z aedon.locean-ipsl.upmc.fr (Darwin)
191	;
192	; * remove old style writing output file
193	; * forget idldoc
194	;
195	; - fplod 20111213T092550Z aedon.locean-ipsl.upmc.fr (Darwin)
196	;
197	; * fix syntax error : rewrite labfile label
198	; * fix nocanal string to integer type
199	; * fix example (missing resol in call)
200	;
201	; - lelod 20111211
202	;
203	; * reintegre les modules de correct_nadir_amsu dans extract
204	; prevu decodage de numch pour connaitre le nom instrument et le
205	; numero du canal (pas fini, au debut) --> a faire aussi dans les
206	; ssprgm appeles (interpol_correc)
207	; espere que c'est le dernier avatar de cette !!!!?XX^! de
208	; chaine de traitement AMSU
209	;
210	; - fplod 20111209T125312Z aedon.locean-ipsl.upmc.fr (Darwin)
211	;
212	; * chgt de terminologie
213	; * ajout parametre numch
214	; * transition pour utilisation de :ref:`mem_to_file_amsu_t2`
215	;
216	; - lelod 20111121
217	;
218	; * introduction du parametre "descending" (orbite descendante 1,
219	; montante 0) dans le fichier de sortie : a tester!
220	;
221	; - fplod 20110818T160435Z aedon.locean-ipsl.upmc.fr (Darwin)
222	;
223	; * no more /bdd in examples
224	; * need to check if missing exist before NaN affection on idl7 cratos
225	;
226	; - pinsard 2011-08-18T10:34:53Z loholt1.ipsl.polytechnique.fr (Linux)
227	;
228	; * correction of examples
229	; * examples with 2 files in list_filea and list_fileb
230	; * match technique from
231	; http://cerebus.as.arizona.edu/~ioannis/teaching/idl/tutorial_strings.html
232	; * revision du cas use_amsua=1 : la liste de fichiers amsu-a n'est ouverte
233	; qu'une seule fois
234	;
235	; - pinsard 2011-08-18T09:10:08Z loholt1.ipsl.polytechnique.fr (Linux)
236	;
237	; * move bad values from -9999.99 to NaN
238	;
239	; - pinsard 2011-08-17T14:58:13Z loholt1.ipsl.polytechnique.fr (Linux)
240	;
241	; * use geolocation_to_string_idl
242	;
243	; - pinsard 2011-08-10T10:46:37Z loholt1.ipsl.polytechnique.fr (Linux)
244	;
245	; * remove blank in output filename
246	; * revision of articifial example with use_amsua=0
247	; * add articifial example with use_amsua=1 and try (but fail) to
248	; make it work on loholt1
249	; * remove length path dependencies using IDL function file_basename
250	; * use same date and box of examples in correct_nadir_amsu.pro
251	;
252	; - fplod 20110810T085644Z aedon.locean-ipsl.upmc.fr (Darwin)
253	;
254	; * revision of header
255	;
256	; - pinsard 2011-05-26T08:44:38Z loholt1.ipsl.polytechnique.fr (Linux)
257	;
258	; * strmid(fich,51) is not anymore pertinent to isolate the correct information
259	;
260	; It was ok when path as /bdd/AMSU-1C/AMSUAN15/L1C/2006/2006_07_28/::
261	;
262	; c='/bdd/AMSU-1C/AMSUAN15/L1C/2006/2006_07_28/NSS.AMAX.NK.D06209.S1323.E1455.B4266162.GC.L1C'
263	; b=strmid(c,51)
264	; print, b
265	; NK.D06209.S1323.E1455.B4266162.GC.L1C
266	;
267	; but now path can be anything (PROJECT_ID environment variable)
268	; so the correct way is to catch the last 36 characters
269	;
270	; test sur skylla: fonctionne au moins pour un fichier!
271	;
272	;-
273	PRO extract_amsu, numch, files_list, yyyy, mm, dd, resol $
274	, lon_min, lon_max, lat_min, lat_max
275	;
276	compile_opt idl2, strictarrsubs
277	;
278	@cm_project
279	@common
280
281	; test if some files to read
282	CASE size(files_list,/DIMENSION) OF
283	0L : BEGIN
284	print, 'pas de fichiers'
285	goto, nofile
286	END
287	ELSE: BEGIN
288	print, 'nb de fichers', size(files_list,/DIMENSION)
289	END
290	ENDCASE
291	nb_file_array = size(files_list,/DIMENSION)
292	nb_file = (nb_file_array)[0]
293	;
294	nomcanal=strmid(numch,0,1)
295	nocanal=0
296	reads, strmid(numch,1,1),nocanal,format='(i1.1)'
297
298	tbmin=100
299	tbmax=350
300	PRINT, 'www : debut programme extract_amsu',SYSTIME()
301
302	jpie = n_elements(xxe)
303	jpje = n_elements(yye)
304
305	; appel au ssprgm qui interpole les fonctions de correction (inutile
306	; de l'appeler a chaque fichier....)
307	interpol_correc,numch,nbpixel,cor_l,cor_s,swath,track
308	; ouverture des fichiers liste (annee, mois, jour, tous satellites) pour
309	; chaque instrument
310	; boucle sur les elements de la liste
311	desc=0
312	print, 'iii : traitement du jour ',yyyy,mm,dd
313	a = STRARR(nb_file)
314	;PRINT,'demarrage boucle sur lichiers', SYSTIME()
315	FOR ifile = 0, nb_file - 1 do begin
316	filea = files_list[ifile]
317	if nomcanal eq 'a' then begin
318	COMMON amsua_header,ama_head
319	COMMON amsua_data ,ama_scan
320	print, 'ouverture et lecture du fichier ', filea, SYSTIME()
321	openr,lu1,filea,Error=erra,/get_lun
322	read_amsua1c,filea, flag1
323	close,lu1
324	free_lun,lu1
325	if (flag1 eq 0) then goto, labfile
326	na=SIZE(reform(ama_scan.btemps[0,,]))
327	na_scan=na[2] ; nb lignes (fauchees) dans l'orbite
328	n_scan=na_scan
329	na_fov=na[1] ; nb points dans la fauchee
330	n=na
331	nbpix= na_fov
332	n_fov=na_fov
333	nosat=ama_head.h_satid
334	; on lit la structure et on extrait les infos: temps, longitude,
335	; latitude, et les Tbs des differents canaux (15 pour AMSUA)
336	; extraction du canal traite dans un domaine enveloppant la zone choisie
337	amalong=REFORM(ama_scan.latlon[1,,]/1.E4)
338	amalati=REFORM(ama_scan.latlon[0,,]/1.E4)
339	amcha=REFORM(ama_scan.btemps[nocanal,,]/100.)
340	amzeni=REFORM(ama_scan.angles[0,,]/100.)
341	ttt=REFORM(ama_scan.scnlintime/3600000.)
342	amaday=REFORM(ama_scan.scnlindy)
343	endif
344	if nomcanal eq 'b' then begin
345	COMMON amsub_header,amb_head
346	COMMON amsub_data ,amb_scan
347	print, 'ouverture et lecture du fichier ', fileb, SYSTIME()
348	openr,lu1,fileb,Error=errb,/get_lun
349	read_amsub1c,fileb, flag2
350	close,lu1
351	free_lun,lu1
352	if (flag2 eq 0) then goto, labfile
353	nb=SIZE(reform(amb_scan.btemps[0,,]))
354	nb_scan=nb[2] ; nb lignes (ou fauchees)
355	n_scan=nb_scan
356	nb_fov=nb[1] ; nb pixels dans la fauchee
357	nbpix = nb_fov ; nb pixels dans la fauchee AMSUB
358	nosat=amb_head.h_satid ; on le lit dans le fichier, donc pas besoin de l'info dans la liste
359	amalong=REFORM(amb_scan.latlon[1,,]/1.E4)
360	amalati=REFORM(amb_scan.latlon[0,,]/1.E4)
361	amcha=REFORM(amb_scan.btemps[nocanal,,]/100.)
362	amzeni=REFORM(amb_scan.angles[0,,]/100.)
363	ttt=REFORM(amb_scan.scnlintime/3600000.)
364	amaday=REFORM(amb_scan.scnlindy)
365	endif
366
367	fov=indgen(nbpix) ; nb pixels dans la fauchee AMSU
368	ntime_1=fltarr(nbpix,n_scan)
369	fovy_1=fltarr(nbpix,n_scan)
370	for j=0L, nbpix-1L do begin
371	fovy_1[j,*]=j+1 ; on remplit un tableau de nscan lignes avec les valeurs de no de pixel
372	endfor
373	midpix=fix(nbpix/2)
374	jnd=where(amalong[midpix,] gt lon_min-15 and amalong[midpix,] lt lon_max+15 and amalati[midpix,] gt lat_min-15 and amalati[midpix,] lt lat_max+15 ,nzon)
375	print,"iii : nb de points du fichier dans le domaine geographique +15deg ",nzon, SYSTIME()
376	if nzon gt 1 then begin
377	amalat=amalati[*,jnd]
378	amalon=amalong[*,jnd]
379	amch=amcha[*,jnd]
380	amzen=amzeni[*,jnd]
381	tt=ttt[jnd]
382	dims = SIZE(amalon, /DIMENSIONS)
383	print,'dimension tableaux extraits',dims
384	n_scan=dims[1]
385	PRINT,'il y a des donnees dans la region '
386	amafov= fovy_1
387	jour_ama=amaday[0] ; jour du passage
388
389	;
390	; correction nadir des donnees
391	ch_nadir=fltarr(nbpix,nzon)
392	landseamask=intarr(nbpix,nzon)+2 ; valeur hors zone selectionnee
393	; PRINT,'boucle sur les points du fichier, correction nadir', SYSTIME()
394	for isc=0L,nzon-1L do begin
395	;recherche de la zone xxe,yye englobant la fauchee
396	utilex=where(xxe ge amalon[midpix,isc]-15 and xxe le amalon[midpix,isc]+15,nxu)
397	if nxu ne 0 then begin
398	xxutile=xxe[utilex]
399	batutilx=bate[utilex,*]
400	utiley=where(yye ge amalat[midpix,isc]-8 and yye le amalat[midpix,isc]+8,nyu)
401	if nyu ne 0 then begin
402	yyutile=yye[utiley]
403	batutil=batutilx[*,utiley]
404
405	for ifo=0,nbpix-1 do begin
406	; approximation: 1km= 1deg/100, et on considere le quart en lon et lat
407	; de la surface du pixel pour optimiser la localisation du centre du
408	; pixel attention approximation brutale des distances en fractions dedegre
409	; recherche de la cote par rapport a la resolution AMSUA
410	xind=where (xxutile ge amalon[ifo,isc]-swath[ifo]/400. and xxutile le amalon[ifo,isc]+swath[ifo]/400.,nxlandsea)
411	yind=where (yyutile ge amalat[ifo,isc]-track[ifo]/400. and yyutile le amalat[ifo,isc]+track[ifo]/400.,nylandsea)
412	if(nxlandsea ne 0 and nylandsea ne 0) then begin
413	bate1=batutil[xind,*]
414	bate2=bate1[*,yind]
415	if (mean(bate2) le 0.5) then begin
416	ch_nadir[ifo,isc]=amch[ifo,isc]-cor_s[ifo]
417	landseamask[ifo,isc]=0 ; mer
418	endif else begin
419	ch_nadir[ifo,isc]=amch[ifo,isc]-cor_l[ifo]
420	landseamask[ifo,isc]=1 ; terre
421	endelse
422	endif
423	endfor
424	endif
425	endif
426
427	endfor
428	; PRINT, 'fin boucle ',SYSTIME()
429	moych=fltarr(nbpix)
430
431	; appel a interpolswath pour ajuster les pixels amsua sur une grille
432	; reguliere dans la fauchee
433	; et selection de la zone conservee
434	; attention test sur la longitude code en dur:
435	; pbs avec les fauchees qui traversent le meridien 180deg!!!!
436	; l'interpolation ajoute des points parasites dans la zone...
437	nbgrid=0
438	cont=0L
439	chint=fltarr(1)
440	for i=0L,nzon-1L do begin
441	tb=ch_nadir[*,i]
442	lon=amalon[*,i]
443	lat=amalat[*,i]
444	mask=landseamask[*,i]
445	ind=where(tb gt tbmin,nbon)
446	; test pour eviter pb d'interpolation de longitude (meridien 180)
447	if nbon eq nbpix then begin
448	cont=cont+1L
449	nbgrid1=nbgrid
450	interpolswath,tb,lat,lon,mask,resol,nbgrid,tbgrid,latgrid,longrid,maskgrid
451	fovgrid=indgen(n_elements(tbgrid))+1
452	;print,'wwwtest: nb points fov interpole ',n_elements(tbgrid)
453	; selection des taches au sol situees dans la zone d'interet
454	zone=where((longrid ge lon_min) and (longrid le lon_max) $
455	and (latgrid ge lat_min) and (latgrid le lat_max) and (tbgrid gt tbmin) and (tbgrid lt tbmax), npt)
456	;print,'wwwtest: nb points zone d interet ',npt
457	if npt ne 0 then begin
458	if (n_elements(chint) LE 1) then begin
459	timeint=replicate(tt[i],npt)
460	chint=tbgrid[zone]
461	latint=latgrid[zone]
462	lonint=longrid[zone]
463	fovint=fovgrid[zone]
464	maskint=maskgrid[zone]
465	endif else begin
466	chint=[chint,tbgrid[zone]]
467	latint=[latint,latgrid[zone]]
468	lonint=[lonint,longrid[zone]]
469	fovint=[fovint,fovgrid[zone]]
470	maskint=[maskint,maskgrid[zone]]
471	timeint=[timeint,replicate(tt[i],npt)]
472	endelse
473	endif else begin
474	;print, 'www : aucun point dans la zone'
475	endelse
476	endif
477	endfor
478	; PRINT,'fin interpolation et selection des points de la zone ', SYSTIME()
479
480	nn=n_elements(chint)
481
482	print,'www : nb points dans la zone en fin de traitement de l''orbite ',nn
483	if nn gt 1 then begin
484	for i=0L, nn-1L do begin
485	if (n_elements(data) EQ 0) then begin
486	data = { desc : desc $
487	, hour : timeint[i] $
488	, fov : fovint[i] $
489	, lon : lonint[i] $
490	, lat: latint[i] $
491	, landseamask: maskint[i] $
492	, tb : chint[i] $
493	}
494	endif else begin
495	data = [data , { desc : desc $
496	, hour : timeint[i] $
497	, fov : fovint[i] $
498	, lon : lonint[i] $
499	, lat: latint[i] $
500	, landseamask: maskint[i] $
501	, tb : chint[i] $
502	}]
503	endelse
504	endfor
505	desc=desc+1 ; incrementation du numero de l'orbite
506	endif
507	;
508	; fin boucle sur les fichiers lus
509	endif
510	PRINT,'www : passage au fichier suivant ', SYSTIME()
511	endfor
512	labfile:
513	if (n_elements(data) NE 0) then begin
514	; write outputfile
515	header1 = { yyyy : yyyy $
516	, mm : mm $
517	, dd : dd $
518	, numch : numch $
519	, nbpix : nbgrid $
520	, resol : resol $
521	}
522	header2= { lon_min: lon_min $
523	, lon_max: lon_max $
524	, lat_min: lat_min $
525	, lat_max: lat_max $
526	}
527
528	amsu_t2={header1: header1, header2: header2, data: data}
529	help, amsu_t2,/structure
530	testfilename=''
531	fileout = mem_to_file_amsu_t2(amsu_t2, testfilename)
532
533	endif else begin
534	print, 'www : no data to write'
535	goto, realend
536	endelse
537
538	nofile:
539	print, 'www : no data to read'
540	goto, realend
541
542	realend:
543	PRINT, 'www : fin programme extract_amsua ',SYSTIME()
544	end

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

Download in other formats:

Original Format