source: trunk/src/amsu2ncdf.pro @ 599

; docformat = 'rst'
;+
; .. _amsu2ncdf.pro:
;
; =============
; amsu2ncdf.pro
; =============
;
; DESCRIPTION
; ===========
;
; +todo+
;
;    :param numch: number of the channel
;    :type numch: string
;    :raise numch: required
;
;    :param yyyy: year
;    :type yyyy: int
;    :raise yyyy: required
;
;    :param mm: month
;    :type mm: int
;    :raise mm: required
;
;    :param dd: day
;    :type dd: int
;    :raise dd: 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`
; produced by :ref:`correct_nadir_amsu.pro`.
;
; read :file:`${PROJECT_ID}/MSG/2006/06/01/msg_IR108_20060601.nc` that will
; be used for regrid.++
;
; write :file:`${PROJECT_ID}/AMSU/yyyy/mm/dd/amsua[_amsub]_{yyyymmdd}_{hhmmss}_{geomin}_{geomax}_nadir.nc`.
;
; .. only:: man
;
;    Figure is visible on PDF and HTML documents only.
;
; .. only:: html or latex
;
;     .. graphviz::
;
;        digraph amsu2ncdf {
;
;           amsu_t3 [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/yyyy/mm/amsua_amsub_yyyymmdd_geomin_geomax_nadir.dat"];
;           amsu_t4 [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/yyyy/mm/amsua_amsub_yyyymmdd_geomin_geomax_nadir.nc"];
;           msg [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/MSG/2006/06/01/msg_IR108_20060601.nc"];
;
;           amsu2ncdf [shape=box,
;           fontname=Courier,
;           color=blue,
;           URL="http://forge.ipsl.jussieu.fr/varamma/browser/trunk/src/amsu2ncdf.pro",
;           label="${PROJECT}/src/amsu2ncdf.pro"];
;
;           {amsu_t3 msg} -> {amsu2ncdf} -> {amsu_t4}
;
;          }
;
; EXAMPLES
; ========
;
; To produce AMSU netCDF file for 20060801::
;
;   numch='a5'
;   yyyy=2006L
;   mm=8
;   dd=01
;   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::
;
;   numch='a5'
;   yyyy=2006L
;   mm=8
;   dd=01
;   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
;
;
; SEE ALSO
; ========
;
; :ref:`data_amsu`
; :ref:`data_msg`
;
; Use: :func:`read_regular_grid`, :ref:`create_amsu_netcdf.pro`, :func:`geolocation_to_string_idl`
;
; Called by :ref:`traite_amsuab.sh`
;
; Previous step : :ref:`extract_amsu.pro`
;
; Next step : ... science !
;
; TODO
; ====
;
;
; EVOLUTIONS
; ==========
;
;
;-
PRO amsu2ncdf, numch, yyyy, mm, dd, lon_min, lon_max, lat_min, lat_max, resol
;
compile_opt idl2, strictarrsubs
;
@cm_project


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)

file='/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'

testfilename=FILE_SEARCH(file)

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 :
; --------------------------------
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 :
  ; ----------------------------

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/' $
   + 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
  ;
  ; Visualisation :
  ; ---------------
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(tbg[*,*,1],Startx,Starty, COMPRESS=1, $
    LATMIN=minlat, LONMIN=minlon, $
    LATMAX=maxlat, LONMAX=maxlon, $
    max_value=280., min_value=235., 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, 'test_b3_05.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:

end