Changeset 457 for trunk

Timestamp:

12/12/11 23:10:30 (12 years ago)

Author:

lelod

Message:

integration correction nadir dans extract_amsua

Location:

trunk/src

Files:

• extract_amsua.pro (modified) (12 diffs)
• interpol_correc.pro (modified) (2 diffs)
• interpolswath.pro (modified) (3 diffs)

trunk/src/extract_amsua.pro

@@ -8 +8 @@
 ; ==================
 ;
-; prgm de lecture des fichiers AMSU (A seulement)
+; prgm de lecture des fichiers AMSU (A seulement, a valider avant de
+; generaliser pour amsub aussi)
 ;
 ; decode les noms des fichiers dans la date choisie, puis
 ; appelle le prgm de lecture
+;
+; appelle interpol_correc, qui fournit les fonctions de correction au
+; nadir et les applique, en fonction du masque terre-mer issu de la
+; carte de bathymetrie de S. Masson
+;
+; applique une interpolation tenant compte de la dimension de la tache
+; au sol selon la position dans la fauchee
 ;
 ; ecrit les donnees pour la zone (lat/long) choisie
@@ -50 +58 @@
 ;   yyyy
 ;   mm
+;   resol
 ;   lon_min
 ;   lon_max
@@ -69 +78 @@
 ;   IDL> mm=8
 ;   IDL> dd=1
+;   IDL> resol=1
 ;   IDL> lon_min=-60.
 ;   IDL> lon_max=50.
@@ -79 +89 @@
 ; TODO
 ;
+; lelod 20111212
+;
+;
+; :History:
+;
+; $Id$
+;
 ; lelod 20111211
-; attention, il faut transformer les donnees 2D en 1D
-; reprendre methode "extract_amsuab", qui utilisait un "where"?
-; ou passer directement dans l'ecriture???
-;
-; :History:
-;
-; $Id$
-;
-; lelod 20111211
-; mis en commentaires la partie interpolation
+;  reintegre les modules de correct_nadir_amsu dans extract
 ; prevu decodage de numch pour connaitre le nom instrument et le
-;numero du canal (pas fini, au debut)
+; numero du canal (pas fini, au debut) --> a faire aussi dans les
+; ssprgm appeles (interpol_correc)
+; espere que c'est le dernier avatar de cette !!!!?XX^! de
+; chaine de traitement AMSU
+;
 ; $URL: svn+ssh://lelod@forge.ipsl.jussieu.fr/ipsl/forge/projets/varamma/svn/trunk/src/extract_amsuab.pro $
 ;
@@ -188 +200 @@
 ;
 ;-
-PRO extract_amsua, numch, yyyy, mm, dd,lon_min, lon_max, lat_min, lat_max
+PRO extract_amsua, numch, yyyy, mm, dd, resol, lon_min, lon_max, lat_min, lat_max
 ;
 compile_opt idl2, strictarrsubs
@@ -199 +211 @@
 
 
+tbmin=100
+tbmax=350
+
+
+; lecture fichier land - sea (S. Masson)
+file=project_id_env+'/MASK/ETOPO1_Ice_g_gmt4.nc'
+domdef,lon_min, lon_max, lat_min, lat_max
+
+bate = ncdf_lec(file, var = 'z') GT 0
+xxe = float(ncdf_lec(file, var = 'lon'))
+yye = float(ncdf_lec(file, var = 'lat'))
+jpie = n_elements(xxe)
+jpje = n_elements(yye)
+
+; extraction zone donnees AMSU 
+domdef,lon_min,lon_max,lat_min,lat_max
+bate = ncdf_lec(file, var = 'z') GT 0
+xxe = float(ncdf_lec(file, var = 'lon'))
+yye = float(ncdf_lec(file, var = 'lat'))
+jpie = n_elements;module lecture fichier a introduire ici
+
+; inutile maintenant: ouverture du fichier de sortie
 look = 'filename'
 scale = 1.
@@ -213 +247 @@
         + geomax $
         + '.dat'
-
 print, 'ouverture pour ecriture de ', fileout
 openw, lun_obs, fileout, /get_lun
