Changeset 486 for trunk

Timestamp:

01/20/12 18:50:41 (12 years ago)

Author:

lelod

Message:

amelioration de l execution

File:

• trunk/src/extract_amsua.pro (modified) (11 diffs)

trunk/src/extract_amsua.pro

@@ -246 +246 @@
 tbmin=100
 tbmax=350
-
+PRINT, 'debut programme',SYSTIME() 
 ; 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
+domdef,lon_min-15, lon_max+15, lat_min-15, lat_max+15
 bate = ncdf_lec(file, var = 'z') GT 0
 xxe = float(ncdf_lec(file, var = 'lon'))
@@ -255 +255 @@
 jpie = n_elements(xxe)
 jpje = n_elements(yye)
-
+PRINT, 'lecture bathy terminee',SYSTIME() 
 ; appel au ssprgm qui interpole les fonctions de correction (inutile
 ; de l'appeler a chaque fichier....)
@@ -266 +266 @@
 
 ; lecture de la liste des fichiers AMSU-A
-list_filea = project_id_env+'list_file'
+list_file = project_id_env+'list_file'
 ;list_filea='list_file'
     ;index_filea= 0
-nb_filea = file_lines(list_filea)
-a = STRARR(nb_filea)
+nb_file = file_lines(list_file)
+a = STRARR(nb_file)
    ; filea = STRARR(nb_filea)
-print, 'ouverture pour lecture de ', list_filea
-openr, lun_a, list_filea, /get_lun
+print, 'ouverture pour lecture de ', list_file
+openr, lun_a, list_file, /get_lun
 filea=''
 nlist=0
 ilist=0
-;loadct,39
-;tek_color
-numwindow=2
-;window,numwindow
-;plot,indgen(30)+1,cor_l+200,psym=2,color=0,yrange=[150,250]
-;oplot,indgen(30)+1,cor_s+200,psym=2,color=100
-
-;while ilist le nlist do begin
+PRINT,'demarrage boucle sur lichiers', SYSTIME() 
 while (not eof(lun_a)) do begin
    ilist=ilist+1
@@ -297 +290 @@
    COMMON amsua_header,ama_head
    COMMON amsua_data  ,ama_scan
-   print, 'ouverture et lecture du fichier ', filea
+   print, 'ouverture et lecture du fichier ', filea, SYSTIME() 
    openr,lu1,filea,Error=erra,/get_lun
    read_amsua1c,filea, flag1
+   close,lu1
+   free_lun,lu1
    if (flag1 eq 0) then goto, labfile
    na=SIZE(reform(ama_scan.btemps[0,*,*]))
@@ -322 +317 @@
 ; latitude, et les Tbs des differents canaux (15 pour AMSUA)
 
-; 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
+; extraction du canal traite dans un domaine enveloppant la zone choisie
+   midpix=fix(nbpix/2)
+   amalong=REFORM(ama_scan.latlon[1,*,*]/1.E4)
+   amalati=REFORM(ama_scan.latlon[0,*,*]/1.E4)
+   amcha=REFORM(ama_scan.btemps[nocanal,*,*]/100.)
+   amzeni=REFORM(ama_scan.angles[0,*,*]/100.)
+   ttt=REFORM(ama_scan.scnlintime/3600000.)
+
+   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)
+   if nzon ne 0 then begin
+      amalat=amalati[*,jnd]
+      amalon=amalong[*,jnd]
+      amch=amcha[*,jnd]
+      amzen=amzeni[*,jnd]
+      tt=ttt[jnd]
+      dims = SIZE(amalon, /DIMENSIONS)
+      print,'dimension tableaux extraits',dims
+      n_scan=dims[1]
+      PRINT,'il y a des donnees dans la region ', SYSTIME() 
+      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
-
+      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)
-   landseamask=intarr(nbpix,n_scan)
-   for isc=0L,n_scan-1L do begin
-      for ifo=0,nbpix-1 do begin
+         ch_nadir=fltarr(nbpix,n_scan)
+      ; ch_nadir=amch[ifo,isc]
+         landseamask=intarr(nbpix,n_scan)+2 ; valeur hors zone selectionnee
+         PRINT,'boucle sur les points du fichier, correction nadir', SYSTIME() 
+
+         for isc=0L,n_scan-1L do begin
+;recherche de la zone xxe,yye englobant la fauchee
+            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[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
-   endfor
-
-moych=fltarr(nbpix)
-
-for ifo=0,nbpix-1 do begin
-   jnd=where(ch_nadir[ifo,*] gt tbmin and ch_nadir[ifo,*] lt tbmax and amalon[ifo,*] ge lon_min and amalon[ifo,*] le lon_max $
-         and amalat[ifo,*] ge lat_min and amalat[ifo,*] le lat_max,npts)
-  if npts ne 0 then begin
-     moych[ifo]= mean(ch_nadir[ifo,jnd])
-  endif
-endfor
+               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[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
+         endfor
+         PRINT, 'fin boucle',SYSTIME() 
+         moych=fltarr(nbpix)
+      
+;for ifo=0,nbpix-1 do begin
+;   jnd=where(ch_nadir[ifo,*] gt tbmin and ch_nadir[ifo,*] lt tbmax and amalon[ifo,*] ge lon_min and amalon[ifo,*] le lon_max $
+;         and amalat[ifo,*] ge lat_min and amalat[ifo,*] le lat_max,npts)
+;  if npts ne 0 then begin
+;     moych[ifo]= mean(ch_nadir[ifo,jnd])
+;  endif
+;endfor
 ;print,'correction',moych
 ;print,'cor_s',cor_s
@@ -381 +392 @@
 ;print,'ecarts max a la valeur moyenne au nadir',min(moych)-mean(moych[12:17]),max(moych)-mean(moych[12:17])
 
-
-numwindow=numwindow+1
 
 ; appel a interpolswath pour ajuster les pixels amsua sur une grille
@@ -390 +399 @@
 ; pbs avec les fauchees qui traversent le meridien 180deg!!!!
 ; l'interpolation ajoute des points parasites dans la zone...
-   nbgrid=0
-cont=0L
-   for i=0L,n_scan-1L do begin
-      tb=ch_nadir[*,i]
-      lon=amalon[*,i]
-      lat=amalat[*,i]
-      mask=landseamask[*,i]
-      ind=where(tb gt tbmin and lon gt lon_min-15 and lon lt lon_max+15,nbon) 
+         nbgrid=0
+         cont=0L
+         chint=fltarr(1)
+         for i=0L,n_scan-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
-         if cont eq 1L then begin
-            fovgrid=indgen(nbgrid)+1
-         endif
+            if nbon eq nbpix then begin
+               cont=cont+1L
+               nbgrid1=nbgrid
+               interpolswath,tb,lat,lon,mask,resol,nbgrid,tbgrid,latgrid,longrid,maskgrid
+               if cont eq 1L then begin
+                  fovgrid=indgen(nbgrid)+1
+               endif
          ;longit=longrid
          ;latit=latgrid
@@ -411 +421 @@
 
 ; 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) and (tbgrid gt tbmin) and (tbgrid lt tbmax), npt)
