PRO example,i

CASE i OF
  1: resu=read_parasol(PARAM='Cloud Cover',FILEIN='/Produits_POLDER/P3L2TRGB050617/P3L2TRGB013044GL',LAT=lat,LON=lon)
  2: resu=read_parasol(PARAM='Cloud Cover',FILEIN='/Produits_POLDER/P3L2TRGB050617/P3L2TRGB013044GL',/DAILY,REGION=[40., 50.,140.,170.],LAT=lat,LON=lon)
  3: resu=read_parasol(PARAM='Cor O2 pressure',FILEIN='/Produits_POLDER/P3L2TRGB050617/P3L2TRGB013044GL',/DAILY,REGION=[40., 50.,140.,170.],LAT=lat,LON=lon)
ELSE : resu=read_parasol(PARAM='Cloud Cover',FILEIN='/Produits_POLDER/P3L2TRGB050617/P3L2TRGB013044GL',/DAILY,REGION=[40., 50.,140.,170.],LAT=lat,LON=lon)
ENDCASE
  
END

;==================================================================================================================
; Cette routine permet de lire facilement un produit Parasol sans se préocuper du format.  Ce n'est pas le plus efficace
;  mais ca peut etre mis entre toutes les mains
;
;resu = read_parasol,FILEIN=filein, PARAM=param, REGION=region, DAILY=daily, LAT=lat,LON=lon,ERROR=err
;
; FILEIN = <chemin>   :  Chemin vers le fichier en entree (LeaderFile ou Datafile)
; PARAM  = 'parametre':  Nom du parametre (defini dans les defines anapol).  Si non valide, le programme imprime
;                            la liste des parametres valide pour ce type de fichier
; REGION=[latmin,latmax,lonmin,lonmax]
;     ou [latcen,loncen]   : Optionnel.  Permet de récupérer les données d'une région réduite.  Si 4 éléments, on va chercher
;                           tous les pixels dans la région.  Si deux éléments, on prend le pixel qui correspond 
;                                            longitude entre -180 et 180
; /DAILY              :  Permet de faire facilement un "daily browse".  Si ce parametre est défini, le code cherche
;                         tous les fichiers de meme type dans le dossier, et genere une image pour tous
; LATITUDE =lat       : Pour récupérer les latitudes
; LONGITUDE=lon       : pour récupérer les longitudes
; NPIX = npix         : Nombre de pixels valides.  0 si aucun.  En cas de probleme a l'exécution, La routine retourne 
;                       un nombre négatif (-1 a -6) qui indique le type de probleme
;==================================================================================================================
FUNCTION read_parasol, FILEIN=filein, PARAM=param, REGION=region, DAILY=daily, LATITUDE=Tlat,LONGITUDE=Tlon,ERROR=err

Npix = 0

IF NOT KEYWORD_SET(param) THEN BEGIN
  PRINT, 'Il faut choisir un parametre grace au keyword "PARAM" '
  RETURN, -1
ENDIF
IF NOT KEYWORD_SET(region) THEN region=[-90.,90.,-180.,180.]

seppos = STRPOS(filein,path_sep(),/REVERSE_SEARCH)  ; Position du separateur des dossiers
pathin = STRMID(filein,0,seppos+1)
orbit_label =STRMID(filein,seppos+1,15)
fileD = filein & STRPUT,fileD,'D',seppos+16  ; Data   file
fileL = filein & STRPUT,fileL,'L',seppos+16  ; Leader file

; Ouverture du fichier
OPENR,lunL,fileL,/get_Lun,ERROR=err
IF err NE 0 THEN BEGIN
  user_response=DIALOG_MESSAGE('The Leader file for this product was not found: '+fileD)
  RETURN, -2
ENDIF
;
; Verification rapide que on lit bien un fichier POLDER
;
POINT_LUN,lunL,8
ref = BYTARR(11) & READU,lunL,ref & ref=STRING(ref)
IF (ref NE 'PAST33131CN' AND ref NE 'P2ST33131CN' AND ref NE 'SPG9N122-31') THEN BEGIN
  user_response=DIALOG_MESSAGE('The file was not recognized as a standard POLDER product')
  RETURN, -3
ENDIF
POINT_LUN,lunL,36
typ_prod = BYTARR(8)  &  READU,lunL,typ_prod & typ_prod = STRING(typ_prod)
parasol = STRMID(typ_prod,1,1) EQ '3'
IF Parasol THEN print,' Parasol data file' ELSE print,'POLDER data file'
PRINT,' Product type : ',typ_prod
IF STRMID(typ_prod,3,1) EQ '3' THEN daily = 0  ;  On ne peut pas faire du daily sur un produit de niveau 3
;
; On lit le nom du fichier (en interne), pour en deduire le type
;
OPENR,lunD,fileD,/get_Lun,ERROR=err
IF err NE 0 THEN BEGIN
 user_response=DIALOG_MESSAGE('Error when trying to open file '+fileD)
 RETURN, -4
ENDIF
; On lit la taille du record (en interne)
POINT_LUN,lunD,58
len_prod = BYTARR(2) & readu,lunD,len_prod
rec_size=len_prod[0]*256l+len_prod[1]

Prod = define_product(typ_prod,rec_size,ERROR=err,PARASOL=parasol)  ; On recupere toutes les infos sur ce produit

