Changeset 597

Timestamp:

07/23/12 18:07:26 (12 years ago)

Author:

llmlod

Message:

creation fichier nc sur grille reguliere MAJ

Location:

trunk/src

Files:

• amsu2ncdf.pro (modified) (6 diffs)
• create_amsu_netcdf.pro (modified) (5 diffs)

trunk/src/amsu2ncdf.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
@@ -28 +28 @@
 ;    :raise dd: required
 ;
-;    :param lonmin: longitude min, W < 0
-;    :units lonmin: deg
-;    :type lonmin: double
-;    :raise lonmin: required
-;
-;    :param lonmax: longitude max, W < 0
-;    :units lonmax: deg
-;    :type lonmax: double
-;    :raise lonmax: required
-;
-;    :param latmin: latitude min, N > 0
-;    :units latmin: deg
-;    :type latmin: double
-;    :raise latmin: required
-;
-;    :param latmax: latitude max, N > 0
-;    :units latmax: deg
-;    :type latmax: double
-;    :raise latmax: required
+;    :param lon_min: longitude min, W < 0
+;    :units lon_min: deg
+;    :type lon_min: double
+;    :raise lon_min: required
+;
+;    :param lon_max: longitude max, W < 0
+;    :units lon_max: deg
+;    :type lon_max: double
+;    :raise lon_max: required
+;
+;    :param lat_min: latitude min, N > 0
+;    :units lat_min: deg
+;    :type lat_min: double
+;    :raise lat_min: required
+;
+;    :param lat_max: latitude max, N > 0
+;    :units lat_max: deg
+;    :type lat_max: double
+;    :raise lat_max: required
 ;
 ; read :file:`${PROJECT_ID}/AMSU/yyyy/mm/dd/amsua[_amsub]_{yyyymmdd}_{geomin}_{geomax}_nadir.dat`
@@ -85 +85 @@
 ; To produce AMSU netCDF file for 20060801::
 ;
-;   use_amsua=1
+;   numch='a5'
 ;   yyyy=2006L
 ;   mm=8
 ;   dd=01
-;   lonmin=-60.
-;   lonmax=50.
-;   latmin=-30.
-;   latmax=45.
-;   amsu2ncdf, use_amsua, yyyy, mm, dd, lonmin, lonmax, latmin, latmax
+;   lon_min=-60.
+;   lon_max=50.
+;   lat_min=-30.
+;   lat_max=45.
+;   amsu2ncdf, use_amsua, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max
 ;
 ; or::
 ;
-;   use_amsua=0
+;   numch='a5'
 ;   yyyy=2006L
 ;   mm=8
 ;   dd=01
-;   lonmin=-60.
-;   lonmax=50.
-;   latmin=-30.
-;   latmax=45.
-;   amsu2ncdf, use_amsua, yyyy, mm, dd, lonmin, lonmax, latmin, latmax
+;   lon_min=-60.
+;   lon_max=50.
+;   lat_min=-30.
+;   lat_max=45.
+;   amsu2ncdf, use_amsua, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max
 ;
 ;
@@ -118 +118 @@
 ; Called by :ref:`traite_amsuab.sh`
 ;
-; Previous step : :ref:`correct_nadir_amsu.pro`
+; Previous step : :ref:`extract_amsu.pro`
 ;
 ; Next step : ... science !
@@ -125 +125 @@
 ; ====
 ;
-; mettre a jour nouvelle structure et nouveau fichier de référence
-;
-; coding rules
-;
-; pourquoi autant de calcul de tmpmin
-;
-; accept uncorrected amsu file
-;
-; ajout parametre pour changer de grille ou de fichier (ie Karim a fourni à Laurence un
-; nouveau fichier MSG sur une nouvelle grille).
-;
-; created daily file  : need to change loop and create_amsu_netcdf.pro !!
 ;
 ; EVOLUTIONS
 ; ==========
 ;
