Changeset 558

Timestamp:

05/22/12 09:59:38 (12 years ago)

Author:

llmlod

Message:

grille reguliere fic plusieurs orbites

File:

• trunk/src/amsu2ncdf_bis.pro (modified) (6 diffs)

trunk/src/amsu2ncdf_bis.pro

@@ -12 +12 @@
 ; +todo+
 ;
-;    :param use_amsua: flag to tell if amsub only or amsub
-;    :type use_amsua: int 0 or 1
-;    :raise use_amsua: required
+;    :param numch: number of the channel
+;    :type numch: string
+;    :raise numch: required
 ;
 ;    :param yyyy: year
@@ -85 +85 @@
 ; To produce AMSU netCDF file for 20060801::
 ;
-;   use_amsua=1
+;   numch='a5'
 ;   yyyy=2006L
 ;   mm=8
@@ -97 +97 @@
 ; or::
 ;
-;   use_amsua=0
+;   numch='a5'
 ;   yyyy=2006L
 ;   mm=8
@@ -131 +131 @@
 ;
 ;-
-PRO amsu2ncdf, use_amsua, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max,resol
+PRO amsu2ncdf, numch, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max, resol
 ;
 compile_opt idl2, strictarrsubs
 ;
 @cm_project
+
+
 testfilename=''
-prefix=numch
 look = 'filename'
 
 geomin = geolocation_to_string_idl(lon_min, lat_min, look,1)
 geomax = geolocation_to_string_idl(lon_max, lat_max, look,1)
+
+testfilename='/homedata/eymard/varamma_d/AMSU/' $
+   + string(yyyy,format='(I4.4)') + '/' $
+   + string(mm,format='(I2.2)')   + '/' $
+   + numch + '_' + string(yyyy,format='(I4.4)') $
+   + string(mm,format='(I2.2)') $
+   + string(dd,format='(I2.2)') + '_' $
+   + geomin + '_' + geomax + '.dat'
 
 result=file_amsu_t2_to_mem( yyyy,mm,dd,numch,lon_min,lon_max,lat_min,lat_max,testfilename)
@@ -156 +165 @@
 ; grille de sortie xaxis, yaxis :
 ; --------------------------------
-;    lon : (indgen((maxlon-minlon)/resolt+1)+long(minlon/resolt))*resolt+resolt/2., $
-;    lat : (indgen((maxlat-minlat)/resolt+1)+long(minlat/resolt))*resolt+resolt/2. $
-xaxis=indgen((lmon_max-lon_min)/resol+1)
-yaxis=indgen((lat_max-lat_min)/resol+1)
-print,xaxis,yaxis
+xaxis=indgen((lon_max-lon_min)/resol+1)*resol+lon_min
+yaxis=indgen((lat_max-lat_min)/resol+1)*resol+lat_min
 nlon=n_elements(xaxis)
 nlat=n_elements(yaxis)
-;minlon=min(xaxis)
-;maxlon=max(xaxis)
-;minlat=min(yaxis)
-;maxlat=max(yaxis)
-;resolt=abs(yaxis[1]-yaxis[0])
-
-norbit=max(desc)
-mask=fltarr(nlon,nlat)
-tb=fltarr(nlon,nlat,nbo)
-maskg[*]=!values.f_nan
-tbg[*]=!values.f_nan
+
+; calcul seuil pour remplir la nouvelle grille
+; --------------------------------------------
+; Ã  adapter pour AMSU-B
+; pour AMU-A aussi car ne rempli pas tout les points de grille avec
+; ces seuils
+if strmid(numch,0,1) eq 'a' then begin
+; AMSU-A
+    ind_fov19=WHERE(fov eq 19)
+    ind_fov20=WHERE(fov eq 20)
+    
+    seuil_lon=ABS(lon[ind_fov19[1]]-lon[ind_fov20[1]])/2.
+    seuil_lat=ABS(lat[ind_fov19[1]]-lat[ind_fov19[2]])/2.
+; arrondi au 0.01 superieur pour lisser un peu mais pas trop
+    seuil_lon=CEIL(seuil_lon*100.)/100.
+    seuil_lat=CEIL(seuil_lat*100.)/100.
+
+;seuil_lon=.3
+;seuil_lat=.25
+    
+    print, seuil_lon
+    print, seuil_lat
+
+endif else begin
+; AMSU-B
+
+endelse
+
+; numero d'orbite
+norb=desc[UNIQ(desc)]
+; nombre d'orbites
+nbo=n_elements(norb)
+
+; initialisation
+maskg=fltarr(nlon,nlat)
+tbg=fltarr(nlon,nlat,nbo)
+temp=fltarr(nlon,nlat,nbo)
 iheur=fltarr(nbo)