+
+
+; appel au ssprgm qui interpole les fonctions de correction (inutile
+; de l'appeler a chaque fichier....)
+if nomcanal eq 'a' then nbpixel=30
+if nomcanal eq 'b' then nbpixel=90
+interpol_correc,nomcanal,nbpixel,cor_l,cor_s,swath,track
+
 
 ; ouverture des fichiers liste (annee, mois, jour, tous satellites) pour
@@ -223 +264 @@
 ; lecture de la liste des fichiers AMSU-A
     ;list_filea = project_id_env+'list_filea'
-     list_filea='list_filea'
+list_filea='list_filea'
     ;index_filea= 0
-    nb_filea = file_lines(list_filea)
-    a = STRARR(nb_filea)
+nb_filea = file_lines(list_filea)
+a = STRARR(nb_filea)
    ; filea = STRARR(nb_filea)
-    print, 'ouverture pour lecture de ', list_filea
-    openr, lun_a, list_filea, /get_lun
-    filea=''
-    nlist=1
-    ilist=1
+print, 'ouverture pour lecture de ', list_filea
+openr, lun_a, list_filea, /get_lun
+filea=''
+nlist=1
+ilist=1
    ; while (not eof(lun_a)) do begin
-    while ilist le nlist do begin
-      ilist=ilist+1
-        onefile = ''
-        readf, lun_a, filea
+while ilist le nlist do begin
+   ilist=ilist+1
+   onefile = ''
+   readf, lun_a, filea
        ; filea[index_filea] = onefile
         ; isolate string independant from satellite
@@ -244 +285 @@
     ;endwhile
 
-        COMMON amsua_header,ama_head
-        COMMON amsua_data  ,ama_scan
-        print, 'ouverture et lecture du fichier ', filea
-        openr,lu1,filea,Error=erra,/get_lun
-        read_amsua1c,filea, flag1
-        if (flag1 eq 0) then goto, labfile
-        na=SIZE(reform(ama_scan.btemps[0,*,*]))
-        na_scan=na[2]           ; nb lignes (fauchees) dans l'orbite
-        na_fov=na[1]            ; nb points dans la fauchee
-        nbpix= na_fov
-        nosat=ama_head.h_satid
-        ; nb pixels dans la fauchee AMSUA
-        fov=indgen(nbpix)
-        n=na
-        n_scan=na_scan
-        n_fov=na_fov
-
-    ntime_1=fltarr(nbpix,n_scan)
-    fovy_1=fltarr(nbpix,n_scan)
-
-    for j=0L, nbpix-1L do begin
-        fovy_1[j,*]=j+1 ; on remplit un tableau de nscan lignes avec les valeurs de no de pixel
-    endfor
+   COMMON amsua_header,ama_head
+   COMMON amsua_data  ,ama_scan
+   print, 'ouverture et lecture du fichier ', filea
+   openr,lu1,filea,Error=erra,/get_lun
+   read_amsua1c,filea, flag1
+   if (flag1 eq 0) then goto, labfile
+   na=SIZE(reform(ama_scan.btemps[0,*,*]))
+   na_scan=na[2]                ; nb lignes (fauchees) dans l'orbite
+   na_fov=na[1]                 ; nb points dans la fauchee
+   nbpix= na_fov
+   nosat=ama_head.h_satid
+                                ; nb pixels dans la fauchee AMSUA
+   fov=indgen(nbpix)
+   n=na
+   n_scan=na_scan
+   n_fov=na_fov
+
+   ntime_1=fltarr(nbpix,n_scan)
+   fovy_1=fltarr(nbpix,n_scan)
+
+   for j=0L, nbpix-1L do begin
+      fovy_1[j,*]=j+1           ; on remplit un tableau de nscan lignes avec les valeurs de no de pixel
+   endfor
 ;
 ; on lit la structure et on extrait les infos: temps, longitude,
@@ -273 +314 @@
 ; extraction du canal traite
 
-    amch=REFORM(ama_scan.btemps[nocanal,*,*]/100.)
-    amzen=REFORM(ama_scan.angles[0,*,*]/100.)
-    amalat=REFORM(ama_scan.latlon[0,*,*]/1.E4)
-    amalon=REFORM(ama_scan.latlon[1,*,*]/1.E4)
-    tt=REFORM(ama_scan.scnlintime/3600000.)
-    amaday=REFORM(ama_scan.scnlindy)
-    amafov= fovy_1
-    jour_ama=amaday[0]          ; jour du passage
+   amch=REFORM(ama_scan.btemps[nocanal,*,*]/100.)
+   amzen=REFORM(ama_scan.angles[0,*,*]/100.)
+   amalat=REFORM(ama_scan.latlon[0,*,*]/1.E4)
+   amalon=REFORM(ama_scan.latlon[1,*,*]/1.E4)
+   tt=REFORM(ama_scan.scnlintime/3600000.)
+   amaday=REFORM(ama_scan.scnlindy)
+   amafov= fovy_1
+   jour_ama=amaday[0]           ; jour du passage
 ; recherche du sens de l'orbite : descendante (1) ou montante (0)
-    lat0=amalat[*,0]
-    lat1=amalat[*,1]
-    if (mean(lat1-lat0) lt 0) then begin
-        desc=1
-    endif else begin
-        desc=0
-    endelse
-; attention tableaux a deux dimensions
-
-
-
-; appel a interpolswatw pour ajuster les pixels amsua sur une grille
+   lat0=amalat[*,0]
+   lat1=amalat[*,1]
+   if (mean(lat1-lat0) lt 0) then begin
+      desc=1
+   endif else begin
+      desc=0
+   endelse
+
+; correction nadir des donnees
+   ch_nadir=fltarr(nbpix,n_scan)
+   lanseamask=intarr(nbpix,n_scan)
+   for isc=0L,n_scan-1L do begin
+      for ifo=0,nbpix-1 do begin
+; approximation: 1km= 1deg/100, et on considere le quart en lon et lat
+; de la surface du pixel pour optimiser la localisation du centre du
+; pixel attention approximation brutale des distances en fractions dedegre
+; recherche de la cote par rapport a la resolution AMSUA
+         xind=where (xxe ge amalon[ifo,isc]-swath[ifo]/400. and xxe le amalon[ifo,isc]+swath[ifo]/400.,nxlandsea)
+         yind=where (yye ge amalat[ifo,isc]-track[ifo]/400. and yye le amalat[ifo,isc]+track[ifo]/400.,nylandsea)
+         if(nxlandsea ne 0 and nylandsea ne 0) then begin
+            bate1=bate[xind,*]
+            bate2=bate1[*,yind]
+            if (mean(bate2) le 0.5) then begin
+               ch_nadir[ifo,isc]=amch[ifo,isc]-cor_s[numcanal,ifo]
+               landseamask[ifo,isc]=0 ; mer
+            endif else begin
+               ch_nadir[ifo,isc]=amch[ifo,isc]-cor_l[numcanal,ifo]
+               landseamask[ifo,isc]=1 ; terre
+            endelse
+         endif
+      endfor
+   endfor
+ 
+
+; appel a interpolswath pour ajuster les pixels amsua sur une grille
 ; reguliere dans la fauchee
 ; et selection de la zone conservee
 ;tek_color
 ;resol=1
-;nbgrid=0
-;chint=fltarr(100)
-;latint=fltarr(100)
-;lonint=fltarr(100)
-;fovint=fltarr(100)
-;timeint=fltarr(100)
-;time=fltarr(100)
-;n_scan-1L
-;for i=0L,10L do begin
-;   tb=amch[*,i]
-;   lon=amalon[*,i]
-;   lat=amalat[*,i]
-  ; nbgrid1=nbgrid
-;   interpolswath,tb,lat,lon,resol,nbgrid,tbgrid,latgrid,longrid
-;     if i eq 0L then begin
+   nbgrid=0
+   for i=0L,n_scan-1L do begin
+      tb=ch_nadir[*,i]
+      lon=amalon[*,i]
+      lat=amalat[*,i]
+      mask=landseamask[*,i]
+      nbgrid1=nbgrid
+      interpolswath,tb,lat,lon,mask,resol,nbgrid,tbgrid,latgrid,longrid,maskgrid
+      if i eq 0L then begin
        ; plot,lon,tb,psym=1
        ; oplot,longrid,tbgrid,psym=4;
-;        fovgrid=indgen(nbgrid)+1
-;     endif
+         fovgrid=indgen(nbgrid)+1
+      endif
   ; if nbgrid ne nbgrid1 then begin
   ;    oplot,lon,tb,psym=1,color=3
    ;   oplot,longrid,tbgrid,psym=4,color=5
    ;   endif
-;  zone=where((longrid ge lon_min) and (longrid le lon_max) $
-;             and (latgrid ge lat_min) and (latgrid le lat_max),npt)
+; selection des taches au sol situees dans la zone d'interet
+      zone=where((longrid ge lon_min) and (longrid le lon_max) $
+                 and (latgrid ge lat_min) and (latgrid le lat_max),npt)
 ;print, 'npt', npt
-;  if npt ne 0 then begin
-;     time[zone]=tt[i]
-;     chint=[chint,tbgrid[zone]]
-;     latint=[latint,latgrid[zone]]
-;     lonint=[lonint,longrid[zone]]
-;     fovint=[fovint,fovgrid[zone]]
-;     timeint=[timeint,time[zone]]
-;  endif else begin
-;     print, 'www : aucun point dans la zone'
-;  endelse
-
-; endfor
-
-
-
-      ; replace missing values -9999.99 by NaN
-  index_missing = where(chint EQ -9999.99, nb_missing)
-  if (nb_missing NE 0) THEN BEGIN
-     tb[index_missing]=!values.f_nan
-  endif
-; param landseamask mis fictivement pour unification du format
-        landseamask=!values.f_nan
-
-print,'nb points dans la zone',nn
-
-      for i=0L, nn-1L do begin
+      if npt ne 0 then begin
+         time[zone]=tt[i]
+ ; alerte: comment initialiser les tableaux qu'on concatene????
+         if (n_elements(chint) EQ 0) then begin
+            chint=tbgrid[zone]
+            latint=latgrid[zone]
+            lonint=longrid[zone]
+            fovint=fovgrid[zone]
+            maskint=maskgrid[zone]
+            timeint= time[zone]
+         endif else begin
+            chint=[chint,tbgrid[zone]]
+            latint=[latint,latgrid[zone]]
+            lonint=[lonint,longrid[zone]]
+            fovint=[fovint,fovgrid[zone]]
+            maskint=[maskint,maskgrid[zone]]
+            timeint=[timeint,time[zone]]
+         endelse
+      endif else begin
+         print, 'www : aucun point dans la zone'
+      endelse
+      
+   endfor
+
+
+   for i=0L, npt-1L do begin
+      if (( chint[i] lt tbmin) or (chint[i] gt tbmax) ) then begin
+         ch_adj[i]=!values.f_nan
+      endif
+   endfor
+
+
+   print,'nb points dans la zone',npt
+
+   for i=0L, npt-1L do begin
 ;format: nosat,mois, day, lon, lat, zen,fov, time, ama_ch1:15,amb_ch1:5
 ; print, 'mois::::', mm,'-----', nab-i
-         printf,lun_obs, mm, jour_ama, desc, nosat, fovint[i], landseamask, lonint[i], latint[i], timeint[i],$
-                chint[i],format='(6(i3,1x),4(f8.4,1x))'
-         if (n_elements(data) EQ 0) then begin
-             data = { desc : desc $
-                    , hour : timeint[i] $
-                    , fov : fovint[i] $
-                    , lon : lonint[i] $
-                    , lat: latint[i] $
-                    , landseamask: 9$
-                    , tb : chint[i] $
-                    }
-         endif else begin
-             data = [data , { desc : desc $
+      printf,lun_obs, mm, jour_ama, desc, nosat, fovint[i], landseamask, lonint[i], latint[i], timeint[i],$
+             chint[i],format='(6(i3,1x),4(f8.4,1x))'
+      if (n_elements(data) EQ 0) then begin
+         data = { desc : desc $
+                  , hour : timeint[i] $
+                  , fov : fovint[i] $
+                  , lon : lonint[i] $
+                  , lat: latint[i] $
+                  , landseamask: maskint[i] $
+                  , tb : chint[i] $
+                }
+      endif else begin
+         data = [data , { desc : desc $
                           , hour : timeint[i] $
                           , fov : fovint[i] $
                           , lon : lonint[i] $
                           , lat: latint[i] $
-                          , landseamask: 9$
+                          , landseamask:  maskint[i] $
                           , tb : chint[i] $
-
-
+                        }]
+      endelse
+   endfor
 ; fin boucle sur les fichiers lus
-labfile:
 endwhile
 print, 'fermeture de ', fileout
@@ -381 +449 @@
 if (n_elements(data) NE 0) then begin
     ; write outputfile
+   ; ajouter resol dans le header
     header = { yyyy : yyyy $
              ,  mm : mm $
@@ -396 +465 @@
     print, 'www : no data to write'
 endelse
-;
+
 end

trunk/src/interpol_correc.pro

@@ -107 +107 @@
 ;
 ;-
-PRO interpol_correc,canal,nbpix,cor_l,cor_s
+PRO interpol_correc,canal,nbpix,cor_l,cor_s,swath,track
 ;
 @cm_project
@@ -115 +115 @@
 
 pixelsize,pixatot,pixbtot,alongatot,alongbtot
+if nomcanal eq 'a' then begin
+   swath=pixatot
+   track=alongatot
+endif else begin
+   swath=pixbtot
+   track=alongbtot
+endelse
 
 ; parametres AMSU et calcul de la fauchee

trunk/src/interpolswath.pro

@@ -41 +41 @@
 ;
 ;-
-pro interpolswath, tb, latlu, lonlu, resol, nbgrid, tbint, latgrid, longrid
+pro interpolswath, tb, latlu, lonlu, masklu, resol, nbgrid, tbint, latgrid, longrid, mask
 
 pixelsize,pixatot,pixbtot,alongatot,alongbtot
@@ -92 +92 @@
 ;print,tb
 tbint=fltarr(nbgrid)
+mask=intarr(nbgrid)
 for i=0,na-1 do begin
    ind=where(abs(grid-fova[i]) le pixatot[i]/2,nii)
@@ -97 +98 @@
    if nii ne 0 then begin
       tbint[ind]=tb[i]
+      mask[ind]=masklu[i]
    endif
    if nii eq 0 then begin