-         if npt ne 0 then begin
-            if (n_elements(chint) EQ 0) 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
-         endif else begin
+               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)
+               if npt ne 0 then begin
+                  if (n_elements(chint) EQ 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
+               endif else begin
             ;print, 'www : aucun point dans la zone'
-         endelse
-      endif
-
-   endfor
-
-
-   nn=n_elements(chint)
-
-   print,'nb points dans la zone',nn
-
-   for i=0L, nn-1L do begin
-      if (n_elements(data) EQ 0) then begin
-         data = { desc : desc $
+               endelse
+            endif
+         endfor
+     ; PRINT,'fin interpolation et selection des points de la zone', SYSTIME() 
+
+         nn=n_elements(chint)
+
+         print,'nb points dans la zone',nn
+
+         for i=0L, nn-1L do begin
+            if (n_elements(data) EQ 0) then begin
+               data = { desc : desc $
                   , hour : timeint[i] $
                   , fov : fovint[i] $
@@ -451 +460 @@
                   , tb : chint[i] $
                 }
-      endif else begin
-         data = [data , { desc : desc $
+            endif else begin
+               data = [data , { desc : desc $
                           , hour : timeint[i] $
                           , fov : fovint[i] $
@@ -460 +469 @@
                           , tb : chint[i] $
                         }]
-      endelse
-   endfor
+            endelse
+         endfor
+ 
 ; fin boucle sur les fichiers lus
+   endif
+  ; PRINT,'passage au fichier suivant', SYSTIME() 
 endwhile
 labfile:
@@ -488 +500 @@
   ;  print, 'www : no data to write'
 endelse
-; loadct,39
 ; window,0
-; plot,lonint,latint,xrange=[lon_min,lon_max],yrange=[lat_min,lat_max],xstyle=1,ystyle=1,/nodata
+loadct,39
+ plot,lonint,latint,xrange=[lon_min,lon_max],yrange=[lat_min,lat_max],xstyle=1,ystyle=1,/nodata
 ; ind=where(maskint eq 1)
 ; ;oplot,lonint[ind],latint[ind],psym=1
-; jnd=where(chint gt tbmin and chint lt tbmax)
+;jnd=where(chint gt tbmin and chint lt tbmax)
 ; tb=chint[jnd]
 ; lolo=lonint[jnd]
 ; lala=latint[jnd]
 ; colscale=intarr(n_elements(tb))
-; for j=0,n_elements(tb)-1 do begin
-;    colscale[j]=(tb[j]-min(tb))*255/(max(tb)-min(tb))
-; endfor
+; for j=0,n_elements(tb)-1 do    colscale[j]=(tb[j]-min(tb))*255/(max(tb)-min(tb))
+ ;endfor
 ; print,min(colscale),max(colscale)
 ;    plots,lolo,lala,psym=5,symsize=1,color=colscale
+tb=chint
+ lolo=lonint
+ lala=latint
+ colscale=intarr(n_elements(tb))
+ for j=0,n_elements(tb)-1 do    colscale[j]=(tb[j]-min(tb))*255/(max(tb)-min(tb))
+ plots,lolo,lala,psym=5,symsize=1,color=colscale,/noerase
+
 end