Changeset 541

Timestamp:

04/25/12 17:50:48 (12 years ago)

Author:

lelod

Message:

evolution ponctuelle

File:

• trunk/src/extract_amsu.pro (modified) (5 diffs)

trunk/src/extract_amsu.pro

@@ -342 +342 @@
       COMMON amsub_header,amb_head
       COMMON amsub_data  ,amb_scan
-      print, 'ouverture et lecture du fichier ', filea, SYSTIME()
-      openr,lu1,filea,Error=erra,/get_lun
+      print, 'ouverture et lecture du fichier ', fileb, SYSTIME()
+      openr,lu1,fileb,Error=errb,/get_lun
       read_amsub1c,fileb, flag2
       close,lu1
@@ -350 +350 @@
       nb=SIZE(reform(amb_scan.btemps[0,*,*]))
       nb_scan=nb[2]             ; nb lignes (ou fauchees)
-     n_scan=na_scan
+      n_scan=nb_scan
       nb_fov=nb[1]              ; nb pixels dans la fauchee
       nbpix = nb_fov          ; nb pixels dans la fauchee AMSUB
       nosat=amb_head.h_satid  ; on le lit dans le fichier, donc pas besoin de l'info dans la liste
-
       amalong=REFORM(amb_scan.latlon[1,*,*]/1.E4)
       amalati=REFORM(amb_scan.latlon[0,*,*]/1.E4)
@@ -361 +360 @@
       ttt=REFORM(amb_scan.scnlintime/3600000.)
       amaday=REFORM(amb_scan.scnlindy)
- 
-
-
-            
-      fov=indgen(nbpix) ; nb pixels dans la fauchee AMSU
-      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
-  midpix=fix(nbpix/2)
+   endif
+
+   fov=indgen(nbpix)            ; nb pixels dans la fauchee AMSU
+   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
+   midpix=fix(nbpix/2)
    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)
    print,"iii : nb de points du fichier dans le domaine geographique +15deg ",nzon, SYSTIME()
@@ -389 +386 @@
 ;
 ; correction nadir des donnees
-         ch_nadir=fltarr(nbpix,nzon)
-         landseamask=intarr(nbpix,nzon)+2 ; valeur hors zone selectionnee
-        ; PRINT,'boucle sur les points du fichier, correction nadir', SYSTIME()
-
-         for isc=0L,nzon-1L do begin
+      ch_nadir=fltarr(nbpix,nzon)
+      landseamask=intarr(nbpix,nzon)+2 ; valeur hors zone selectionnee
+; PRINT,'boucle sur les points du fichier, correction nadir', SYSTIME()
+      for isc=0L,nzon-1L do begin
 ;recherche de la zone xxe,yye englobant la fauchee
-            utilex=where(xxe ge amalon[midpix,isc]-15 and xxe le amalon[midpix,isc]+15,nxu)
-            if nxu ne 0 then begin
-               xxutile=xxe[utilex]
-               batutilx=bate[utilex,*]
-               utiley=where(yye ge amalat[midpix,isc]-8 and yye le amalat[midpix,isc]+8,nyu)
-               if nyu ne 0 then begin
-                  yyutile=yye[utiley]
-                  batutil=batutilx[*,utiley]
-
-                  for ifo=0,nbpix-1 do begin
+         utilex=where(xxe ge amalon[midpix,isc]-15 and xxe le amalon[midpix,isc]+15,nxu)
+         if nxu ne 0 then begin
+            xxutile=xxe[utilex]
+            batutilx=bate[utilex,*]
+            utiley=where(yye ge amalat[midpix,isc]-8 and yye le amalat[midpix,isc]+8,nyu)
+            if nyu ne 0 then begin
+               yyutile=yye[utiley]
+               batutil=batutilx[*,utiley]
+
+               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 (xxutile ge amalon[ifo,isc]-swath[ifo]/400. and xxutile le amalon[ifo,isc]+swath[ifo]/400.,nxlandsea)
-                     yind=where (yyutile ge amalat[ifo,isc]-track[ifo]/400. and yyutile le amalat[ifo,isc]+track[ifo]/400.,nylandsea)
-                     if(nxlandsea ne 0 and nylandsea ne 0) then begin
-                        bate1=batutil[xind,*]
-                        bate2=bate1[*,yind]
-                        if (mean(bate2) le 0.5) then begin
-                           ch_nadir[ifo,isc]=amch[ifo,isc]-cor_s[ifo]
-                           landseamask[ifo,isc]=0 ; mer
-                        endif else begin
-                           ch_nadir[ifo,isc]=amch[ifo,isc]-cor_l[ifo]
-                           landseamask[ifo,isc]=1 ; terre
-                        endelse
-                     endif
-                  endfor
-               endif
+                  xind=where (xxutile ge amalon[ifo,isc]-swath[ifo]/400. and xxutile le amalon[ifo,isc]+swath[ifo]/400.,nxlandsea)
+                  yind=where (yyutile ge amalat[ifo,isc]-track[ifo]/400. and yyutile le amalat[ifo,isc]+track[ifo]/400.,nylandsea)
+                  if(nxlandsea ne 0 and nylandsea ne 0) then begin
+                     bate1=batutil[xind,*]
+                     bate2=bate1[*,yind]
+                     if (mean(bate2) le 0.5) then begin
+                        ch_nadir[ifo,isc]=amch[ifo,isc]-cor_s[ifo]
+                        landseamask[ifo,isc]=0 ; mer
+                     endif else begin
+                        ch_nadir[ifo,isc]=amch[ifo,isc]-cor_l[ifo]
+                        landseamask[ifo,isc]=1 ; terre
+                     endelse
+                  endif
+               endfor
             endif