-; $Id: amsu2ncdf.pro 336 2011-08-10 12:07:19Z pinsard $
-;
-; $URL$
-;
-; - fplod 20110818T155126Z aedon.locean-ipsl.upmc.fr (Darwin)
-;
-;   * fix array syntax (not detected on idl8 climserv !)
-;
-; - pinsard 2011-08-18T09:40:48Z loholt1.ipsl.polytechnique.fr (Linux)
-;
-;   * avoid amsub__ch5 file name
-;
-; - pinsard 2011-08-17T14:58:13Z loholt1.ipsl.polytechnique.fr (Linux)
-;
-;   * use geolocation_to_string_idl
-;
-; - pinsard 2011-08-10T14:11:58Z loholt1.ipsl.polytechnique.fr (Linux)
-;
-;   * complete description with parameters
-;   * add lonmin, lonmax, latmin, latmax parameter
-;   * align parameters order to correct_nadir_amsu
-;   * change output filename
-;
-; - lelod 2011-08-09
-;
-;   * mise a jour lecture des fichiers issus de
-;     extract_amsuab.pro et correct_nadir_amsu.pro
-;   * modification des formats et adaptation aux cas use_amsua=0 ou 1
-;   * remplacement des -999. par !values.f_nan
-;   * /MSG/2006/08/01/200608010000_msg-tb108_map_15min.nc have been replaced by
-;     project_id_env + '/MSG/2006/06/01/msg_IR108_20060601.nc'
-;
-; - pinsard 2011-08-10T11:58:43Z loholt1.ipsl.polytechnique.fr (Linux)
-;
-;   * remove blank in input filename
-;   * complete description
-;
-; - pinsard 2011-05-30T15:15:52Z loholt1.ipsl.polytechnique.fr (Linux)"
-;
-;   * remove noaa and satellite number from filename to avoid
-;     difference between alphanumeric and chronologic order
-;
-; - fplod 20110530T125608Z aedon.locean-ipsl.upmc.fr (Darwin)
-;
-;   * new parameters for :ref:`create_amsu_netcdf.pro`
-;   * limit stdout if debug=1 (more if debug=2)
-;
-; - pinsard 2011-05-27T08:10:27Z loholt1.ipsl.polytechnique.fr
-;
-;   * replace dirout by proetc_id_env
-;   * add compile_opt
-;
-; - pinsard 2011-05-26T15:17:36Z loholt1.ipsl.polytechnique.fr
-;
-;   * add dd, use_amsua, use_amsub parameters
-;
-; - fplod 20110422T160635Z aedon.locean-ipsl.upmc.fr (Darwin)
-;
-;   * use PROJECT
-;
-; - fplod 20110121T150550Z aedon.locean-ipsl.upmc.fr (Darwin)
-;
-;   * externalize function and pro
-;
-; - pinsard 2010-07-23T10:42:49Z loholt1.ipsl.polytechnique.fr (Linux)
-;
-;   * add graphic
-;   * file save in ${PROJECT_OD}/AMSU/yyyy/mm/
-;   * remove condition yrday eq 241 and hour gt 21
 ;
 ;-
-PRO amsu2ncdf, use_amsua, yyyy, mm, dd, lonmin, lonmax, latmin, latmax
+PRO amsu2ncdf, numch, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max, resol
 ;
 compile_opt idl2, strictarrsubs
 ;
 @cm_project