FREE_LUN, lunL,lunD
;
; Recherche dans la liste du parametre choisi
;
vu = WHERE(prod.U_par.name EQ param,ccu)
vd = WHERE(prod.D_par.name EQ param,ccd)
IF ccu+ccd EQ 0 THEN BEGIN
  PRINT, 'Parametre non reconnu: ',param
  PRINT, 'Parametres valides:'
  FOR ip=0, prod.Nu_par-1 DO print," >"+prod.U_par[ip].name+"<"
  FOR ip=0, prod.Nd_par-1 DO print," >"+prod.D_par[ip].name+"<"
  RETURN, -5
ENDIF
IF ccu GT 0 THEN BEGIN
  ip = vu[0]
  par = prod.U_par[ip]
  Ndir = 0
  pos = par.fi_by
ENDIF ELSE BEGIN
  ip = vd[0]
  par = prod.D_par[ip]
  IF parasol THEN Ndir=16 ELSE Ndir=14
  pos = Nbyte_nd + INDGEN(Ndir)*Nbyte_d + par.fi_by
ENDELSE
;
; On cherche la zone
;
IF N_ELEMENTS(region) EQ 4 THEN BEGIN
  linmin = ROUND((90.-region[1])*prod.resol) > 0
  linmax = ROUND((90.-region[0])*prod.resol) < (180*prod.resol -1)
ENDIF ELSE BEGIN
  linmin = ROUND((90.-region[0])*prod.resol)
  linmax = linmin
  lat = 90. - (linmin+0.5)/prod.resol
  NcolM = ROUND(180.*prod.resol*COS(lat*!DTOR))
  colpix = FIX(region[1]*NcolM/180. + 180.*prod.resol) + 1
ENDELSE
;help,par,/struc

IF KEYWORD_SET(daily) THEN liste = FILE_SEARCH(pathin,STRMID(orbit_label,0,8)+'*D') $
                      ELSE liste = [filein]

FOR ifile = 0, N_ELEMENTS(liste)-1 DO BEGIN
  fileName = STRMID(liste[ifile],0,STRLEN(liste[ifile])-1)
  OPENR,lunL,fileName+'L',/GET_LUN
  OPENR,lunD,fileName+'D',/GET_LUN
;-----------------------------------------------------------
;  Lectures
  GetNpix,lunL,Prod,NPIXLINE=Npix,NPIXSUM=NpixSum,PARASOL=parasol
;
;  Boucle sur les lignes
;
  FOR ilin = linmin,linmax DO IF Npix[ilin] GT 0 THEN BEGIN
    point_lun,lunD,180+NpixSum[ilin]*rec_size
    data = BYTARR(rec_size,  Npix[ilin])
    readu,lunD,data
    col    = REFORM(data[ 8,*]*256s + data[ 9,*])
    IF N_ELEMENTS(region) EQ 4 THEN BEGIN
      lat = 90. - (ilin+0.5)/prod.resol
      NcolM = ROUND(180.*prod.resol*COS(lat*!DTOR))
      lon = (col-180*prod.resol-0.5)*180./NcolM
      valid = WHERE(lon GE region[2] AND lon LE region[3],cc)
    ENDIF ELSE BEGIN
      valid = WHERE(col EQ colpix,cc)
    ENDELSE
    IF cc EQ 0 THEN CONTINUE
    CASE par.size OF
      1 : resu = REFORM(data[pos,valid])
      2 : resu = REFORM(data[pos,valid])*256us + REFORM(data[pos+1,valid])
     -2 : resu = REFORM(data[pos,valid])*256s  + REFORM(data[pos+1,valid])
    ENDCASE
    IF N_ELEMENTS(Tresu) EQ 0 THEN BEGIN
      Tlat = REPLICATE(lat,cc)
      Tlon = lon[valid]
      Tresu = resu
    ENDIF ELSE BEGIN
      Tlat  = [Tlat,REPLICATE(lat,cc)]
      Tlon  = [Tlon,lon[valid]]
      Tresu = [Tresu,resu]
    ENDELSE
  ENDIF   ; Fin de la boucle sur les lignes
  FREE_LUN,lunL,lunD
ENDFOR ; Fin de la boucle sur les fichiers

Npix = N_ELEMENTS(Tlat) & IF Npix EQ 0 THEN RETURN,-6
IF Ndir GT 0 THEN Tresu = REFORM(Tresu,Ndir,Npix)

RETURN, Tresu*par.slope + par.offset
END

;==================================================================
;  Cette routine retourne Npix et NpixSum qui permettent de lire rapidement le fichier ligne par ligne
;==================================================================
PRO GetNpix,lunL,Prod,NPIXLINE=Npix,NPIXSUM=NpixSum,PARASOL=parasol

; On lit le nombre de pixels par ligne
Nlin = 180 * Prod.resol  &  Npix=INTARR(Nlin)
CASE STRMID(Prod.type,1,1) OF
 '1' : point_lun,lunL,180+360+1620+180+166320l+720+13140+204
 '2' : point_lun,lunL,180+360+1620+180        +720+13140+204
 '3' : point_lun,lunL,180+360                 +720+13140+204
ENDCASE
READF,lunL,Npix,format='(3240I4)'
;
; On trouve le nombre de pixels par ligne
;
NpixSum = LONARR(Nlin)
IF parasol AND prod.type EQ 'L1TBG1' THEN BEGIN
  NpixSum[Nlin-1] = 0
  FOR i = Nlin-2, 0, -1 DO NpixSum[i] = NpixSum[i+1]+ Npix[i+1]
ENDIF ELSE BEGIN
  NpixSum[0] = 0
  FOR i = 1, Nlin-1 DO NpixSum[i] = NpixSum[i-1]+ Npix[i-1]
ENDELSE
END