-
-         endfor
+         endif
+
+      endfor
         ; PRINT, 'fin boucle ',SYSTIME()
-         moych=fltarr(nbpix)
+      moych=fltarr(nbpix)
 
 ; appel a interpolswath pour ajuster les pixels amsua sur une grille
@@ -436 +432 @@
 ; pbs avec les fauchees qui traversent le meridien 180deg!!!!
 ; l'interpolation ajoute des points parasites dans la zone...
-         nbgrid=0
-         cont=0L
-         chint=fltarr(1)
-         for i=0L,nzon-1L do begin
-            tb=ch_nadir[*,i]
-            lon=amalon[*,i]
-            lat=amalat[*,i]
-            mask=landseamask[*,i]
-            ind=where(tb gt tbmin,nbon)
+      nbgrid=0
+      cont=0L
+      chint=fltarr(1)
+      for i=0L,nzon-1L do begin
+         tb=ch_nadir[*,i]
+         lon=amalon[*,i]
+         lat=amalat[*,i]
+         mask=landseamask[*,i]
+         ind=where(tb gt tbmin,nbon)
 ; test pour eviter pb d'interpolation de longitude (meridien 180)
-            if nbon eq nbpix then begin
-               cont=cont+1L
-               nbgrid1=nbgrid
-               interpolswath,tb,lat,lon,mask,resol,nbgrid,tbgrid,latgrid,longrid,maskgrid                  
-               fovgrid=indgen(n_elements(tbgrid))+1
+         if nbon eq nbpix then begin
+            cont=cont+1L
+            nbgrid1=nbgrid
+            interpolswath,tb,lat,lon,mask,resol,nbgrid,tbgrid,latgrid,longrid,maskgrid                  
+            fovgrid=indgen(n_elements(tbgrid))+1
 ;print,'wwwtest: nb points fov interpole ',n_elements(tbgrid)
 ; selection des taches au sol situees dans la zone d'interet
-               zone=where((longrid ge lon_min) and (longrid le lon_max) $
+            zone=where((longrid ge lon_min) and (longrid le lon_max) $
                           and (latgrid ge lat_min) and (latgrid le lat_max) and (tbgrid gt tbmin) and (tbgrid lt tbmax), npt)
 ;print,'wwwtest: nb points zone d interet ',npt
-               if npt ne 0 then begin
-                  if (n_elements(chint) LE 1) then begin
-                     timeint=replicate(tt[i],npt)
-                     chint=tbgrid[zone]
-                     latint=latgrid[zone]
-                     lonint=longrid[zone]
-                     fovint=fovgrid[zone]
-                     maskint=maskgrid[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,replicate(tt[i],npt)]
-                  endelse
+            if npt ne 0 then begin
+               if (n_elements(chint) LE 1) then begin
+                  timeint=replicate(tt[i],npt)
+                  chint=tbgrid[zone]
+                  latint=latgrid[zone]
+                  lonint=longrid[zone]
+                  fovint=fovgrid[zone]
+                  maskint=maskgrid[zone]
                endif else begin
-            ;print, 'www : aucun point dans la zone'
+                  chint=[chint,tbgrid[zone]]
+                  latint=[latint,latgrid[zone]]
+                  lonint=[lonint,longrid[zone]]
+                  fovint=[fovint,fovgrid[zone]]
+                  maskint=[maskint,maskgrid[zone]]
+                  timeint=[timeint,replicate(tt[i],npt)]
                endelse
-            endif
+            endif else begin
+                                ;print, 'www : aucun point dans la zone'
+            endelse
+         endif
+      endfor
+ ; PRINT,'fin interpolation et selection des points de la zone ', SYSTIME()
+
+      nn=n_elements(chint)
+
+      print,'www : nb points dans la zone en fin de traitement de l''orbite ',nn
+      if nn gt 1 then begin
+         for i=0L, nn-1L do begin
+            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:  maskint[i] $
+                                , tb : chint[i] $
+                              }]
+            endelse
          endfor
-     ; PRINT,'fin interpolation et selection des points de la zone ', SYSTIME()
-
-         nn=n_elements(chint)
-
-         print,'www : nb points dans la zone en fin de traitement de l''orbite ',nn
-         if nn gt 1 then begin
-            for i=0L, nn-1L do begin
-               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:  maskint[i] $
-                                   , tb : chint[i] $
-                                 }]
-               endelse
-            endfor
-            desc=desc+1; incrementation du numero de l'orbite
-         endif
+         desc=desc+1            ; incrementation du numero de l'orbite
+      endif
 ;
 ; fin boucle sur les fichiers lus