-;
-;set debug to 1 to print info
-debug=0
-;
-;
-s = size(SCOPE_TRACEBACK(/STRUCTURE),/DIMENSION)
-routine = (SCOPE_TRACEBACK(/STRUCTURE))[s-1].Routine
-if keyword_set(key_performance) EQ 1 THEN BEGIN
-    time1 = SYSTIME(1)
-ENDIF
-; Initialisation des variables :
-; ------------------------------
-if (use_amsua eq 1) then begin
-   prefix='amsua_amsub_'
-endif else begin
-   prefix='amsub_'
-endelse
-
-look='filename'
-scale = 1.
-geomin=geolocation_to_string_idl(lonmin,latmin, look, scale)
-geomax=geolocation_to_string_idl(lonmax,latmax, look, scale)
-filename=project_id_env+ 'AMSU/' $
-        + string(yyyy,format='(I4.4)') + '/' $
-        + string(mm,format='(I2.2)') + '/' $
-        + prefix+string(yyyy,format='(I4.4)') $
-        + string(mm,format='(I2.2)') $
-        + string(dd,format='(I2.2)') + '_' $
-        + geomin + '_' $
-        + geomax + '_' $
-        + 'nadir.dat'
-;
-; Lecture de la grille réguliere :
+
+
+testfilename=''
+look = 'filename'
+
+; ne pas utiliser les bornes lon et lat fournies en entree car elles
+; peuvent etre differentes de celle du fichier .dat puisqu'on peut
+; creer un fichier .nc extrait du .dat
+;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)') + '_' $
+   + '*.dat'
+;   + geomin + '_' + geomax + '.dat'
+
+print, testfilename
+
+result=file_amsu_t2_to_mem( yyyy,mm,dd,numch,lon_min,lon_max,lat_min,lat_max,testfilename)
+         desc=result.data.desc
+         hour=result.data.hour
+         fov=result.data.fov
+         lon=result.data.lon
+         lat=result.data.lat
+         mask=result.data.landseamask
+         tb1=result.data.tb
+         ; decodage du nb de points (nn)
+         nn=n_elements(tb1)
+
+; numero d'orbite
+norb=desc[UNIQ(desc)]
+; nombre d'orbites
+nbo=n_elements(norb)
+
+; grille de sortie xaxis, yaxis :
 ; --------------------------------