-temp=fltarr(nlon,nlat)
-
-norbit=0
-for i=0,n_elements(tb1)-1 do begin
-  while desc[i] eq norbit do begin
-     ind=where(abs(lon[i]-xaxis) lt resol/2 and (abs(lat[i]-yaxis) lt resol/2,nbpt)
-     if nbpt ne 0 then begin
-        tbg[ind,norbit]=tb1[i]
-        maskg[ind]=mask[i]
-        temp[ind]=hour[i]
-     endif 
-  endwhile
-  iheur[norbit]=mean(temp)
-  norbit=norbit+1
+
+; valeur erreur -999.
+val_err=-999.
+maskg[*]=val_err
+tbg[*]=val_err
+temp[*]=val_err
+
+;;; pour chaque orbite ;;;
+for k=0,nbo-1 do begin
+    ; trouve tous les indices de l'orbite
+    ind_orb=WHERE(desc eq norb[k])
+;;; pour chaque latitude ;;;
+    for j=0,nlat-1 do begin
+        ; trouve la latitude d'origine la plus proche de la nouvelle latitude
+        ; en fonction du seuil choisi
+        ind_lat=WHERE(ABS(yaxis[j]-lat[ind_orb]) lt seuil_lat)
+        ; teste si un unique indice est trouvé
+        if (n_elements(ind_lat) eq 1) then begin
+        ; teste si l'indice trouvé n'est pas vide
+            if (ind_lat ne -1) then begin
+;;; pour chaque longitude ;;;
+                for i=0,nlon-1 do begin
+                    ind_lon=WHERE(ABS(xaxis[i]-lon[ind_orb[ind_lat]]) lt seuil_lon)
+                    ; si plusieurs indices sont trouvés, calcule la moyenne
+                    if (n_elements(ind_lon) gt 1) then begin
+                        tbg[i,j,k]=MEAN(tb1[ind_orb[ind_lat[ind_lon]]])
+                        temp[i,j,k]=MEAN(hour[ind_orb[ind_lat[ind_lon]]])
+                    endif else begin
+                        if (ind_lon ne -1) then begin
+                            tbg[i,j,k]=tb1[ind_orb[ind_lat[ind_lon]]]
+                            temp[i,j,k]=hour[ind_orb[ind_lat[ind_lon]]]
+                            maskg[i,j]=mask[ind_orb[ind_lat[ind_lon]]]
+                        endif
+                    endelse
+                endfor
+            endif
+        endif else begin
+;;; pour chaque longitude ;;;
+            for i=0,nlon-1 do begin
+                ind_lon=WHERE(ABS(xaxis[i]-lon[ind_orb[ind_lat]]) lt seuil_lon)
+                if (n_elements(ind_lon) gt 1) then begin
+                    tbg[i,j,k]=MEAN(tb1[ind_orb[ind_lat[ind_lon]]])
+                    temp[i,j,k]=MEAN(hour[ind_orb[ind_lat[ind_lon]]])
+                endif else begin
+                    if (ind_lon ne -1) then begin
+                        tbg[i,j,k]=tb1[ind_orb[ind_lat[ind_lon]]]
+                        temp[i,j,k]=hour[ind_orb[ind_lat[ind_lon]]]
+                        maskg[i,j]=mask[ind_orb[ind_lat[ind_lon]]]
+                    endif
+                endelse
+            endfor
+        endelse
+    endfor
+tps=temp[*,*,k]
+ind = WHERE(tps eq val_err)
+tps[ind]=!VALUES.F_NAN
+iheur[k]=MEAN(tps,/NaN)
 endfor
+
+; decomposition de l'heure de l'orbite en heure, minute, seconde
+h_orb=fix(iheur)
+m_orb=fix((iheur-h_orb)*60)
+s_orb=fix((iheur-h_orb-m_orb/60.)*3600)
 
   ;
   ; Creation du fichier NetCDF :
   ; ----------------------------
-  htime=fix(iheur)
-  mtime=fix((iheur-htime)*60)
-  stime=fix((iheur-htime-mtime/60.)*3600)
-  caldat, julday(1,1,yyyy)+iday-1,mttime,dtime,ytime
-  file_ncdf=project_id_env+ 'AMSU/' $
-        + string(yyyy,format='(I4.4)') + '/' $
-        + string(mm,format='(I2.2)')   + '/' $
-        + prefix+'ch5_' $
-        + string(yyyy,format='(I4.4)') $
-        + string(mm,format='(I2.2)') $
-        + string(dd,format='(I2.2)') $
-        +'_' $
-        + string(htime,format='(I2.2)') $
-        + string(mtime,format='(I2.2)') $
-        + string(stime,format='(I2.2)') + '_' $
-        + geomin + '_' $
-        + geomax + '_' $
-        + 'nadir.nc'
-  print, 'iii : avant ', 'create_amsu_netcdf pour ' + file_ncdf
-  create_amsu_netcdf, file_ncdf, ytime, mttime, dtime, htime, mtime, stime $
-      ,xaxis,yaxis, tbg, maskg
+
+file_ncdf=project_id_env+ 'AMSU/' $
+  + string(yyyy,format='(I4.4)') + '/' $
+  + string(mm,format='(I2.2)')   + '/' $
+  + numch + '_' $
+  + string(yyyy,format='(I4.4)') $
+  + string(mm,format='(I2.2)') $
+  + string(dd,format='(I2.2)') + '_' $
+  + geomin + '_' $
+  + geomax + '_' $
+  + 'nadir.nc'
+
+create_amsu_netcdf, file_ncdf, yyyy, mm, dd, h_orb, m_orb, s_orb $
+  ,xaxis,yaxis, tbg, maskg, val_err, numch  
+
   ; comment the next line to produce an image
-  goto, ret0
+;  goto, ret0
   ;
   ; Visualisation :
   ; ---------------
-  window
+minlon=min(xaxis)
+maxlon=max(xaxis)
+minlat=min(yaxis)
+maxlat=max(yaxis)
+  window, 0
   loadct, 39
   map_set, limit=[minlat,minlon,maxlat,maxlon]
 
-  tabcolor = MAP_IMAGE(bch5,Startx,Starty, COMPRESS=1, $
+  tabcolor = MAP_IMAGE(tbg[*,*,1],Startx,Starty, COMPRESS=1, $
     LATMIN=minlat, LONMIN=minlon, $
     LATMAX=maxlat, LONMAX=maxlon, $
-    min_value=0, missing=-999.)
-
-  lewh=where(tabcolor ne -999., cnt, complement=lewhcmpl)
+    max_value=280., min_value=241., missing=val_err)
+
+  lewh=where(tabcolor ne val_err, cnt, complement=lewhcmpl)
   if (cnt ne n_elements(tabcolor)) then begin
     tabcolor[lewhcmpl]=255
@@ -246 +324 @@
   MAP_CONTINENTS, /coasts, color=0
 
+;saveimage, '20060802_ficdat_orb2_sajust.png', PNG=PNG
+
+   window, 1
+   loadct, 39
+   map_set, limit=[minlat,minlon,maxlat,maxlon]
+
+   tabcolor = MAP_IMAGE(maskg,Startx,Starty, COMPRESS=1, $
+     LATMIN=minlat, LONMIN=minlon, $
+     LATMAX=maxlat, LONMAX=maxlon, $
+     max_value=2., min_value=-1., missing=val_err)
+
+   lewh=where(tabcolor ne val_err, cnt, complement=lewhcmpl)
+   if (cnt ne n_elements(tabcolor)) then begin
+     tabcolor[lewhcmpl]=255
+     tabcolor[lewh]=bytscl(tabcolor[lewh], top=254)
+   endif else begin
+     tabcolor=bytscl(tabcolor, top=254)
+   endelse
+
+;   ; Display the warped image on the map at the proper position:
+   TV, tabcolor, Startx, Starty
+
+;   ; Draw gridlines over the map and image:
+   MAP_GRID, latdel=10, londel=10, /LABEL, /HORIZON  , color=0
+
+;   ; Draw continent outlines:
+   MAP_CONTINENTS, /coasts, color=0
+
+;saveimage, '20060802_mask_orb2_sajust.png', PNG=PNG
+
   ret0:
-endwhile
-if (debug eq 1) then begin
-   print, 'fin grand while'
-endif
-
-close, lun
-free_lun, lun
 
 end
+