-;remplacement des fichiers msg initiaux (2006 seulement) par une
-;nopuvelle version, avec grille reellement reguliere, crees par
-;K. Ramage l'ete 2011. Ici un exemple, copie dans varamma_d
-;verifier la compatibilite avec les outils appeles ensuite
-;domaine -25 0 W, 0 - 20 N
-; old grid filemsgref=project_id_env + '/MSG/2006/08/01/200608010000_msg-tb108_map_15min.nc'
-filemsgref=project_id_env + '/EPSAT/extracted-surf-rr_epsat-sg_multi-sat_010d_30min_20060801_v3.1-02c.nc'
-print, 'Lecture de la grille reguliere ',filemsgref
-reggrid=read_regular_grid(filemsgref)
-minlon=min(reggrid.lon)
-maxlon=max(reggrid.lon)
-minlat=min(reggrid.lat)
-maxlat=max(reggrid.lat)
-resolt=abs(reggrid.lon[1]-reggrid.lon[0])
-;resolt=0.1
-;reggrid={ $
-;    lon : (indgen((maxlon-minlon)/resolt+1)+long(minlon/resolt))*resolt+resolt/2., $
-;    lat : (indgen((maxlat-minlat)/resolt+1)+long(minlat/resolt))*resolt+resolt/2. $
-;}
-if (debug eq 1) then begin
-    print,'bords du domaine, lon, lat',minlon,maxlon,minlat,maxlat
-endif
-;
-; Déclaration du tableau de reprojection :
-; ----------------------------------------
-;
-; Lecture des données AMSU :
-; --------------------------
-print, 'Lecture du fichier de donnees ',filename
-openr, lun, filename, /GET_LUN, ERROR = error
-IF (error NE 0) then begin
-   ras = report(['eee : can not open for reading'$
-         + '!C' $
-         + 'code : ' + !ERROR_STATE.MSG $
-         + filename])
-   STOP
-ENDIF
-
-;
-; Boucle sur les lignes du fichier :
-; ----------------------------------
-print, '--> Saut des donnees deja traitees'
-
-kont1=0L
-repeat begin
-   if (use_amsua eq 1) then begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                     ach1_adj,ach2_adj,ach3,ach4_adj,ach5_adj,ach6_adj,ach7_adj,$
-                     ach8_adj,ach15_adj,bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,20(1x,f8.3))'
-   endif else begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                    bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,8(1x,f8.3))'
-   endelse
-
-  kont1=kont1+1L
-  if (debug eq 1) then begin
-      condition1=(lon ge minlon)
-      condition2=(lon le maxlon)
-      condition3=(lat ge minlat)
-      condition4=(lat le maxlat)
-      condition5=(yrday eq 241)
-      condition6=(hour gt 21)
-      condition7=eof(lun)
-      condition=(lon ge minlon and $
-                 lon le maxlon and $
-                 lat ge minlat and $
-                 lat le maxlat and yrday eq 241 and hour gt 21) or eof(lun)
-      print,'lecture 1 lu',kont1, lon,lat,yrday,hour,condition
-      print,'lecture 1',condition,condition1,condition2,condition3,condition4,condition5,condition6,condition7
-      print,'lecture 1 seuil repeat', minlon,maxlon,minlat,maxlat,'241','21'
-  endif
-
-endrep until ((lon ge minlon and $
-                 lon le maxlon and $
-                 lat ge minlat and $
-                 lat le maxlat) or eof(lun) )
-
-if (debug eq 1) then begin
-    print,'fin lecture 1', kont1, lon,lat, eof(lun)
-endif
-while (not eof (lun)) do begin
-
-  if (debug eq 1) then begin
-     print,'entree grand while'
-  endif
-  mask=fltarr(n_elements(reggrid.lon), n_elements(reggrid.lat))
-  bch5=fltarr(n_elements(reggrid.lon), n_elements(reggrid.lat))
-  mask[*]=!values.f_nan
-  bch5[*]=!values.f_nan
-  print, '--> Recherche des premiers pixels dans la fenetre'
-  kont2=0L
-  if (debug eq 1) then begin
-     print, 'dernier lu ',lon,lat,yrday,hour
-  endif
-  repeat begin
-   if (use_amsua eq 1) then begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                     ach1_adj,ach2_adj,ach3,ach4_adj,ach5_adj,ach6_adj,ach7_adj,$
-                     ach8_adj,ach15_adj,bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,20(1x,f8.3))'
-   endif else begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                    bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,8(1x,f8.3))'
-   endelse
-    if (debug eq 1) then begin
-        kont2=kont2+1L
-        print,'lecture 2 lu',kont2, lon,lat,yrday,hour
-        print,'lecture 2 seuil repeat', minlon,maxlon,minlat,maxlat
-    endif
-  endrep until ((lon ge minlon and $
-                 lon le maxlon and $
-                 lat ge minlat and $
-                 lat le maxlat) or eof(lun) )
-
-  if (debug eq 1) then begin
-      print,'fin lecture 2', kont2, lon,lat
-  endif
-  print, "--> Lecture de donnees de l'orbite"
-  iheur = hour
-  iday = yrday
-  numsat = nosat
-  print,'orbite',numsat,iday,iheur
-  lontab=[!values.f_nan]
-  lattab=[!values.f_nan]
-  masktab=[!values.f_nan]
-  bch5tab=[!values.f_nan]
-  kont3=0L
-  while (hour eq iheur and not eof(lun)) do begin
-    lontab=[lontab,lon]
-    lattab=[lattab,lat]
-    masktab=[masktab,land_sea]
-    bch5tab=[bch5tab,bch5_adj]
-   if (use_amsua eq 1) then begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                     ach1_adj,ach2_adj,ach3,ach4_adj,ach5_adj,ach6_adj,ach7_adj,$
-                     ach8_adj,ach15_adj,bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,20(1x,f8.3))'
-   endif else begin
-
-       readf,lun, month,land_sea,yrday,hour,fov,nosat,lon,lat,zen,$
-                    bch1_adj,bch2_adj,bch3_adj,bch4_adj,bch5_adj,$
-                    format='(3(i3,1x),f8.4,2i3,8(1x,f8.3))'
-   endelse
-    kont3=kont3+1
-    if (debug gt 1) then begin
-        print,'lecture 3 lu',kont2, lon,lat,yrday,hour
-    endif
-  endwhile
-    if (debug gt 1) then begin
-       print, 'tableau donnee', bch5tab
-    endif
-  ;
-  ; Extraction des données AMSU dans la fenetre :
-  ; ---------------------------------------------
-  lewh=where(lontab ge minlon and $
-           lontab le maxlon and $
-           lattab ge minlat and $
-           lattab le maxlat and $
-           bch5tab ne -999., cnt)
-  if (debug eq 1) then begin
-      print,'debut extraction ',cnt
-  endif
-  if (cnt gt 1) then begin
-    lontab=reform(lontab[lewh])
-    lattab=reform(lattab[lewh])
-    masktab=reform(masktab[lewh])
-    bch5tab=reform(bch5tab[lewh])
-  endif else begin
-    goto, ret0
-  endelse
-  ;
-  ; Reprojection dans une grille reguliere :
-  ; ----------------------------------------
-  print, "Reprojection des donnees dans grille reguliere"
-  print, "--> Recherche de la fenetre encadrant l'orbite"
-; ++ pourquoi 5 calculs de tmpmin
-  tmpmin = min(abs(reggrid.lon-min(lontab)), imin)
-  tmpmin = min(abs(reggrid.lon-max(lontab)), imax)
-  tmpmin = min(abs(reggrid.lat-min(lattab)), jmin)
-  tmpmin = min(abs(reggrid.lat-max(lattab)), jmax)
-
-  tmpmin = min((reggrid.lon[imin]-lontab)^2+(reggrid.lat[jmin]-lattab)^2,iwh)
-  print, "--> Recherche pixel par pixel du plus proche voisin",imin,imax,jmin,jmax
-  for i=imin,imax do begin
-    for j=jmin,jmax do begin
-; ++ encore un calcul de tmpmin !!
-      tmpmin = min((reggrid.lon[i]-lontab)^2+(reggrid.lat[j]-lattab)^2,iwh)
-      tmptab=(lontab-lontab[iwh])^2+(lattab-lattab[iwh])^2
-      minvrf = min(tmptab[where(tmptab ne 0)])
-      if (tmpmin le minvrf) then begin
-        mask[i,j] = masktab[iwh]
-        bch5[i,j] = bch5tab[iwh]
-        if (debug gt 1) then begin
-            print, ' remplissage bch5',i,j
-        endif
-     endif
-    endfor
-  endfor
-  status=where(bch5 ne !values.f_nan,nb_bch5_novm)
-  if (debug eq 1) then begin
-      print, 'compteur de bch5 non vm',nb_bch5_novm
-  endif
-  if (nb_bch5_novm eq 0) then begin
-     print, 'eee : aucun point non vm'
-     print, 'eee : we should stop here'
-  endif
+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)
+
+; calcul seuil pour remplir la nouvelle grille
+; --------------------------------------------
+;seuil latitude = MAX((variation de latitude sur un meme fov)/2)
+;pour une orbite quelconque (on choisie la plus longue)
+;moyenne pour fov sur le bord (=min(fov)) et au milieu (=fix(max(fov)/2))
+torb=LONARR(nbo)
+FOR i=0,nbo-1 DO BEGIN
+    inutile=WHERE(desc EQ norb[i],nb)
+    torb[i]=nb
+ENDFOR
+longorb=norb[WHERE(torb EQ MAX(torb))]
+
+ind_bor=WHERE(fov EQ MIN(fov) AND desc EQ longorb[0])
+dif_bor=ABS(lat[ind_bor[1:N_ELEMENTS(ind_bor)-1]]-lat[ind_bor[0:N_ELEMENTS(ind_bor)-2]])
+; arrondi au 0.1 superieur pour lisser un peu mais pas trop
+seuil_lat=CEIL(MAX([dif_bor])/2*10)/10.
+;print, seuil_lat
+
+;seuil longitude = MAX((difference de longitude sur une ligne)/2)
+;pour une orbite quelconque (on choisit la plus longue)
+;pour une ligne quelconque (on choisit la plus longue)
+ind_orb=WHERE(desc EQ longorb[0])
+dif_lon=lon[ind_orb[1:N_ELEMENTS(ind_orb)-1]]-lon[ind_orb[0:N_ELEMENTS(ind_orb)-2]]
+; arrondi au 0.1 superieur pour lisser un peu mais pas trop
+seuil_lon=CEIL(MAX(dif_lon[WHERE(dif_lon LT 1)])/2*10)/10.
+;print, seuil_lon
+
+; initialisation
+maskg=fltarr(nlon,nlat)
+tbg=fltarr(nlon,nlat,nbo)
+temp=fltarr(nlon,nlat,nbo)
+iheur=fltarr(nbo)
+
+; valeur erreur -999.
+val_err=-999.
+maskg[*]=val_err
+tbg[*]=val_err
+temp[*]=val_err
+
+;;; pour chaque orbite ;;;
+rT=6371 ; rayon de la Terre
+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
+;;; pour chaque longitude ;;;
+        for i=0,nlon-1 do begin
+            ind_pt=WHERE(ABS(yaxis[j]-lat[ind_orb]) LT seuil_lat AND ABS(xaxis[i]-lon[ind_orb]) LT seuil_lon)
+            IF ind_pt[0] NE -1 THEN BEGIN
+;;;cherhce le plus proche en km;;;
+; lon/lat en radian
+                lat_rad=lat[ind_orb[ind_pt]]*!Pi/180.
+                lon_rad=lon[ind_orb[ind_pt]]*!Pi/180.
+                y_rad=yaxis[j]*!Pi/180.
+                x_rad=xaxis[i]*!Pi/180.
+; calcul distance
+                dist = 2*rT*ASIN(SQRT((SIN((y_rad-lat_rad)/2.)*SIN((y_rad-lat_rad)/2.))+(COS(x_rad)*COS(lon_rad)*SIN((x_rad-lon_rad)/2.)*SIN((x_rad-lon_rad)/2.))))
+; selectionne le plus proche
+                dist_min=WHERE(dist EQ MIN(dist))
+; attribut le point le plus proche
+                tbg[i,j,k]=tb1[ind_orb[ind_pt[dist_min]]]
+                temp[i,j,k]=hour[ind_orb[ind_pt[dist_min]]] 
+                maskg[i,j]=mask[ind_orb[ind_pt[dist_min]]]
+            ENDIF 
+        ENDFOR
+    ENDFOR
+    tps=temp[*,*,k]
+    ind = WHERE(tps eq val_err)
+    tps[ind]=!VALUES.F_NAN
+    iheur[k]=MEAN(tps,/NaN)
+;print, 'orbite '+string(norb[k])+' finie'
+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 $
-      , reggrid.lon, reggrid.lat, bch5, mask
+
+geomin = geolocation_to_string_idl(lon_min, lat_min, look,1)
+geomax = geolocation_to_string_idl(lon_max, lat_max, look,1)
+
+ 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
@@ -520 +316 @@
   MAP_CONTINENTS, /coasts, color=0
 
+saveimage, 'test.png', PNG=PNG
+;saveimage, '20060802_ficdat_orb2_sajust.png', PNG=PNG
+
+   window, 3
+   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, 'essai_mask.png', PNG=PNG
+;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
-
-IF key_performance EQ 1 THEN BEGIN
-    msg = report(['ppp : ' + routine + ' : durée totale ' $
-                   + string(SYSTIME(1)-time1,format='(F12.6)')])
-ENDIF
 
 end
+

trunk/src/create_amsu_netcdf.pro

@@ -128 +128 @@
 ; ------------------
 time_id = NCDF_VARDEF(id, 'time', [dtime_id], /double)
+NCDF_ATTPUT, id, time_id, 'long_name', 'Time'
 NCDF_ATTPUT, id, time_id, 'units', $
     STRING(format = '("hours since ",I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2," UTC")', $
     ytime_org, mtime_org, dtime_org, htime_org, mmtime_org, stime_org)
 NCDF_ATTPUT, id, time_id, 'calendar', 'standard'
-NCDF_ATTPUT, id, time_id, 'long_name', 'Time'
 NCDF_ATTPUT, id, time_id, 'actual_range', [min(timetab),max(timetab)], /double
 ;
@@ -138 +138 @@
 ; ---------------------
 latd_id = NCDF_VARDEF(id, 'lat', [dlatd_id], /float)
+NCDF_ATTPUT, id, latd_id, 'long_name', 'Latitude'
 NCDF_ATTPUT, id, latd_id, 'units', 'degrees_north'
-NCDF_ATTPUT, id, latd_id, 'long_name', 'Latitude'
 NCDF_ATTPUT, id, latd_id, 'actual_range', [min(lats),max(lats)]
 ;
@@ -145 +145 @@
 ; ----------------------
 long_id = NCDF_VARDEF(id, 'lon', [dlong_id], /float)
+NCDF_ATTPUT, id, long_id, 'long_name', 'Longitude'
 NCDF_ATTPUT, id, long_id, 'units', 'degrees_east'
-NCDF_ATTPUT, id, long_id, 'long_name', 'Longitude'
 NCDF_ATTPUT, id, long_id, 'actual_range', [min(lons),max(lons)], /float
 ;
@@ -152 +152 @@
 ; -----------------
 data_id = NCDF_VARDEF(id, 'tb'+numch, [dlong_id,dlatd_id,dtime_id], /float)
+NCDF_ATTPUT, id, data_id, 'long_name', "AMSU-"+strupcase(strmid(numch,0,1))+" Channel "+strmid(numch,1,1)+" Brightness Temperature"
 NCDF_ATTPUT, id, data_id, 'units', "Kelvin"
-NCDF_ATTPUT, id, data_id, 'long_name', "AMSU-"+strupcase(strmid(numch,0,1))+" Channel "+strmid(numch,1,1)+" Brightness Temperature"
-NCDF_ATTPUT, id, data_id, 'valid_range', [100,500]
+;NCDF_ATTPUT, id, data_id, 'valid_range', [100,500]
 if (debug eq 1) then begin
    print, 'dans create_amsu_netcdf'
@@ -167 +167 @@
 ; Variable : Mask :
 ; -----------------
-mask_id = NCDF_VARDEF(id, 'mask', [dlong_id,dlatd_id,dtime_id], /float)
-NCDF_ATTPUT, id, mask_id, 'units', '1'
-NCDF_ATTPUT, id, mask_id, 'long_name', "land-sea mask, 1 all land, 0 all sea"
+mask_id = NCDF_VARDEF(id, 'mask', [dlong_id,dlatd_id], /float)
+;mask_id = NCDF_VARDEF(id, 'mask', [dlong_id,dlatd_id,dtime_id], /float)
+NCDF_ATTPUT, id, mask_id, 'units', '1 = land, 0 = sea'
+NCDF_ATTPUT, id, mask_id, 'long_name', "land-sea mask"
 NCDF_ATTPUT, id, mask_id, 'valid_range', [0,1]
 NCDF_ATTPUT, id, mask_id, 'missing_value', val_err