; ; Ce code permet de comparer tres rapidement les estimations aerosols obtenus par differents satellites avec les mesures aeronet. ; En entree, on specifie le repertoire ou sont ranges les mesures aeronet, celui ou sont les produits satellite ; et une distance maximale pour la comparaison. ; Le code genere des plots et un fichiers de resultats facilement exploitables ; Si on a deja fait les extractions Aeronet, utiliser l'option /RESTORE (plus rapide) ; Si on a deja fait les coincidences, prendre l'option /PLOT (Beaucoup plus rapide) ; ; Exemple : valid_parasol_aeronet,/RESTORE ; permet de faire les stats sur les pixels "Terres", dans le cas ou on a deja fait les extractions Aeronet ; L'option /VERBOSE permet de faire des prints pour controler ce qui se passe ; SATELLITE = ParasolTE ou ParasolOC ou MODIS_OC ou MODIS_TE ou MERIS_OC ou MERIS_TE ou SEVIRI ou CALIPSO ; SAVE_STATS= permet d'enregistrer les resultats statistiques dans un fichier ; /NOPLOTS pour ne pas faire de plots (utile pour lancer en batch) ; /VERBOSE pour des affichages intermediaires (debug) ; DISTMAX=val pour specifier la distance maximale (km, defaut=50) ; DTMAX=val pour specifier la distance temporelle maximale (minutes, defaut=60) ; /MAPS pour tracer des cartes des résultats, en plus des scatter ; plots ; ; Les resultats numeriques (coincidences entre aeronet et satellite) ; sont conserves dans une structure accessible par RESTORE dans le ; fichier Coincidences.sav ; ; Il est necessaire de specifier la periode. Voir exemple ci-dessous (Period=) ; ; Francois-Marie Breon, LSCE, Juillet 2010 ; fmbreon@cea.fr ; PRO valid_satellite_aeronet,RESTORE=restore,PLOTS=plots,VERBOSE=verbose,SATELLITE=satellite,SAVE_STATS=savestats,NOPLOTS=noplots, $ DISTMAX=distmax, DTMAX=dtmax,MAPS=maps FORWARD_FUNCTION TOstat,StringPos,read_aeronet,readMODIS, readCALIPSO,readMERIS,readSEVIRI,readAATSR,readParasol,readATSR,STMIN,Get_Location,get_range COMMON COM_MASQUE,masque,mult ; !EXCEPT=1 IF NOT KEYWORD_SET(satellite) THEN satellite = 'CALIPSO' IF NOT KEYWORD_SET(distmax) THEN distmax = 50. ; Distance maximale acceptable en km IF NOT KEYWORD_SET( dTmax) THEN dTMax = 60. ; Difference temporelle acceptable, en minutes ;HomePath = '/Users/fmbreon/Desktop/GEOMON/' HomePath = '/home/breon/GEOMON/' ;HomePath = '/scratch2/vermeulen/GEOMON/' ;Period = '2004/2004_1*/' ; Periode choisie pour les fichiers satellite Period = '200*/*' ; Periode choisie pour les fichiers satellite ;Period = '2006/2006_*' ; Periode choisie pour les fichiers satellite ;pathAeronet = '/scratch2/vermeulen/DATA_AERONET/AOTCF_LEV20/SINCE_2002/' ; Chemin vers fichiers aeronet pathAeronet = '/scratch2/vermeulen/DATA_AERONET/AOTCF_LEV20_20100626/SINCE_2002/' ; Chemin vers fichiers aeronet ;pathAeronet = '/scratch2/vermeulen/DATA_AERONET/AOT/LEV20/TEST2/' ; Chemin vers fichiers aeronet pathSave = HomePath + 'SAVE_FILES/' + satellite +'/' ; Chemin vers les fichiers intermediaires pathOut = HomePath + 'PLOTS/' + satellite +'/' ; Chemin vers les plots pathStat = HomePath + 'STATS/' + satellite +'/' ; Chemin vers les statistiques fileAeronet = pathSave + 'Aeronet.sav' ; Fichier contenant les obs aeronet fileCoincidences = pathSave + 'Coincidences.sav' ; Fichiers contenant toutes les coindicences fileTempo = HomePath + 'tempo.dat' print,'SATELLITE: ',satellite satname = '' te = 0 CASE satellite OF 'ParasolTE' : BEGIN pathSatellite = '/DATA/LIENS/PARASOL/LS2/'+Period ; Chemin vers les fichiers N2 Parasol TE satname = 'parasol' te = 1 ytit = 'Parasol Land' END 'ParasolOC' : BEGIN pathSatellite = '/DATA/LIENS/PARASOL/OC2/'+Period ; Chemin vers les fichiers N2 Parasol OC satname = 'parasol' ytit = 'Parasol Ocean' END 'MODIS_OC' : BEGIN pathSatellite = '/DATA/LIENS/MODIS/MYD04_L2/'+Period ; Chemin vers les fichiers N2 Modis satname = 'modis' ytit = 'MODIS Ocean' END 'MODIS_TE' : BEGIN pathSatellite = '/DATA/LIENS/MODIS/MYD04_L2/'+Period ; Chemin vers les fichiers N2 Modis satname = 'modis' te = 1 ytit = 'MODIS Land' END 'MERIS_OC' : BEGIN pathSatellite = '/DATA/LIENS/MERIS/MER_RR__2P/'+Period ; Chemin vers les fichiers N2 Meris satname = 'meris' ytit = 'MERIS Ocean' EPR_INIT_API END 'MERIS_TE' : BEGIN pathSatellite = '/DATA/LIENS/MERIS/MER_RR__2P/'+Period ; Chemin vers les fichiers N2 Meris satname = 'meris' te = 1 ytit = 'MERIS Land' EPR_INIT_API END 'SEVIRI' : BEGIN pathSatellite = '/DATA/LIENS/SEVIRI/AER_OC/'+Period ; Chemin vers les fichiers SEVIRI satname = 'seviri' ytit = 'SEVIRI Ocean' Define_Seviri_latlon END 'SEVIRI_TE': BEGIN pathSatellite = '/DATA/LIENS/SEVIRI/SMAOL_AOT.v1.3.2/'+Period ; Chemin vers les fichiers SEVIRI satname = 'seviri' ytit = 'SEVIRI Land' te = 1 Define_Seviri_latlon END 'CALIPSO' : BEGIN pathSatellite = '/DATA/LIENS/CALIOP/05kmALay/'+Period ; Chemin vers les fichiers CALIPSO satname = 'calipso' ytit = 'CALIPSO' distmax = 150. END 'AATSR' : BEGIN pathSatellite = '/home/breon/AATSR/' ; Chemin vers les fichiers AATSR satname = 'aatsr' ytit = 'AATSR' te = 1 END ELSE : BEGIN PRINT,'SATELLITE must be "ParasolOC" or "ParasolTE" or "MODIS_OC" or "MODIS_TE" '+$ 'or "MERIS_TE" or "MERIS_OC" or "SEVIRI" or "CALIPSO" or "AATSR"' RETURN END ENDCASE IF NOT FILE_TEST(pathSave,/DIR) THEN FILE_MKDIR,pathSave IF NOT FILE_TEST(pathout ,/DIR) THEN FILE_MKDIR,pathout IF NOT FILE_TEST(pathStat,/DIR) THEN FILE_MKDIR,pathStat IF KEYWORD_SET(verbose) THEN BEGIN PRINT,Distmax,FORMAT='("DistMax:",F5.0," km ")' PRINT,dtMax ,FORMAT='("D_T_Max:",F5.0," min")' PRINT,' pathAeronet: ',pathAeronet PRINT,'pathSatellite: ',pathSatellite ENDIF IF KEYWORD_SET(plots) THEN GOTO, jump_plots IF KEYWORD_SET(maps ) THEN GOTO, jump_maps ;==================================================================================================================== ; On fait la liste des fichiers Parasol disponibles. Permet de definir le debut et la fin de la periode d'analyse ;==================================================================================================================== CASE satellite OF 'ParasolTE' : listeSatellite = FILE_SEARCH(pathSatellite,'P3L2TLGC*D', COUNT=Norb) 'ParasolOC' : listeSatellite = FILE_SEARCH(pathSatellite,'P3L2TOGC*D', COUNT=Norb) 'MODIS_OC' : listeSatellite = FILE_SEARCH(pathSatellite,'MYD04*.hdf', COUNT=Norb) 'MODIS_TE' : listeSatellite = FILE_SEARCH(pathSatellite,'MYD04*.hdf', COUNT=Norb) 'MERIS_OC' : listeSatellite = FILE_SEARCH(pathSatellite,'MER*' , COUNT=Norb) 'MERIS_TE' : listeSatellite = FILE_SEARCH(pathSatellite,'MER*' , COUNT=Norb) 'SEVIRI' : listeSatellite = FILE_SEARCH(pathSatellite,'SEVIRI_AER_OC_*.*.hdf', COUNT=Norb) 'SEVIRI_TE': listeSatellite = FILE_SEARCH(pathSatellite,'SMAOL-AOT*.hdf', COUNT=Norb) 'CALIPSO' : listeSatellite = FILE_SEARCH(pathSatellite, 'CAL_LID_L2_05kmALay*.hdf', COUNT=Norb) 'AATSR' : listeSatellite = FILE_SEARCH(pathSatellite, '*.hdf', COUNT=Norb) ENDCASE IF KEYWORD_SET(verbose) or Norb EQ 0 THEN print,'Nombre de fichiers Satellite trouves:',Norb IF Norb EQ 0 THEN RETURN ; ; On va lire la date du premier et du dernier pour definir la periode ; IF satname EQ 'aatsr' THEN deb=1 ELSE deb=-10 pp = STRPOS(listeSatellite[ 0],PATH_SEP(),/REVERSE_SEARCH) & date=STRMID(listeSatellite[ 0],pp+deb,10) date = STRSPLIT(date,'_-',/EXTRACT) & J0 = JULDAY(date[1],date[2],date[0], 0., 0., 0.) pp = STRPOS(listeSatellite[Norb-1],PATH_SEP(),/REVERSE_SEARCH) & date=STRMID(listeSatellite[Norb-1],pp+deb,10) date = STRSPLIT(date,'_-',/EXTRACT) & Jlast = JULDAY(date[1],date[2],date[0],23.,59.,59.) J0 = J0 - (30+dTmax)/60./24. ; On prend un peu de marge pour l'extraction des donnÃ©es Aeronet Jlast = Jlast + (30+dTmax)/60./24. CALDAT,J0 ,minmon,minday,minyear ; Conversion en jour/mois/annee CALDAT,Jlast,maxmon,maxday,maxyear PRINT,minday,minmon,minyear,maxday,maxmon,maxyear,Norb,$ FORMAT='("Satellite: Du ",I2,"/",I2.2,"/",I4," au ",I2,"/",I2.2,"/",I4," Nb orbites:",I6)' IF KEYWORD_SET(restore) THEN GOTO,jump_restore ;==================================================================================================================== ; Lecture des donnees Aeronet ;==================================================================================================================== liste = FILE_SEARCH(pathAeronet+'*',COUNT=Nstations) ; On cherche la liste des fichiers aeronet IF KEYWORD_SET(verbose) THEN PRINT,'Nombre de stations aeronet:',Nstations Nstations = 0 FOR i=0, N_ELEMENTS(liste)-1 DO BEGIN ; Boucle sur les fichiers aeronet infos=read_aeronet(liste[i],DATERANGE=[J0,Jlast],VARNAMES=varnames,VERBOSE=verbose) ; Lecture du fichier aeronet if i eq 0 AND KEYWORD_SET(verbose) THEN print,varnames ; Pour Ã©crire les variabes lues dans le premier fichier IF infos.Nobs EQ 0 THEN CONTINUE IF Nstations EQ 0 THEN stations = REPLICATE(infos,N_ELEMENTS(liste)) ELSE stations[Nstations] = infos Nstations = Nstations+1 ENDFOR print,'Number of Aeronet Stations with measurements within time range: ', Nstations IF Nstations EQ 0 THEN RETURN stations = stations[0:Nstations-1] PRINT,'Now build masque' build_masque,stations,distmax,VERBOSE=verbose PRINT,'masque is computed.' ;==================================================================================================================== ; Les lectures Aeronet sont terminees. On a donc range dans une structure tout ce qui est potentiellement utile ;==================================================================================================================== SAVE,stations,varnames,masque,mult,FILENAME=fileAeronet jump_restore: IF KEYWORD_SET(restore) THEN RESTORE,FILENAME=fileAeronet IF KEYWORD_SET(verbose) THEN PRINT,'now search coincidences' ;============================================================================================================ ; Les coincidences sont un tableau de structure. Dans chaque structure, on a le jour/heure, le numÃ©ro de la station, ; le dÃ©calage temporel satellite/sol, la distance spatiale, les Ã©paisseurs optiques aeronet (Azzz) et les Ã©paisseurs ; optiques Parasol (Szzz) ; Chaque "obs" est en fait une structure qui donne d'une part la valeur la plus proche en spatial/temporel, et d'autre ; part la moyenne et ecart-type des valeurs dans la boite spatio-temporelle. ; obs = {obs, Nobs:0, one:0., Mean:0., Sigma:0.} coincidence = {coincidence,Jloc:0l,ista:0,deltaT:0.,dist:0.,Nsat:0,Naeronet:0,$ ATo500:obs, ATo675:obs,ATo870:obs,AToF55:obs,AToF5M:obs,$ STo500:obs,Sto670:obs,Sto865:obs,StoF55:obs,Squal:obs} Ncoi = 0l ; Ncoi est le nombre de coincidences ito870 = WHERE(varnames EQ 'AOT_870') ; Position des donnees interessantes dans le tableau stations.Mes ito675 = WHERE(varnames EQ 'AOT_675') ito500 = WHERE(varnames EQ 'AOT_500') itoFine= WHERE(varnames EQ 'AOT_F550') ; Estimation de l'ep optique Fine, avec Algo Loeb itoFinM= WHERE(varnames EQ 'AOT_Fine') ; Estimation de l'ep optique Fine avec mon algo ;============================================================================================================ ; On attaque maintenant la lecture des donnees Satellite ; FOR ifile = 0l, N_ELEMENTS(listeSatellite)-1 DO BEGIN IF FILE_TEST(listeSatellite[ifile],/ZERO_LENGTH) THEN CONTINUE CASE satname OF 'parasol' : TSat = readParasol(FILEIN=listeSatellite[ifile],Nv=Nv, VERBOSE=verbose) 'modis' : Tsat = readMODIS (FILEIN=listeSatellite[ifile],Nv=Nv,TE=te,VERBOSE=verbose) 'meris' : Tsat = readMERIS (FILEIN=listeSatellite[ifile],Nv=Nv,TE=te,VERBOSE=verbose,fileTempo=fileTempo) 'seviri' : Tsat = readSEVIRI (FILEIN=listeSatellite[ifile],Nv=Nv ,VERBOSE=verbose) 'calipso' : Tsat = readCALIPSO(FILEIN=listeSatellite[ifile],Nv=Nv ,VERBOSE=verbose) 'aatsr' : Tsat = readAATSR (FILEIN=listeSatellite[ifile],Nv=Nv ,VERBOSE=verbose) ENDCASE IF Nv EQ 0 THEN CONTINUE Duplicate_SatObs,TSat,Nv ; On duplique les obs satellite vues par plusieurs stations sta_indice = Tsat.ista ista = sta_indice[UNIQ(sta_indice,SORT(sta_indice))] ; Donne le numeros des stations vus par cette orbite IF KEYWORD_SET(verbose) THEN PRINT, N_ELEMENTS(ista),ista+1,FORMAT='(I3," stations vues par cette orbite. Numeros:",1000I4)' FOR is = 0, N_ELEMENTS(ista)-1 DO BEGIN v = WHERE(TSat.ista EQ ista[is]) ; Analyse des coincidences. On commence par etudier les mesures aeronet, pour verifier si il y a au moins ; une obs dans l'intervalle Jobs = TSat[v[0]].J DeltaT = ABS(Jobs - *(stations[ista[is]].Jdays))*24.*60 ; Difference temporelle, en minutes valid = WHERE(DeltaT LT DTmax,Naeronet) IF KEYWORD_SET(verbose) THEN PRINT,stations[ista[is]].location,Naeronet,STMIN(MIN(deltaT)), $ format='(A15,":",I3," mesures dans l intervalle temporel. ",A)' IF Naeronet EQ 0 THEN CONTINUE s = SORT(DeltaT[valid]) dTmin = deltaT[valid[s[0]]] Aobs = (*stations[ista[is]].Mes)[*,valid[s]] ; Aobs donne les observations aeronet, classees par ecart en temps croissant ; On etudie maintenant les mesures Satellite CosLat = COS(stations[ista[is]].lat*!DTOR) dist = SQRT((stations[ista[is]].lat-TSat[v].lat)^2 + ((stations[ista[is]].long-TSat[v].lon)*CosLat)^2)*40000./360. s = SORT(dist) ; v donne les indices des observations parasol qui sont proches de la station Sat = TSat[v[s]] ; classees par distance croissante MinDist = dist[s[0]] Jloc = ROUND(Jobs + stations[ista[is]].long/360.) ; Calcul du jour "local" coincidence={coincidence,Jloc:Jloc,ista:ista[is],deltaT:dTmin,dist:MinDist,$ ; On remplit le tableau des coincidences Nsat:N_ELEMENTS(v),Naeronet:Naeronet,$ ATo500:TOstat(Aobs[ito500,*]) ,$ Ato675:TOstat(Aobs[ito675,*]) ,$ Ato870:TOstat(Aobs[ito870,*]) ,$ AtoF55:TOstat(Aobs[itoFine,*]),$ AtoF5M:TOstat(Aobs[itoFinM,*]),$ STo500:TOstat(Sat.to500 ,QUAL=Sat.qual),$ STo670:TOstat(Sat.to670 ,QUAL=Sat.qual),$ STo865:TOstat(Sat.to865 ,QUAL=Sat.qual),$ StoF55:TOstat(Sat.toF550,QUAL=Sat.qual),$ Squal :TOstat(Sat.qual)} IF Ncoi EQ 0 THEN coincidences = coincidence ELSE coincidences = [coincidences, coincidence] ; print,coincidences[Ncoi] Ncoi = Ncoi+1 ENDFOR ; Fin de la boucle sur les stations aeronet ENDFOR ; Fin de la boucle sur les fichiers Parasol print,'Nombre de coincidences:', Ncoi IF Ncoi EQ 0 THEN RETURN aeronet = {aeronet,name:'',lat:0.,lon:0.,elev:0.} aeronet = REPLICATE(aeronet,N_ELEMENTS(stations)) aeronet.name = stations.Location & aeronet.lat = stations.lat & aeronet.lon=stations.long & aeronet.elev=stations.elev SAVE,coincidences,aeronet,FILE=fileCoincidences ; On enregistre le fichier de coincidences pour pouvoir refaire des plots ; tres rapidement IF KEYWORD_SET(noplots) THEN RETURN ; Si on ne fait pas de plots (utile en mode batch), on a termine jump_plots: IF KEYWORD_SET(plots) THEN BEGIN RESTORE,FILE=fileCoincidences valid = WHERE(coincidences.deltaT LE dTmax AND coincidences.dist LE distMax) coincidences = coincidences[valid] ENDIF IF KEYWORD_SET(savestats) THEN BEGIN s55F=21 & OPENW,s55F,pathStat+'statsM_550F.txt' s500=22 & OPENW,s500,pathStat+'statsM_500.txt' s670=23 & OPENW,s670,pathStat+'statsM_670.txt' s865=24 & OPENW,s865,pathStat+'statsM_865.txt' PRINTF,s55F,' Num Name Nobs Corr Bias RMS slope intercept Frac_valid' PRINTF,s500,' Num Name Nobs Corr Bias RMS slope intercept Frac_valid' PRINTF,s670,' Num Name Nobs Corr Bias RMS slope intercept Frac_valid' PRINTF,s865,' Num Name Nobs Corr Bias RMS slope intercept Frac_valid' ENDIF !P.CHARSIZE=1. mypalette,NC=NC Np=50 vmin=0.0001 & vmax=1. siz = 250 ;WINDOW,1,XSIZE=3*siz,YSIZE=2*siz,RETAIN=2 WINDOW,1,XSIZE=3*siz,YSIZE=siz,RETAIN=2 FOR ista = 0, N_ELEMENTS(aeronet) DO BEGIN ; Boucle sur les stations IF ista EQ 0 THEN BEGIN ; On commence par faire un plot pour toutes les stations lat = aeronet[coincidences.ista].lat lon = aeronet[coincidences.ista].lon desert = (lat GT 0. AND lat LT 25. AND lon GT -10. AND lon LT 50.) ; On fait un cas spÃ©cial pour les dÃ©serts v=WHERE(coincidences.Squal.one GT 0.01 AND coincidences.AtoF55.one GT 0. AND NOT desert) a=WHERE(coincidences.AtoF55.Mean GT 0. ) tit = 'AAll_Stations' & out=tit stn = '0000 All_Stations ' ENDIF ELSE BEGIN ; puis station par station v = WHERE(coincidences.ista EQ (ista-1) AND coincidences.Squal.one GT 0.01 AND coincidences.AtoF55.one GT 0. ,cc) a = WHERE(coincidences.ista EQ (ista-1) ,cc) IF cc EQ 0 THEN CONTINUE tit = aeronet[ista-1].name+' '+StringPos(aeronet[ista-1].lat,aeronet[ista-1].lon,aeronet[ista-1].elev) out = aeronet[ista-1].name stn = STRING(ista,aeronet[ista-1].name,FORMAT='(I4,x,A20)') ENDELSE IF N_ELEMENTS(v) EQ 1 THEN IF v EQ -1 THEN v=[0] ; Pour eviter les problemes en cas de pas de points IF N_ELEMENTS(a) EQ 1 THEN IF a EQ -1 THEN a=[0] DEVICE,DECOMPOSE=0 & ERASE ; scatter,coincidences[a].AtoF55.Mean ,coincidences[a].StoF55.Mean ,XTIT='AerFine 550 nm',YTIT=ytit,NC=NC,POS=[0.00,0.00,0.00,0.00],FSAVE=s55F,STN=stn,TE=te ; scatter,coincidences[a].Ato500.Mean ,coincidences[a].Sto500.Mean ,XTIT='Aeronet 500 nm', NC=NC,POS=[0.00,0.00,0.00,0.00],FSAVE=s500,STN=stn,TE=te ; scatter,coincidences[a].Ato870.Mean ,coincidences[a].Sto865.Mean ,XTIT='Aeronet 865 nm', NC=NC,POS=[0.00,0.00,0.00,0.00],FSAVE=s865,STN=stn,TE=te ; scatter,coincidences[a].Ato675.Mean ,coincidences[a].Sto670.Mean ,XTIT='Aeronet 670 nm', NC=NC,POS=[0.00,0.00,0.00,0.00],FSAVE=s670,STN=stn,TE=te scatter,coincidences[v].AtoF55.one ,coincidences[v].StoF55.one ,XTIT='AerFine 550 nm',YTIT=ytit,NC=NC,POS=[0.05,0.15,0.33,0.99],FSAVE=s55F,STN=stn,TE=te scatter,coincidences[v].Ato500.one ,coincidences[v].Sto500.one ,XTIT='Aeronet 500 nm', NC=NC,POS=[0.38,0.15,0.66,0.99],FSAVE=s500,STN=stn,TE=te scatter,coincidences[v].Ato870.one ,coincidences[v].Sto865.one ,XTIT='Aeronet 865 nm', NC=NC,POS=[0.71,0.15,0.99,0.99],FSAVE=s865,STN=stn,TE=te scatter,coincidences[v].Ato675.one ,coincidences[v].Sto670.one ,XTIT='Aeronet 670 nm', NC=NC,POS=[0.00,0.00,0.00,0.00],FSAVE=s670,STN=stn,TE=te ; scatter,coincidences[v].AtoF55.one ,coincidences[v].StoF55.one , YTIT=ytit,NC=NC,POS=[0.05,0.57,0.33,0.99],FSAVE=s55F,STN=stn,TE=te ; scatter,coincidences[v].Ato500.one ,coincidences[v].Sto500.one , NC=NC,POS=[0.38,0.57,0.66,0.99],FSAVE=s500,STN=stn,TE=te ; scatter,coincidences[v].Ato870.one ,coincidences[v].Sto865.one , NC=NC,POS=[0.71,0.57,0.99,0.99],FSAVE=s865,STN=stn,TE=te ; scatter,coincidences[a].AtoF55.Mean,coincidences[a].StoF55.Mean,XTIT='AerFine 550 nm',YTIT=ytit,NC=NC,POS=[0.05,0.07,0.33,0.49], TE=te ; scatter,coincidences[a].Ato500.Mean,coincidences[a].Sto500.Mean,XTIT='Aeronet 500 nm', NC=NC,POS=[0.38,0.07,0.66,0.49], TE=te ; scatter,coincidences[a].Ato870.Mean,coincidences[a].Sto865.Mean,XTIT='Aeronet 865 nm', NC=NC,POS=[0.71,0.07,0.99,0.49], TE=te ; xyouts,0.50,0.51,tit,ALIGN=0.5,/NORMAL,CHARSIZE=1.5,COL=NC-3 IF KEYWORD_SET(pathout) THEN wpng,pathout+out+'.PNG' ELSE stop ENDFOR IF KEYWORD_SET(savestats) THEN CLOSE,s55F,s500,s670,s865 IF NOT KEYWORD_SET(maps) THEN RETURN ; ; On va maintenant dessiner des cartes du RMS et de la fraction de ; mesures valides ; jump_maps: ang = INDGEN(9)*45.*!DTOR DEVICE,DECOMPOSE=0 USERSYM,COS(ang),SIN(ang),/FILL mypalette RESTORE,FILE=fileCoincidences res = {num:0,name:'',NObs:0l,Corr:0.,Bias:0.,RMS:0.,slope:0.,inter:0.,frac:0.} & fmt = '(I4,x,A20,I6,F7.3,2F6.2,2x,F6.3,F6.2,F6.3)' ; ; On commence par les stats sur l'ep optique totale ; OPENR,lun,pathstat + 'stats_500.txt',/GET_LUN line='' & READF,lun,line IF NOT EOF(lun) THEN BEGIN READF,lun,res,FORMAT=fmt & Tres=res WHILE NOT EOF(lun) DO BEGIN READF,lun,res,FORMAT=fmt & Tres=[Tres,res] ENDWHILE val = WHERE(Tres.Nobs GT 30, Nstat) Tres = Tres[val] Tres.name = STRTRIM(Tres.name,2) lat = FLTARR(Nstat) & lon = FLTARR(Nstat) names = aeronet.name FOR i=1, Nstat-1 DO BEGIN val = WHERE(Tres[i].name EQ names,cc) IF cc EQ 0 THEN STOP, 'Not Found' IF cc GT 1 THEN STOP,' Probleme' lat[i] = aeronet[val[0]].lat lon[i] = aeronet[val[0]].lon PRINT,Tres[i].name,lat[i],lon[i],format='(A20,2F8.2)' ENDFOR draw_map,lat,lon,Tres.frac*100.,TITLE='% of Accurate Retrievals',FILENAME=pathStat+'map_fracValid.PNG',SATELLITE=satellite draw_map,lat,lon,Tres.RMS ,TITLE='RMS Error' ,FILENAME=pathStat+ 'map_RMS.PNG',SATELLITE=satellite ENDIF FREE_LUN,lun ; ; Meme chose avec le Fine Mode ; OPENR,lun,pathstat + 'stats_550F.txt',/GET_LUN line='' & READF,lun,line IF NOT EOF(lun) THEN BEGIN READF,lun,res,FORMAT=fmt & Tres=res WHILE NOT EOF(lun) DO BEGIN READF,lun,res,FORMAT=fmt & Tres=[Tres,res] ENDWHILE val = WHERE(Tres.Nobs GT 30, Nstat) Tres = Tres[val] Tres.name = STRTRIM(Tres.name,2) lat = FLTARR(Nstat) & lon = FLTARR(Nstat) names = aeronet.name FOR i=1, Nstat-1 DO BEGIN val = WHERE(Tres[i].name EQ names,cc) IF cc EQ 0 THEN STOP, 'Not Found' IF cc GT 1 THEN STOP,' Probleme' lat[i] = aeronet[val[0]].lat lon[i] = aeronet[val[0]].lon ENDFOR draw_map,lat,lon,Tres.frac*100.,TITLE='% of Accurate Retrievals',FILENAME=pathStat+'map_fracValid_FM.PNG',SATELLITE=satellite+' Fine Mode' draw_map,lat,lon,Tres.RMS ,TITLE='RMS Error' ,FILENAME=pathStat+ 'map_RMS_FM.PNG',SATELLITE=satellite+' Fine Mode' ENDIF FREE_LUN,lun END ;============================================================================================================ ; Dessin d'une carte avec un point sur chaque site aeronet et un code couleur ;------------------------------------------------------------------------------------------------------------ PRO draw_map,lat,lon,tab,TITLE=tit,FILENAME=fpath,SATELLITE=satellite mypalette,NC=NC WINDOW,1,XSIZE=600,ysize=250,TIT=tit,RETAIN=2 MAP_SET,0,0,/continents,/NOBORDER,POS=[0,0,1,1],LIMIT=[-60.,-170.,80.,180.] res = get_range(tab*1000., PERCENT=5)/1000. vmin=0. & vmax=res[1] visu = ROUND(((((tab+vmin)>0.)/(vmax-vmin))<1.)*(NC-3)) + 1 FOR i=0, N_ELEMENTS(lat)-1 DO oplot,[lon[i]],[lat[i]],PSYM=8,COL=visu[i],SYMSIZE=1.5 XYOUTS, 0.13,0.05,tit, ALIGN=0.5,/NORMAL XYOUTS, 0.5,0.05,satellite,CHARSIZE=1.5, ALIGN=0.5, /NORMAL drawctbl,pos=[0.12,0.1,0.15,0.55],/VERT,VMIN=vmin,VMAX=vmax,FI=1,NC=nc-3 wpng,fpath END ;============================================================================================================ ; Lecture des produits POLDER/Parasol ; FILEIN : Chemin vers le fichier Data ; MASQUE : Taleau masque qui indique le numÃ©ro de/des station proches (plate carrÃ©, rÃ©sol [2*Ngr,Ngr]) ; MULT : Tableau qui permet de gÃ©rer les coincidences multiples ; En sortie, on Ã "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ;------------------------------------------------------------------------------------------------------------ FUNCTION readParasol,FILEIN=filein,Nv=Nv,VERBOSE=verbose COMMON COM_MASQUE,masque,mult Nv=0 OC = STRPOS(filein,'P3L2TOGC') GT 0 TE = STRPOS(filein,'P3L2TLGC') GT 0 IF OC+TE NE 1 THEN STOP,' Probleme dans la reconnaissance du type de fichier' Get_Location, filein, Jtu=Jeq IF KEYWORD_SET(verbose) THEN BEGIN CALDAT,Jeq,mm,dd,yy,hh,mi,ss PRINT, filein+STRING(dd,mm,yy,hh,mi,ss,FORMAT='(3x,2(I2.2,"/"),I4,I4,2(":",I2.2))') ENDIF OPENR,lun,filein,/GET_LUN,/SWAP_IF_LITTLE ; Ouverture du fichier header = LONARR(180/4) READU,lun, header ; Lecture du header pour lire la taille et le nombre de records Npix = header[13] rec_size = header[14] IF Npix EQ 0 THEN FREE_LUN,lun IF Npix EQ 0 THEN RETURN,0 dat = BYTARR(rec_size,Npix) ; Lecture des donnÃ©es et fermeture du fichier READU,lun,dat FREE_LUN,lun lin = REFORM(dat[ 6,*]*256 + dat[ 7,*]) ; Numero de ligne et colonne dans la grille POLDER col = REFORM(dat[ 8,*]*256 + dat[ 9,*]) lat = -(lin- 0.5)*180./1080. + 90. lon = (col-1080.5)*180/ROUND(1080.*cos(lat*!DTOR)) s = SIZE(masque) & Ngr=s[2] cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels parasol sont au voisinage de stations aeronet valid = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station ; print,'Nombre de pixels proches:',Nc IF Nc EQ 0 THEN RETURN,0 ; Si aucun pixel proche, on passe au fichier suivant dat = dat[*,valid] cc = cc[valid] ll = ll[valid] sta_indice = sta_indice[valid] lin = REFORM(dat[ 6,*]*256 + dat[ 7,*]) col = REFORM(dat[ 8,*]*256 + dat[ 9,*]) Nc = N_ELEMENTS(lin) Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Nc) TSat.ista = sta_indice TSat.lat = -(lin- 0.5)*180./1080. + 90. TSat.lon = (col-1080.5)*180/ROUND(1080.*cos(TSat.lat*!DTOR)) TSat.J = Jeq + TSat.lat/360.*16./233 ; Formule approchÃ©e pour une orbite ascendante IF OC THEN BEGIN ; On decode les parametres en epaisseur optique et coef d'angstrom slopeOC = 2.E-3 TSat.qual = REFORM(dat[17,*])*0.01 TSat.to865 = REFORM(dat[20,*]*256 + dat[21,*])*slopeOC TSat.to670 = REFORM(dat[22,*]*256 + dat[23,*])*slopeOC Ang = REFORM(dat[24,*]*256 + dat[25,*])*0.01 - 0.5 TSat.to500 = TSat.to670 * (670./500.)^Ang toF = REFORM(dat[35,*]*256 + dat[36,*])*slopeOC ; Fine mode optical thickness 670 nm AngF = REFORM(dat[37,*]*256 + dat[38,*])*0.01 - 0.5 ; Fine mode Angstrom exponent TSat.toF550 = toF * (670./550.)^AngF ENDIF ELSE BEGIN slopeTE = 2.E-3 TSat.qual = REFORM(dat[27,*])*0.01 AngF = REFORM(dat[20,*])*0.014 toF = REFORM(dat[17,*]*256 + dat[18,*])*slopeTE ; Fine mode optical thickness 865 nm TSat.toF550 = toF * (865./550.)^AngF ENDELSE Nv = N_ELEMENTS(TSat) RETURN,TSat END ;============================================================================================================ ; Lecture des produits SEVIRI ; FILEIN : Chemin vers le fichier Data ; MASQUE : Taleau masque qui indique le numero de/des station proches (plate carrÃ©, rÃ©sol [2*Ngr,Ngr]) ; MULT : Tableau qui permet de gÃ©rer les coincidences multiples ; En sortie, on Ã "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ;------------------------------------------------------------------------------------------------------------ FUNCTION readSEVIRI,FILEIN=filein,Nv=Nv,VERBOSE=verbose FORWARD_FUNCTION GETHDF COMMON Seviri_latlon,Sev_lat,Sev_lon COMMON COM_MASQUE,masque,mult to_unvalid = 252 ; This value and higher for unvalid fac_to = 2. /251 ; Facteur pour convertir les ep.optiques en valeurs physiques fac_ang = 3.1/251 ; Facteur pour convertir les coef d'angstrom en valeurs physiques ang_offset = -0.1 Lambda0 = 550. Nv=0 IF NOT HDF_ISHDF(filein) THEN RETURN,0 ; On vérifie que c'est bien un fichier HDF datetime = STRMID(filein,STRLEN(filein)-17,12) yy=0 & mm=0 & dd=0 & hh=0 & mi=0 READS,datetime,yy,mm,dd,hh,mi,FORMAT='(I4,2I2,x,2I2)' JulObs = JULDAY(mm,dd,yy,hh,mi,0.) to = GETHDF(filein=filein,PARAMETER='Aerosol_Optical_Depth',/BINARY) ang = GETHDF(filein=filein,PARAMETER='Angstrom_Exponent' ,/BINARY) ;to = BYTARR(3712,3712) ;ang = BYTARR(3712,3712) ;OPENR,lun,filein,/GET_LUN ;READU,lun,to ;FREE_LUN,lun ;STRPUT,filein,'ang_over_sea',STRPOS(filein,'tau_over_sea') ;OPENR,lun,filein,/GET_LUN ;READU,lun,ang ;FREE_LUN,lun valid = WHERE(to LT to_unvalid,nn) IF nn EQ 0 THEN RETURN,0 lat = Sev_lat[valid] ; Donner ici les lats et lons SEVIRI lon = Sev_lon[valid] s = SIZE(masque) & Ngr=s[2] cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels parasol sont au voisinage de stations aeronet vv = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station IF KEYWORD_SET(verbose) THEN print,'Nombre de pixels proches:',Nc IF Nc EQ 0 THEN RETURN,0 ; Si aucun pixel proche, on passe au fichier suivant to = to[valid[vv]] ang = ang[valid[vv]] cc = cc[vv] ll = ll[vv] lat = lat[vv] lon = lon[vv] sta_indice = sta_indice[vv] Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:1.} Tsat = REPLICATE(Sat,N_ELEMENTS(to)) Tsat.ista = sta_indice Tsat.lat = lat Tsat.lon = lon Tsat.J = JulObs Tsat.to500 = (to*fac_to)*(Lambda0/500.)^(ang*fac_ang+ang_offset) Tsat.to670 = (to*fac_to)*(Lambda0/670.)^(ang*fac_ang+ang_offset) Tsat.to865 = (to*fac_to)*(Lambda0/865.)^(ang*fac_ang+ang_offset) Nv = N_ELEMENTS(TSat) RETURN,TSat END ;============================================================================================================ ; Lecture des produits SEVIRI Terres-Emergées ; FILEIN : Chemin vers le fichier hdf ; MASQUE : Taleau masque qui indique le numero de/des station proches (plate carre, resol [2*Ngr,Ngr]) ; MULT : Tableau qui permet de gerer les coincidences multiples ; En sortie, on a "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ;------------------------------------------------------------------------------------------------------------ FUNCTION readSEVIRI_TE,FILEIN=filein,Nv=Nv,VERBOSE=verbose FORWARD_FUNCTION GETHDF COMMON Seviri_latlon,Sev_lat,Sev_lon COMMON COM_MASQUE,masque,mult IF NOT HDF_ISHDF(filein) THEN RETURN,0 ; On vérifie que c'est bien un fichier HDF Nv=0 s = SIZE(masque) & Ngr=s[2] datetime = STRMID(filein,STRLEN(filein)-16,12) READS,datetime,yy,mm,dd,hh,mi,FORMAT='(I4,4I2)' JulObs = JULDAY(mm,dd,yy,hh,mi,0.) IF KEYWORD_SET(verbose) THEN BEGIN CALDAT,JulObs,mm,dd,yy,hh,mi,ss PRINT, filein+STRING(dd,mm,yy,hh,mi,ss,FORMAT='(3x,2(I2.2,"/"),I4,I4,2(":",I2.2))') ENDIF to630 = GETHDF(filein=filein,PARAMETER='AOT630_best_model') model = GETHDF(filein=filein,PARAMETER='MODELE') mask = GETHDF(filein=filein,PARAMETER='MASK',/BINARY) Qual = GETHDF(filein=filein,PARAMETER='Quality',/BINARY) mb1 = ( mask MOD 2) EQ 0 mb2 = (ISHIFT(mask,- 1) MOD 2) EQ 0 mb6 = (ISHIFT(mask,- 5) MOD 2) EQ 0 ; Exclure déserts mb8 = (ISHIFT(mask,- 7) MOD 2) EQ 0 mb10= (ISHIFT(mask,- 9) MOD 2) EQ 0 mb11= (ISHIFT(mask,-10) MOD 2) EQ 0 mb13= (ISHIFT(mask,-12) MOD 2) EQ 0 valid = WHERE(mb1 AND mb2 AND mb8 AND mb10 AND mb11 AND mb13 AND (Qual GE 10)) lat = Sev_lat[valid] ; Donner ici les lats et lons SEVIRI lon = Sev_lon[valid] cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels MODIS sont au voisinage de stations aeronet vv = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station IF KEYWORD_SET(verbose) THEN print,'Nombre de pixels proches:',Nc IF Nc EQ 0 THEN RETURN,0 ; Si aucun pixel proche, on passe au fichier suivant ; Continental WMO (alpha = 1.23) ; Moderately Absorbing (1.85) ; Urban Industrial (1.92) ; Smoke (2.11) ; Spherical dust (1.96) ; Maritime WMO (seulement au dessus des océans, alpha = 0.2) angs = [ 1.23, 1.85, 1.92, 2.11, 1.49, 0.2] ; Ang coef for the 6 models to630 = to630[valid[vv]] ang = angs[model[valid[vv]]-1] qual = qual[valid[vv]] cc = cc[vv] ll = ll[vv] lat = lat[vv] lon = lon[vv] sta_indice = sta_indice[vv] Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:1.} Tsat = REPLICATE(Sat,N_ELEMENTS(to630)) Tsat.ista = sta_indice Tsat.lat = lat Tsat.lon = lon Tsat.J = JulObs Lambda0 = 630. Tsat.to500 = to630*(Lambda0/500.)^ang Tsat.to670 = to630*(Lambda0/670.)^ang Tsat.to865 = to630*(Lambda0/865.)^ang Tsat.qual = qual Nv = N_ELEMENTS(TSat) RETURN,TSat END ;============================================================================================================ ; Lecture des produits MODIS ; FILEIN : Chemin vers le fichier hdf ; MASQUE : Taleau masque qui indique le numero de/des station proches (plate carre, resol [2*Ngr,Ngr]) ; MULT : Tableau qui permet de gerer les coincidences multiples ; En sortie, on a "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ; TE si on veut les inversions "terres". TE=0 si on veut les "mer" ;------------------------------------------------------------------------------------------------------------ FUNCTION readMODIS,FILEIN=filein,Nv=Nv,TE=te,VERBOSE=verbose FORWARD_FUNCTION GETHDF COMMON COM_MASQUE,masque,mult Nv=0 s = SIZE(masque) & Ngr=s[2] Get_Location, filein, Jtu=Jtu, DayNightFlag=DNFlag IF KEYWORD_SET(verbose) THEN BEGIN CALDAT,Jtu,mm,dd,yy,hh,mi,ss PRINT, filein+STRING(dd,mm,yy,hh,mi,ss,FORMAT='(3x,2(I2.2,"/"),I4,I4,2(":",I2.2))') IF DNFlag NE 'Day' THEN PRINT,'Night Time Granule' ENDIF IF DNFlag NE 'Day' THEN RETURN,0 lat = GETHDF(filein=filein,PARAMETER='Latitude' ) lon = GETHDF(filein=filein,PARAMETER='Longitude') cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels MODIS sont au voisinage de stations aeronet valid = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station ; print,'Nombre de pixels proches:',Nc IF Nc EQ 0 THEN RETURN,0 ; Si aucun pixel proche, on passe au fichier suivant cc = cc[valid] ll = ll[valid] sta_indice = sta_indice[valid] ; IF TE THEN BEGIN ;Corrected_Optical_Depth_Land Corrected optical thickness at 0.47, 0.55,0.66 micron. Fill -9.99 ;Optical_Depth_Small_Land optical thickness for small Mode at 0.47, 0.55,0.66 & 2.13 [ 135, 203, 4] ;Angstrom_Exponent_Land Angstrom exponent for 0.47 and 0.67 micron ang = GETHDF(filein=filein,PARAMETER='Angstrom_Exponent_Land') & ang = REFORM(ang[valid]) satval = WHERE(ang GT -1.,Npix) IF Npix EQ 0 THEN RETURN,0 ; Si aucun pixel avec des estimations satellite valides ang = ang[satval] pix = valid[satval] od = GETHDF(filein=filein,PARAMETER='Corrected_Optical_Depth_Land') s = SIZE(od ) & od = REFORM(od ,s[1]*LONG(s[2]),s[3]) & od = od [pix,*] odf = GETHDF(filein=filein,PARAMETER='Optical_Depth_Small_Land' ) s = SIZE(odf) & odf = REFORM(odf,s[1]*LONG(s[2]),s[3]) & odf = odf[pix,*] qa = GETHDF(filein=filein,PARAMETER='Quality_Assurance_Land' ) s = SIZE(qa ) & qa = REFORM(qa ,s[1],LONG(s[2])*s[3]) & qa = qa [*,pix] Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Npix) TSat.to500 = od[*,0]*(470./500.)^ang TSat.to670 = od[*,2]*(660./670.)^ang TSat.to865 = od[*,2]*(660./865.)^ang TSat.toF550 =odf[*,1] TSat.qual = REFORM( qa[0,*]/32)/3. ; Vérifié par rapport à la doc. Correspond au qa (0 à 3) pour l'ep optique à 0.66 ENDIF ELSE BEGIN ;Effective_Optical_Depth_Best_Ocean AOT at 7 bands for best solution for 0.47, 0.55,0.66,0.86,1. [ 135, 203, 7] ;Optical_Depth_Small_Best_Ocean AOT at 7 bands for small mode of best solution for 0.47, 0.55, ... od = GETHDF(filein=filein,PARAMETER='Effective_Optical_Depth_Average_Ocean') s = SIZE(od ) & od = REFORM(od ,s[1]*LONG(s[2]),s[3]) & od = od [valid,*] satval = WHERE(od[*,1] GT -1.,Npix) IF Npix EQ 0 THEN RETURN,0 ; Si aucun pixel avec des estimations satellite valides od = od [satval,*] pix = valid[satval] ang = GETHDF(filein=filein,PARAMETER='Angstrom_Exponent_1_Ocean' ) s = SIZE(ang) & ang = REFORM(ang,s[1]*LONG(s[2]),s[3]) & ang = ang[pix,*] odf = GETHDF(filein=filein,PARAMETER='Optical_Depth_Small_Average_Ocean' ) s = SIZE(odf) & odf = REFORM(odf,s[1]*LONG(s[2]),s[3]) & odf = odf[pix,*] qa = GETHDF(filein=filein,PARAMETER='Quality_Assurance_Ocean' ) s = SIZE(qa ) & qa = REFORM(qa ,s[1],LONG(s[2])*s[3]) & qa = qa [*,pix] Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Npix) TSat.to500 = od[*,0]*(470./500.)^ang TSat.to670 = od[*,2]*(660./670.)^ang TSat.to865 = od[*,3]*(860./865.)^ang TSat.toF550 =odf[*,1] TSat.qual = REFORM( qa[0,*]/32)/3. ; Vérifié par rapport a la doc. Correspond au quality for averaged solution ENDELSE tim = GETHDF(filein=filein,PARAMETER='Scan_Start_Time' ) sta_indice = sta_indice[satval] TSat.ista = sta_indice TSat.lat = lat[pix] TSat.lon = lon[pix] TSat.J = tim[pix]/86400. + JULDAY(1,1,1993,0.,0.,0.) ; Nv = N_ELEMENTS(TSat) RETURN,TSat END ;============================================================================================================ ; Lecture des produits MERIS ; FILEIN : Chemin vers le fichier hdf ; En sortie, on Ã "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ; TE si on veut les inversions "terres". TE=0 si on veut les "mer" ; fileTempo permet de copier les fichiers sur un disque local pour ; pouvoir les decompresser ;------------------------------------------------------------------------------------------------------------ FUNCTION readMERIS,FILEIN=filein,Nv=Nv,TE=te,VERBOSE=verbose,fileTempo=fileTempo FORWARD_FUNCTION get_meris_param COMMON COM_MASQUE,masque,mult COMMON MERIS_Param,prod_id,nbcols, nblines, offx, offy Nv=0 s = SIZE(masque) & Ngr=s[2] IF STRPOS(filein,'N1.gz') GT 0 THEN BEGIN cmd = 'cp ' + filein + ' ' + fileTempo + '.gz' ; On copie le fichier sur un disque local pour aller plus vite IF KEYWORD_SET(verbose) THEN PRINT,cmd SPAWN,cmd cmd = 'gunzip -f ' + fileTempo + '.gz' ; On copie le fichier sur un disque local pour aller plus vite IF KEYWORD_SET(verbose) THEN PRINT,cmd SPAWN,cmd ENDIF ELSE fileTempo = filein prod_id = EPR_OPEN_PRODUCT(fileTempo) ;-- Open the product ; Get dimensions of scene nbcols = EPR_GET_SCENE_WIDTH (prod_id) nblines = EPR_GET_SCENE_HEIGHT(prod_id) offx = 0 & offy = 0 ;-- Retrieve acquisition Julian date (begin/end) mph=EPR_GET_MPH(prod_id) beg_date = EPR_GET_FIELD_DATA(EPR_GET_FIELD(mph, 'SENSING_START')) dd=0 & mmstr = '' & yy=0 & hh=0 & mi=0 & ss=0. READS,beg_date, dd, mmstr, yy, hh, mi, ss, format='(I2,1X,a3,1X,I4,1X,I2,1X,I2,1X,F9.6)' ok = WHERE(STRUPCASE(mmstr) EQ ['JAN','FEB','MAR','APR', 'MAY', 'JUN', 'JUL', 'AUG','SEP','OCT','NOV','DEC']) & mm=ok[0]+1 Jbeg=JULDAY(mm,dd,yy,hh,mi,ss) end_date = EPR_GET_FIELD_DATA(EPR_GET_FIELD(mph, 'SENSING_STOP')) READS,end_date, dd, mmstr, yy, hh, mi, ss, format='(I2,1X,a3,1X,I4,1X,I2,1X,I2,1X,F9.6)' ok = WHERE(STRUPCASE(mmstr) EQ ['JAN','FEB','MAR','APR', 'MAY', 'JUN', 'JUL', 'AUG','SEP','OCT','NOV','DEC']) & mm=ok[0]+1 Jend=JULDAY(mm,dd,yy,hh,mi,ss) ; ; Lecture des données ; lat = get_meris_param('latitude') lon = get_meris_param('longitude') flag = get_meris_param('l2_flags') to = get_meris_param('aero_opt_thick') ang = get_meris_param('aero_alpha') EPR_CLOSE_PRODUCT, prod_id ; Les lectures sont terminées. On ferme le fichier Npixels = N_ELEMENTS(flag) IF KEYWORD_SET(te) THEN BEGIN ; On cherche les pixels clairs "valides" clair = WHERE( ISHFT(flag,-23) MOD 2,cc) IF cc EQ 0 THEN RETURN,0 vv = WHERE(to[clair] GE 0.,cc) IF cc EQ 0 THEN RETURN,0 vv = clair[vv] Lambda0 = 443. ENDIF ELSE BEGIN ; On cherche les pixels mer valides clair = WHERE( ISHFT(flag,-21) MOD 2,cc) IF cc EQ 0 THEN RETURN,0 vv = WHERE(to[clair] GE 0.,cc) IF cc EQ 0 THEN RETURN,0 vv = clair[vv] Lambda0 = 865. ENDELSE IF KEYWORD_SET(verbose) THEN BEGIN tempo = MIN(lon,poslon) Jtu = Jbeg + (Jend-Jbeg)*FLOAT(poslon)/Npixels CALDAT,Jtu,mm,dd,yy,hh,mi,ss PRINT, filein+STRING(dd,mm,yy,hh,mi,ss,FORMAT='(3x,2(I2.2,"/"),I4,I4,2(":",I2.2))') ENDIF cc = FIX((lon[vv]+180.)/180.*Ngr) ll = FIX(( 90.-lat[vv])/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels MERIS sont au voisinage de stations aeronet valid = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station ; print,'Nombre de pixels proches:',Nc IF Nc EQ 0 THEN RETURN,0 ; Si aucun pixel proche, on passe au fichier suivant ; lat = lat[vv[valid]] ; On ne garde que les pixels valides lon = lon[vv[valid]] ang = ang[vv[valid]] to = to[vv[valid]] cc = cc[ valid ] ll = ll[ valid ] sta_indice = sta_indice[valid] J = Jbeg + (Jend-Jbeg)*vv[valid]/FLOAT(Npixels) ; Calcul de la date par interpolation Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Nc) TSat.ista = sta_indice TSat.lat = lat TSat.lon = lon TSat.J = J TSat.to500 = to*(Lambda0/500.)^ang TSat.to670 = to*(Lambda0/670.)^ang TSat.to865 = to*(Lambda0/865.)^ang TSat.qual = 1. ; Nv = N_ELEMENTS(TSat) RETURN,TSat END ;============================================================================ ; Lecture d'un parametre d'un fichier MERIS ;============================================================================ FUNCTION get_meris_param,param COMMON MERIS_Param,prod_id,nbcols, nblines, offx, offy band_id = EPR_GET_BAND_ID(prod_id, param) raster = EPR_CREATE_COMPATIBLE_RASTER(band_id, nbcols, nblines, 1,1) status = EPR_READ_BAND_RASTER(band_id, offx, offy, raster) dat = EPR_GET_RASTER_DATA(raster) EPR_FREE_RASTER, raster RETURN,dat END ;============================================================================================================ ; Lecture des produits CALIPSO ; FILEIN : Chemin vers le fichier hdf ; En sortie, on Ã "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ;------------------------------------------------------------------------------------------------------------ FUNCTION readCALIPSO,FILEIN=filein,Nv=Nv,VERBOSE=verbose COMMON COM_MASQUE,masque,mult Nv=0 s = SIZE(masque) & Ngr=s[2] lat = GETHDF(filein=filein,PARAMETER='Latitude' ) & lat = lat[2,*] ; Dimension(3,N); le troisieme est pour la valeur centrale lon = GETHDF(filein=filein,PARAMETER='Longitude') & lon = lon[2,*] ; Dimension(3,N) cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels MERIS sont au voisinage de stations aeronet valid = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station IF KEYWORD_SET(verbose) THEN PRINT,NC,' ',filein IF Nc EQ 0 THEN RETURN,0 cl_freq = GetHDF(FILEIN=filein,PARAM='Surface_Elevation_Detection_Frequency') & cl_freq = cl_freq[ valid] vv = WHERE(cl_freq EQ 163,Nc) ; Cases when surface is seen through the cloud and aerosol layers IF Nc EQ 0 THEN RETURN,0 valid = valid[vv] Time = GetHDF(FILEIN=filein,PARAM='Profile_Time') & Time=REFORM(Time[2,valid]) & J = JULDAY(1,1,1993,0,0,0) + Time/86400. ; Conversion, to Julian ANlay = GetHDF(FILEIN=filein,PARAM='Number_Layers_Found') & ANlay = ANlay [ valid] AOD_532 = GetHDF(FILEIN=filein,PARAM='Feature_Optical_Depth_532') & AOD_532 = AOD_532 [*,valid] AOD_1064= GetHDF(FILEIN=filein,PARAM='Feature_Optical_Depth_1064') & AOD_1064 = AOD_1064[*,valid] to532 = FLTARR(Nc) & FOR i=0,Nc-1 DO IF ANlay[i] NE 0 THEN to532 [i]=TOTAL(AOD_532 [0:ANlay[i]-1,i]) ; Sum the various aerosol layers to1064= FLTARR(Nc) & FOR i=0,Nc-1 DO IF ANlay[i] NE 0 THEN to1064[i]=TOTAL(AOD_1064[0:ANlay[i]-1,i]) ang = ALOG(to532/to1064) / ALOG(2.) Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Nc) TSat.ista = sta_indice[valid] TSat.lat = lat [valid] TSat.lon = lon [valid] TSat.J = J TSat.to500 = to532*(532./500.)^ang TSat.to670 = to532*(532./670.)^ang TSat.to865 = to532*(532./865.)^ang TSat.qual = 1. Nv = Nc RETURN,TSat END ;============================================================================================================ ; Lecture des produits AATSR ; FILEIN : Chemin vers le fichier hdf. Attention. C'est un fichier HDF5 ; En sortie, on Ã "0" si le fichier n'a pas vu de station aeronet, et un tableau de structure ; Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} ; si il y a des stations vues ;------------------------------------------------------------------------------------------------------------ FUNCTION readAATSR,FILEIN=filein,Nv=Nv,VERBOSE=verbose COMMON COM_MASQUE,masque,mult Nv=0 s = SIZE(masque) & Ngr=s[2] file_id = H5F_OPEN(filein) ds_id = H5D_OPEN(file_id,'/dataset1.5/Latitude') lat = H5D_READ(ds_id) / 100. H5D_CLOSE, ds_id ds_id = H5D_OPEN(file_id,'/dataset1.6/Longitude') lon = H5D_READ(ds_id) / 100. H5D_CLOSE, ds_id lon = (lon+540.) MOD 360 - 180. ds_id = H5D_OPEN(file_id,'/dataset1.1/AATSR AOD at 555 nm') aod555 = H5D_READ(ds_id) / 1000. H5D_CLOSE, ds_id ds_id = H5D_OPEN(file_id,'/dataset1.2/AATSR AOD at 659 nm') aod659 = H5D_READ(ds_id) / 1000. H5D_CLOSE, ds_id H5F_CLOSE,file_id cc = FIX((lon+180.)/180.*Ngr) ll = FIX(( 90.-lat)/180.*Ngr) sta_indice = masque[cc,ll] ; On regarde si les pixels MERIS sont au voisinage de stations aeronet valid = WHERE(sta_indice NE 9999,Nc) ; si le masque vaut 9999, le pixel correspondant n'est pas proche d'une station IF KEYWORD_SET(verbose) THEN PRINT,NC,' ',filein IF Nc EQ 0 THEN RETURN,0 pos = STRPOS(filein,'.hdf') date = STRMID(filein,pos-16,16) READS,date,year,mon,day,hr,mi,FORMAT='(I4,4(x,I2))' J = JULDAY(mon,day,year,hr,mi) J=DBLARR(Nc) + J aod555 = aod555[valid] aod659 = aod659[valid] ang = ALOG(aod555/aod659) / ALOG(659./555.) Sat = {ista:0,lat:0.,lon:0.,J:0.d,to500:-1.,to670:-1.,to865:-1.,toF550:-1.,qual:-1.} TSat = REPLICATE(Sat,Nc) TSat.ista = sta_indice[valid] TSat.lat = lat [valid] TSat.lon = lon [valid] TSat.J = J TSat.to500 = aod555*(555./500.)^ang TSat.to670 = aod659*(659./670.)^ang TSat.to865 = aod659*(659./865.)^ang TSat.qual = 1. Nv = Nc RETURN,TSat END ;============================================================================================================ ; Lecture des donnees Aeronet ; INPUT : Chemin vers le fichier ; OUTPUT: station : structure contenant des infos sur la station et les mesures (par des pointeurs) ; VARNAMES : Tableau contenant la liste des variables ; ; DATERANGE = [jday_min,jday_max] Pour specifier une plage de temps utile ; OPTION /JUSTINFO permet de ne lire que le header du fichier (Beaucoup plus rapide) ; /VERBOSE pour faire des impressions de controle ;------------------------------------------------------------------------------------------------------------ FUNCTION read_aeronet,filein,VARNAMES=varnames,JUSTINFO=justinfo,DATERANGE=daterange,VERBOSE=verbose FORWARD_FUNCTION tocoarsefine IF NOT KEYWORD_SET(filein) THEN filein=DIALOG_PICKFILE() tab = FLTARR(2) station = {station,Location:'',long:0.,lat:0.,elev:0,Nmeas:0,PI:'',Email:'',Nobs:0l,Jdays:PTR_NEW(tab),Mes:PTR_NEW(tab)} line='' OPENR,lun,filein,/GET_LUN READF,lun,line ;& print,line READF,lun,line ;& print,line READF,lun,line ;& print,line infos = STRSPLIT(line,",",/EXTRACT) FOR i=0,N_ELEMENTS(infos)-1 DO BEGIN res = STRSPLIT(infos[i],"=",/EXTRACT) IF N_ELEMENTS(res) LT 2 THEN CONTINUE CASE res[0] OF 'Location' : station.Location = res[1] 'long' : station.long = res[1] 'lat' : station.lat = res[1] 'elev' : station.elev = res[1] 'Nmeas' : station.Nmeas = res[1] 'PI' : station.PI = res[1] 'Email' : station.Email = res[1] ENDCASE ENDFOR READF,lun,line READF,lun,line & varnames = STRSPLIT(line,",",/EXTRACT) varnames = varnames[2:*] Nvar = N_ELEMENTS(varnames) choix = WHERE(varnames EQ 'AOT_1020' OR varnames EQ 'AOT_870' OR varnames EQ 'AOT_675' OR $ varnames EQ 'AOT_500' OR varnames EQ 'AOT_440' OR varnames EQ 'AOT_F550',Nwl) varnames = [varnames[choix],'AOT_Coarse','AOT_Fine','Ang_Coarse','Ang_Fine'] IF KEYWORD_SET(justinfo) THEN BEGIN FREE_LUN,lun RETURN,station ENDIF Nmax = 200000l & N=0l res = STRARR(Nvar,Nmax) Jday = DBLARR(Nmax) fmt = '(I2,x,I2,x,I4,x,3(I2,x),A)' WHILE NOT EOF(lun) DO BEGIN READF,lun,day,mon,year,hh,mi,ss,line,FORMAT=fmt jul = JULDAY(mon,day,year,hh,mi,ss) IF KEYWORD_SET(daterange) THEN IF Jul LT daterange[0] OR Jul GT daterange[1] THEN CONTINUE res[*,N] = STRSPLIT(line,",",/EXTRACT) Jday[N] = jul N = N+1 IF N EQ Nmax THEN STOP,'Il faut augmenter Nmax' ENDWHILE FREE_LUN,lun IF KEYWORD_SET(verbose) AND N EQ 0 THEN PRINT,Station.location,FORMAT='(A12," Aucune mesure dans la pÃ©riode temporelle choisie")' station.Nobs = N IF N EQ 0 THEN RETURN,station res = res[*,0:N-1] Jday = Jday[0:N-1] choix = [choix,0,0,0,0] res = res[choix,*] unvalid = WHERE(res EQ 'N/A',cc) IF cc NE 0 THEN res[unvalid] = -9.99 res = FLOAT(res) ;vmax = FLTARR(16) & FOR i=1,16 DO vmax[i-1] = MAX(res[i,*]) ;print,vmax,format='(16F6.2)' ;print,varnames lambda = FLTARR(Nwl) & FOR i=0, Nwl-1 DO lambda[i]=FLOAT(STRMID(varnames[i],4,4)) FOR i=0l, station.Nobs-1 DO res[Nwl:Nwl+3,i]=tocoarseFine(lambda,res[0:Nwl-1,i]) ; ; Lignes ajoutÃ©es pour le cas ou il y a un probleme dans DFPMIN (CATCH, voir ci dessous). ; Si probleme, on a mis le resultat de tocoarseFine a -999 valid = WHERE(res[Nwl,*] GT -990.,cc) IF cc NE station.Nobs THEN BEGIN res = res[*,valid] Jday = Jday[valid] Station.Nobs = cc ENDIF PTR_FREE,station.Jdays,station.Mes Station.Jdays = PTR_NEW(Jday) Station.Mes = PTR_NEW(res) IF KEYWORD_SET(verbose) THEN BEGIN ; Cette ligne et les suivantes pour Ã©crire des informations sur le fichier aeronet Jmin = MIN(*station.Jdays,MAX=Jmax) CALDAT,Jmin,minmon,minday,minyear CALDAT,Jmax,maxmon,maxday,maxyear PRINT,Station.location,minday,minmon,minyear,maxday,maxmon,maxyear,Station.Nobs,$ format='(A12,2(4x,I2,"/",I2.2,"/",I4),2X,I6)' ENDIF RETURN,Station END ;=================================================================================================================== ; PrÃ©paration d'un masque qui permet de repÃ©rer facilement les pixels qui sont "vus" par les stations aeronet ; ;; masque est un tableau de dimension [2*Ngr, Ngr], ou Ngr est par exemple 2000 pour nue grille Ã 10 km de rÃ©solution ; qui indique de quelle station le pixel est proche. Si aucune station, on a 9999. Si plusieurs stations ; on a 7777 ; De plus, le tableau "mult" donne les coordonnees des pixels "multiples" et les stations correspondantes ; mult[0,*] Numero de ligne (commence a zero) ; mult[1,*] Numero de colonne (commence a zero) ; mult[2,*] Nombre de stations (normalement, au mois 2) ; mult[3...,*] Numero des stations qui "voient" se pixel ;=================================================================================================================== PRO build_masque,stations,distmax,VERBOSE=verbose COMMON COM_MASQUE,masque,mult Ngr = 2000 masque = INTARR(2*Ngr,Ngr) + 9999 Nx = 20 mult = 0 & Nmult=0 & NmaxSt=0 lat = 90.-(INDGEN(Ngr)+0.5)*180./Ngr FOR i=0, Ngr-1 DO BEGIN lon = -180. + (INDGEN(2*Ngr)+0.5)*180./Ngr CosLat = COS(lat[i]*!DTOR) FOR j=0, 2*Ngr-1 DO BEGIN dist = SQRT((stations.lat-lat[i])^2 + ((stations.long-lon[j])*CosLat)^2)*40000./360. v = WHERE(dist LT distmax,cc) IF cc EQ 0 THEN CONTINUE IF cc EQ 1 THEN BEGIN masque[j,i] = v[0] CONTINUE ENDIF masque[j,i] = 7777 ; PRINT,cc,stations[v].location pix = INTARR(Nx) & pix[0:2+cc] = [i,j,cc,v] mult = [mult,pix] Nmult = Nmult+1 NmaxSt = NmaxSt > cc ENDFOR ENDFOR mult = REFORM(mult[1:Nx*LONG(Nmult)],Nx,Nmult) mult = mult[0:NmaxSt+2,*] IF KEYWORD_SET(verbose) THEN BEGIN s = INTARR(N_ELEMENTS(stations)) FOR i=0, N_ELEMENTS(stations)-1 DO s[i] = ROUND(TOTAL(masque EQ i)) FOR i=0, N_ELEMENTS(stations)-1 DO print,i,stations[i].location, ROUND(TOTAL(masque EQ i)),format='(I3,A15,I8)' print,777,ROUND(TOTAL(masque EQ 7777)) ENDIF END ;====================================================================================================== ; Routine pour gÃ©rer les pixels satellites "vus" par plusieurs stations ; On duplique les obs et on les affecte aux diffÃ©rentes stations Aeronet ;====================================================================================================== PRO Duplicate_SatObs,TSat,Nv COMMON COM_MASQUE,masque,mult s = SIZE(masque) & Ngr=s[2] im = WHERE(TSat.ista EQ 7777,Nmult) ; Pour gÃ©rer les pixels vus par plusieurs stations IF Nmult EQ 0 THEN RETURN cc = FIX((TSat[im].lon+180.)/180.*Ngr) ll = FIX(( 90.-TSat[im].lat)/180.*Ngr) current = INDGEN(N_ELEMENTS(TSat)) sta_indice = TSat.ista FOR i=0, Nmult-1 DO BEGIN k = WHERE(mult[0,*] EQ ll[i] AND mult[1,*] EQ cc[i],Nverif) & IF Nverif NE 1 THEN STOP, "Probleme" Np = REFORM(mult[2,k]) current = [current,INTARR(Np-1)+im[i]] ; On indique quel est l'indice des tableaux sta_indice[im[i]] = mult[3,k] sta_indice = [sta_indice,mult[4:2+Np,k]] ENDFOR TSat = TSat[current] TSat.ista = sta_indice Nv = N_ELEMENTS(TSat) END ;====================================================================================================== ; On construit la grille de latitudes/longitudes pour SEVIRI ;====================================================================================================== PRO Define_Seviri_latlon COMMON Seviri_latlon,Sev_lat,Sev_lon ; ; Parametres definissant le rayon de la Terre ; req = 6378.140 ; Rayon Ã l'Ã©quateur flat = 1.d0 / 298.25722 ; Applatissement au pole rpo = ( 1.d0 - flat ) * req ; Rayon au pole re2_rp2 = req*req / (rpo*rpo) rfalt = 35785.831 ; Altitude du satellite MSG, en mÃ¨tres ouvradpx = 17.831962 ; Ouverture, en degrÃ©s ouvradln = 17.831962 ; Ouverture, en degrÃ©s nbpix = 3712 nblin = 3712 rsat=req+rfalt angpix = (INDGEN(nbpix)+0.5 -nbpix/2.)/nbpix *ouvradpx*!DTOR anglig = (INDGEN(nblin)+0.5 -nblin/2.)/nblin *ouvradln*!DTOR ; ; cosinus directeurs du vecteur SM satellite-point vise dans repere satellite ;; axes Z, X: plan equatorial, Z : satellite -> centre de la terre ; X : vers l'ouest ; Y: parallele a l'axe des poles, positif vers le nord ; smcx = SIN(angpix) # COS(anglig) smcy = (FLTARR(nbpix)+1.) # SIN(anglig) smcz = COS(angpix) # COS(anglig) ; ; calcul de la distance satellite-point vise ; a = smcx * smcx + re2_rp2 * smcy * smcy + smcz * smcz b = rsat * smcz c = rsat * rsat - req*req delta = b * b - a * c distSM = ( b - SQRT(delta>0.) ) / a ; calcul du vecteur OM, centre de la terre - point vise dans le repere terrestre ; axes Z, X: plan equatorial, Z : centre de la terre -> satellite ; X : vers l'est ; Y: axe des poles, positif vers le nord omx = distSM * smcx omy = - distSM * smcy omz = rsat - distSM * smcz Sev_lat = ATAN ( re2_rp2 * omy / SQRT(omx*omx + omz*omz) ) * !RADEG Sev_lon = ATAN(omx,omz) * !RADEG Sev_lat[WHERE(delta LT 0)] = -999. Sev_lon[WHERE(delta LT 0)] = -999. END ;============================================================================================================ ; Cette fonction retourne une coordonnee lat,lon en une chaine de caracteres "propre" ;------------------------------------------------------------------------------------------------------------ FUNCTION stringPos,lat,lon,alt IF lat GE 0. THEN NS='N' ELSE NS='S' IF lon GE 0. THEN WE='E' ELSE WE='W' RETURN, STRING(ABS(lat),NS,ABS(lon),WE,alt>0.,format='("(",F4.1,A1,",",F5.1,A1,")",I4,"m")') END ;============================================================================================================ ; Procedure de plot avec ajout des statistiques ;------------------------------------------------------------------------------------------------------------ PRO scatter,X,Y,XTIT=xtit,YTIT=ytit,NC=nc,POS=pos,FSAVE=snumb,STN=stn,TE=te FORWARD_FUNCTION densityhisto vmin=0.0001 & vmax=1. & Np=50 IF KEYWORD_SET(te) THEN vsl=0.15 ELSE vsl=0.08 IF KEYWORD_SET(te) THEN vin=0.05 ELSE vin=0.03 IF N_ELEMENTS(X) LE 1 THEN RETURN ; Calcul des statistiques vMinValid = -0.1 v = WHERE(x GT vMinValid AND y GT vMinValid,N) IF N LT 2 THEN RETURN bias = TOTAL (Y[v]-X[v]) /FLOAT(N) RMS = SQRT(TOTAl((Y[v]-X[v])^2)/FLOAT(N)) Corr = CORRELATE (Y[v],X[v]) tempo = WHERE(Y[v] LT (X[v]*(1.+vsl)+vin) AND Y[v] GT (X[v]*(1.-vsl)-vin),Nin) frac = FLOAT(Nin)/FLOAT(N) res = LINFIT(X[v],Y[v]) res[1] = (res[1]>(-9.999))<9.999 IF KEYWORD_SET(snumb) THEN PRINTF,snumb,stn,N,Corr,bias,RMS,res[1],res[0],frac,FORMAT='(A,I6,F7.3,2F6.2,2x,F6.3,F6.2,F6.3)' IF POS[3] EQ 0. THEN RETURN ; C'est le cas dans lequel on veut juste les stats ; ; On fait maintenant le trace ; IF N_ELEMENTS(X) LT 200 THEN $ PLOT,X<(vmax*0.99),Y<(vmax*0.99),PSYM=6,SYMSIZE=0.5,XRANGE=[vmin,vmax],YRANGE=[vmin,vmax],XSTYLE=1,YSTYLE=1,XTIT=xtit,YTIT=ytit,POS=pos,/NOERASE $ ELSE $ PLOTARRAY,DENSITYHISTO(X<(vmax*0.99),Y<(vmax*0.99),XMIN=vmin,XMAX=vmax,YMIN=vmin,YMAX=vmax,NP=Np),$ XRANGE=[vmin,vmax],YRANGE=[vmin,vmax],XTIT=xtit,YTIT=ytit,NC=NC,POS=pos oplot,!X.CRANGE,!X.CRANGE , COL=NC-3 oplot,!X.CRANGE,!X.CRANGE*(1.+vsl)+vin,LINESTYLE=2, COL=NC-3 oplot,!X.CRANGE,!X.CRANGE*(1.-vsl)-vin,LINESTYLE=2, COL=NC-3 XYOUTS,0.02,0.92,STRING(N ,format='("Npts:",I5 )') XYOUTS,0.02,0.85,STRING(Corr,format='("Corr:",F5.2)') XYOUTS,0.98,0.25,STRING(bias,format='("Bias :",F5.2)'),ALIGN=1 XYOUTS,0.98,0.15,STRING(RMS ,format='("RMS :",F5.2)'),ALIGN=1 XYOUTS,0.98,0.05,STRING(ROUND(Frac*100)<99,format='("Valid: ",I2,"%")'),ALIGN=1 END ;============================================================================================================ ; Fonction qui retourne une structure avec quelques informations statistiques ; ; En sortie, on a une structure avec ; - Le nombre de pixel POLDER proches de la station aeronet ; - La valeur Parasol pour le pixel le plus proche ; - La moyenne (ponderee par l'indice de qualite) des valeurs parasol ; - L'ecart type des valeurs parasol ;------------------------------------------------------------------------------------------------------------ FUNCTION TOstat,tab,QUAL=q IF NOT KEYWORD_SET(q) THEN q = FLTARR(N_ELEMENTS(tab))+1. valid = WHERE(tab GT -0.5 AND tab LT 5.,cc) IF cc EQ 0 THEN RETURN,{obs,0,-1.,-1.,-1.} tot = TOTAL(q[valid]) IF tot LT 0.001 THEN RETURN,{obs,cc,tab[valid[0]],-1.,-1.} x = TOTAL(q[valid]*tab[valid] )/FLOAT(tot) x2 = TOTAL(q[valid]*tab[valid]^2)/FLOAT(tot) RETURN,{obs,cc,tab[valid[0]],x,SQRT((x2-x^2)>0.)} END ;======================================================================================================== ; Cette procedure permet de recuperer le jour Julien LOCAL et la longitude de passage a ; l'equateur d'un fichier POLDER ou MODIS ; INPUT : Le chemin vers le fichier (Leader ou Data file pour POLDER/Parasol) ; SORTIE : LONEQ (float) et Jloc (entier long) ;------------------------------------------------------------------------------------------------------------ PRO Get_Location, filename, LONEQ = LonEq, JLOC=Jloc, Jtu=J, DayNightFlag=DayNight FORWARD_FUNCTION get_val IF (STRPOS(filename,'MYD') GE 0 AND STRPOS(filename,'.hdf') GT 0) THEN BEGIN ; On reconnait un fichier MODIS file_handle = HDF_SD_START(filename,/READ) ; SD ID of theHDF file HDF_SD_FILEINFO,file_handle,Ndatasets,Nattribs FOR ia = 0, Nattribs-1 DO BEGIN ; Lecture des attributs pour recuperer des informations sur l'orbite HDF_SD_ATTRINFO,file_handle,ia,COUNT=count, NAME=name,DATA=data,TYPE=type IF name EQ 'CoreMetadata.0' THEN BEGIN EqCrossDate = get_val(data,'EQUATORCROSSINGDATE',/DATETIME) ;& PRINT,EqCrossDate EqCrossTime = get_val(data,'EQUATORCROSSINGTIME',/DATETIME) ;& PRINT,EqCrossTime EqCrossLong = get_val(data,'EQUATORCROSSINGLONGITUDE') ;& PRINT,EqCrossLong BegDate = get_val(data,'RANGEBEGINNINGDATE',/DATETIME) ;& PRINT,BegDate BegTime = get_val(data,'RANGEBEGINNINGTIME',/DATETIME) ;& PRINT,BegTime EndDate = get_val(data,'RANGEENDINGDATE',/DATETIME) ;& PRINT,EndDate EndTime = get_val(data,'RANGEENDINGTIME',/DATETIME) ;& PRINT,EndTime DayNight = get_val(data,'DAYNIGHTFLAG',/DAYNIGHT) ;& PRINT,DayNight ENDIF ENDFOR HDF_SD_END, file_handle ; Fermeture du fichier HDF LonEq = EqCrossLong J = JULDAY(BegDate[1],BegDate[2],BegDate[0],BegTime[0],BegTime[1],BegTime[2]) Jeq = JULDAY(EqCrossDate[1],EqCrossDate[2],EqCrossDate[0],$ EqCrossTime[0],EqCrossTime[1],EqCrossTime[2]) Jloc = ROUND(Jeq + lonEq/360.) ; Calcul du jour "local ENDIF ELSE BEGIN ; C'est un fichier POLDER/Parasol slonEq = BYTARR(8) sdate = BYTARR(16) ll = STRLEN(filename) fileL = STRMID(filename,0,ll-1) + 'L' ; Leaderfile name OPENR,1, fileL POINT_LUN,1,180+360+50 READU,1,slonEq,sdate ; On recupere date et passage a l'equateur CLOSE,1 lonEq = ((FLOAT(STRING(slonEq)) + 180.) MOD 360.) - 180. date = STRING(sdate) READS,date,yy,mm,dd,hh,mi,ss,cc,FORMAT='(I4,6I2)' J = JULDAY(mm,dd,yy,hh,mi,ss+cc/100.) ; Determination jour julien (*.5 a 12h TU) Jloc = ROUND(J + lonEq/360.) ; Calcul du jour "local" ENDELSE END ;=================================================================================================== ; Cette fonction retourne une estimation des epaisseurs optique "Coarse" et "Fine" ; a partir d'un ensemble de mesures de ep optiques pour diverses longueurs d'onde (en nm). ; ; On cherche a fiter les mesures par la somme de deux modes dont les coefficient d'anstrom ; sont proches de 0.2 et 2.5 ;--------------------------------------------------------------------------------------------------- FUNCTION toCoarseFine,lambdaA,toA FORWARD_FUNCTION Func COMMON minim, lambda,to valid = WHERE(toA GE 0. AND toA LT 10.,cc) ; Permet de virer les valeurs non significant IF cc LT 4. THEN RETURN,[-1.,-1,-1.,-1.] s = sort(lambdaA[valid]) ; On classe les mesures par longueur d'onde croissante lambda = lambdaA[valid[s]] to = toA[valid[s]] alphaC = 0.2 ; Premiere estimation d'une epaisseur optique "coarse" et "fine" a partir de valeurs alphaF = 2.5 ; a-priori de leur coefficient d'angstrom to550c = to[cc-1]*(lambda[cc-1]/550.)^alphaC to550f = ((to[0]-to550c*(550./lambda[0])^alphaC)>0.)*(lambda[0]/550.)^alphaF A=[to550c,to550f,alphaC,alphaF] ;measure_errors = 0.02 * to + 0.01 ; Erreurs sur les epaiseurs optiques mesurees ;to_fited = LMFIT(lambda, to, A, MEASURE_ERRORS=measure_error,FUNCTION_NAME = 'tolambda', ; CONVERGENCE=converg,ITER=iter,SIGMA=sigma) errno = 0 ; Procedure CATCH ajoutÃ©e pour contourner les erreurs dans DFPMIN. Si erreur, on se retrouve lÃ , CATCH, errno ; et on peut continuer le traitement IF errno ne 0 then begin CATCH, /CANCEL ; Stop IDL from catching more errors PRINT, !error_state.MSG ; help, !error_state, /STRUCT ; help, /TRACEBACK RETURN,A*0. -999. ENDIF DFPMIN, A, 1.E-5, Fmin, 'Func', 'Dfunc',/DOUBLE,ITMAX=300 ; Ici, on fait une inversion plus precise ;IF cc EQ 4 THEN fmt='(2F6.2,2F6.2,3x,4F6.2,3x,4F6.2,3X,5F6.2)' $ ; ELSE fmt='(2F6.2,2F6.2,3x,5F6.2,3x,5F6.2,3X,5F6.2)' ;print,A,to,FUNC(A,/GETDIF),format=fmt ; On imprime les resultats, les ep optiques mesurÃ©es, les ecarts ;print,dfunc(a) ; au modele, et 5 composantes de la fonction de cout RETURN, A END ;======================================================================================= ; Cette fonction est utilisee pour determiner la repartition Coarse/Fine mode. ; Elle retourne la valeur de la fonction que on cherche a minimiser. ; Avec le keyword GETDIF, elle retourne les differences modele-mesure ;-------------------------------------------------------------------------------------- FUNCTION FUNC,A,GETDIF=getdif COMMON minim, lambda,to lambda0 = 550. RapLamb=lambda0/lambda toC = A[0]*RapLamb^A[2] toF = A[1]*RapLamb^A[3] dif = toC+toF-to err = to*0.02 + 0.02 weight = [TOTAL((dif/err)^2) , exp(((-A[0]/0.001)<50.)>(-50.)) , exp(((-A[1]/0.001)<50.)>(-50.)) , (A[2]-0.2)^2 , (A[3]-2.5)^2] IF KEYWORD_SET(getdif) THEN RETURN, [dif,weight] RETURN,TOTAL(weight) END ;======================================================================================= ; Cette fonction est utilisee pour determiner la repartition Coarse/Fine mode. ; Elle retourne les derivees de la fonction que on cherche a minimiser ;-------------------------------------------------------------------------------------- FUNCTION DFUNC,A COMMON minim, lambda,to lambda0 = 550. RapLamb=lambda0/lambda toC = A[0]*RapLamb^A[2] toF = A[1]*RapLamb^A[3] dif = toC+toF-to err = to*0.02 + 0.02 ;RETURN,TOTAL((dif/err)^2) + exp(-A[0]/0.001) + exp(-A[1]/0.001) + (A[2]-0.5)^2 + (A[3]-2.5)^2 RETURN,[2.*TOTAL(dif/(err^2)*RapLamb^A[2]) - 1000.*exp(((-A[0]/0.001)<50.)>(-50.)), $ 2.*TOTAL(dif/(err^2)*RapLamb^A[3]) - 1000.*exp(((-A[1]/0.001)<50.)>(-50.)), $ 2.*TOTAL(dif/(err^2)*toC*ALOG(RapLamb)) + 2*(A[2]-0.2) , $ 2.*TOTAL(dif/(err^2)*toF*ALOG(RapLamb)) + 2*(A[3]-2.5) ] END ;=========================================================================== FUNCTION get_val,string, name,DATETIME=datetime,DAYNIGHT=daynight ;=========================================================================== posB = STRPOS(string,name) posE = STRPOS(string,name,/REVERSE_SEARCH) substring = STRMID(string,posB,posE-posB-1) val = STRSPLIT(substring,' ',/EXTRACT) ip = WHERE(val EQ 'VALUE') ;FOR i=0, N_ELEMENTS(val)-1 DO print,i,' ',val[i] val = val[ip+2] IF KEYWORD_SET(daynight) THEN RETURN, STRMID(val,1,3) IF KEYWORD_SET(datetime) THEN RETURN, FLOAT((STRSPLIT(val,'"-:',/EXTRACT))[0:2]) ELSE RETURN, FLOAT(val[0]) END ;=========================================================================== ; Retourne une chaine de caractÃ¨re pour indiquer un temps en minute, heure ou jours ;=========================================================================== FUNCTION STMIN,deltaT dT = ROUND(deltaT) IF dT LT 180 THEN RETURN,'Mesure la plus proche: '+STRING(dT,format='(I3)')+' minutes' IF dT LT 24*60 THEN RETURN,'Mesure la plus proche: '+STRING(ROUND(dT/60.),format='(I3)')+' heures' RETURN,'Mesure la plus proche: '+STRING(ROUND(dT/60./24.),format='(I3)')+' jours' END ;========================================================================================================== ; gethdf.pro ;========================================================================================================== ; AUTEUR : F.M. Breon, Septembre 2006 ;========================================================================================================== ; Cette fonction permet de recuperer facilement un Scientific dataset d'un fichier HDF. ; Plusieurs autres fonctions d'information si les parametres ne sont pas donnes ; ; USAGE: resu = gethdf(FILEIN=,PARAMETER=,/INFO,/BINARY_VALUE) ; Tous les parametres sont optionnels ; ; FILEIN=filein : Le nom du fichier HDF a lire. Si absent, un dialogue apparait ; PARAMETER = param : Nom court du parametre a recuperer. ; /INFO : Pour obtenir des informations supplementaires (attributs) sur le parametre en question ; /BINARY_VALUE : Pour retourner les valeurs binaires (par defaut, valeurs physiques) ; ; EXEMPLES ; liste_param = getHDF() Ouvre un dialogue pour choisir le fichier, donne la liste des parametres ; accessibles, Retourne la liste des parametres accessibles ; liste_param = getHDF(FILEIN=filein) Retourne la liste des parametres accessibles ; ; data = getHDF(FILEIN=filein,PARAMETER=param). Retourne le parametre "param" en valeurs physiques ; ; data = getHDF(FILEIN=filein,PARAMETER=param,/INFO). Retourne le parametre "param" et imprime ses attributs ; ; data = getHDF(FILEIN=filein,PARAMETER=param,/BINARY). Retourne le parametre "param" en valeurs de codage ; ; Lorsque le fichier n'est pas trouve, la fonction retourne -1 ; Lorsque le parametre choisi n'est pas disponible, la fonction retourne -2 ; FUNCTION getHDF,FILEIN=filein,PARAMETER=param,INFO=info,BINARY_VALUE=binary IF NOT KEYWORD_SET(filein) THEN filein=DIALOG_PICKFILE(FILTER='*.hdf') file_handle = HDF_SD_START(filein,/READ) IF file_handle LE 0 THEN BEGIN print,' File not found' RETURN,-1 ENDIF IF NOT KEYWORD_SET(param) THEN BEGIN PRINT,'Available parameters:' HDF_SD_FILEINFO,file_handle,Nsds,Nattribs liste = STRARR(Nsds) FOR id=0, Nsds-1 DO BEGIN sd_id = HDF_SD_SELECT(file_handle,id) HDF_SD_GETINFO,sd_id,NAME=name liste[id] = name PRINT,name ENDFOR PRINT,'More info with the command: HDFdump,FILEIN=filein,/NOSTOP' HDF_SD_END, file_handle RETURN,liste ENDIF id = HDF_SD_NAMETOINDEX(file_handle,param) IF id LT 0 THEN BEGIN PRINT,'This parameter is not available. Available parameters:' HDF_SD_FILEINFO,file_handle,Nsds,Nattribs FOR id=0, Nsds-1 DO BEGIN sd_id = HDF_SD_SELECT(file_handle,id) HDF_SD_GETINFO,sd_id,NAME=name PRINT,name ENDFOR HDF_SD_END, file_handle RETURN,-2 ENDIF sd_id = HDF_SD_SELECT(file_handle,id) HDF_SD_GETDATA,sd_id,data IF KEYWORD_SET(info) THEN BEGIN HDF_SD_GETINFO,sd_id,DIMS=dims,NAME=sdsname,NDIMS=ndims,TYPE=type,UNIT=unit, $ NATTS=Natts,COORDSYS=coordsys,LABEL=label,RANGE=range,FILL=fill;,CALDATA=caldata CASE Ndims OF 1 : print,id,sdsname,type,dims,unit,format='(I3,x,A30,A8," [",I4,"]", 11x,A10)' 2 : print,id,sdsname,type,dims,unit,format='(I3,x,A30,A8," [",I4,",",I4,"]", 6x,A10)' 3 : print,id,sdsname,type,dims,unit,format='(I3,x,A30,A8," [",2(I4,","),I4,"]", x,A10)' 4 : print,id,sdsname,type,dims,unit,format='(I3,x,A30,A8," [",3(I4,","),I4,"]", x,A10)' ENDCASE ; print,name,type,unit,ndims,dims,format='(/,A25,A8,A10,I3,5I5)' print,"--------------------------------------------------------------------------------------------------" FOR ia=0, natts-1 DO BEGIN HDF_SD_ATTRINFO,sd_id,ia,COUNT=count, NAME=name,DATA=adata,TYPE=type CASE type OF 'STRING' : PRINT,ia,name,adata,format='(10x,I2,2x,A20,x,A60)' 'FLOAT' : PRINT,ia,name,count,adata[0:(count-1)<10],format='(10x,I2,2x,A20,2x,I4,6E12.3)' 'DOUBLE' : PRINT,ia,name,count,adata[0:(count-1)<10],format='(10x,I2,2x,A20,2x,I4,6E12.3)' 'UINT' : PRINT,ia,name,count,adata[0:(count-1)<10],format='(10x,I2,2x,A20,2x,I4,11I8)' 'INT' : PRINT,ia,name,count,adata[0:(count-1)<10],format='(10x,I2,2x,A20,2x,I4,11I8)' 'BYTE' : PRINT,ia,name,count,adata[0:(count-1)<19],format='(10x,I2,2x,A20,2x,I4,20I4)' 'LONG' : PRINT,ia,name,count,adata[0:(count-1)<19],format='(10x,I2,2x,A20,2x,I4,20I8)' 'ULONG' : PRINT,ia,name,count,adata[0:(count-1)<19],format='(10x,I2,2x,A20,2x,I4,20I8)' ELSE: PRINT,type,name,count,adata ENDCASE ENDFOR ENDIF ; ; If the type of data is integer, it may need to be scaled ; v = WHERE([2,3,12,113,14,15] EQ SIZE(data,/TYPE),cc) IF NOT KEYWORD_SET(binary) AND cc NE 0 THEN BEGIN HDF_SD_GETINFO,sd_id,NATTS=Natts scale = -999. & offset=-999. FOR ia=0, natts-1 DO BEGIN HDF_SD_ATTRINFO,sd_id,ia,NAME=name,DATA=cdata IF name EQ 'scale_factor' THEN scale = cdata[0] IF name EQ 'add_offset' THEN offset = cdata[0] ENDFOR IF scale EQ -999. OR offset EQ -999. THEN PRINT,'Unable to get calibration Information' $ ELSE data = FLOAT((data-offset)*scale) ENDIF HDF_SD_END, file_handle RETURN,data END ;=================================================================================================================== ; Cette routine permet de tracer ou on veut un tableau bidim. C'est donc une extension de la ; commande tv qui plote un tableau dont la taille est determinee par la dimension du tableau ; En plus, plotarray permet de mettre des axes (sauf si "nobord" est active). ; En entree, le tableau tab. Si tab est de type "byte", on le trace tel quel. Sinon, on applique ; un scaling entre les valeurs min (code couleur 1) et max (code couleur ncol-2). ;=================================================================================================================== PRO plotarray,tab,POS=pos,XTIT=xtit,XRANGE=xrange,YTIT=ytit,YRANGE=yrange,NOBORD=nobord,NC=Nc,TIT=tit,BADDATA=baddata FORWARD_FUNCTION get_range s=size(tab) IF (s[0] NE 2) THEN RETURN IF NOT KEYWORD_SET(NC) THEN NC=256 IF s[3] EQ 1 THEN visu=tab ELSE BEGIN res=get_range(tab, PERCENT=1, UNVALID=baddata) vmin = res[0] & vmax = res[1] ; vmin=MIN(tab,MAX=vmax) visu=BYTE(1+(tab-vmin)*(NC-3.)/(vmax-vmin) ) IF KEYWORD_SET(baddata) THEN visu(where(tab EQ baddata)) = !P.BACKGROUND ENDELSE IF NOT KEYWORD_SET(pos) THEN pos=[0.05,0.05,0.95,0.95] IF NOT KEYWORD_SET(xrange) THEN xrange=[0,s(1)] IF NOT KEYWORD_SET(yrange) THEN yrange=[0,s(2)] IF !D.NAME EQ 'PS' THEN BEGIN tv,visu,pos(0),pos(1),/NORMAL,xsize=pos(2)-pos(0),ysize=pos(3)-pos(1) ENDIF ELSE BEGIN pxsize = LONG((POS(2)-POS(0))*!D.X_SIZE + 0.5) pysize = LONG((POS(3)-POS(1))*!D.Y_SIZE + 0.5) V = BYTARR(pxsize,pysize) Px = FIX(INDGEN(pxsize)*s(1)/float(pxsize)) Py = FIX(INDGEN(pysize)*s(2)/float(pysize)) FOR i = 0, Pysize-1 DO V(*,i)=visu(Px,Py(i)) tv,V,pos(0),pos(1),/NORMAL ENDELSE IF NOT KEYWORD_SET(nobord) THEN $ plot,xrange,yrange,XRANGE=xrange,YRANGE=yrange,XSTYLE=1,YSTYLE=1,$ XTIT=xtit,YTIT=ytit,/NODATA,/NOERASE,/NORMAL,POS=pos,TIT=tit END ;======================================================================== ; NAME: ; densityhisto ; ; PURPOSE: ; Builds a density (2D) histogram of two vectors X and Y ; Most usefull to analyze the relationship between two variables ; ; KEYWORDS: ; ; ; CALLING SEQUENCE: ; result = densityhisto(X,Y,XMIN=xmin,XMAX=xmax,YMIN=ymin,YMAX=ymax,NP=npts) ; Can be used with fmb_plotarray ; ; INPUTS: ; X and Y, two vectors of same size ; ; OPTIONAL INPUT PARAMETERS: ; XMIN, XMAX, YMIN, YMAX : To set the range of "usefull" values for both X and Y ; NP : Number of bins in each dimension. The default is 100 ; NC : Value to be used in the byte array for the bin with the largest number of points. ; Default is 247 (saves "reserved" color values) ; ; OUTPUTS: ; A 2D arrays of bytes. Can be plotted with fmb_plotarray ; ; RESTRICTIONS: ; ; ;======================================================================== FUNCTION densityhisto,X,Y,XMIN=xmin,XMAX=xmax,YMIN=ymin,YMAX=ymax,NP=npts,NC=nc IF N_ELEMENTS(xmin) EQ 0 THEN xmin = MIN(X) IF N_ELEMENTS(xmax) EQ 0 THEN xmax = MAX(X) IF N_ELEMENTS(ymin) EQ 0 THEN ymin = MIN(Y) IF N_ELEMENTS(ymax) EQ 0 THEN ymax = MAX(Y) valid =WHERE(X GE xmin AND X LT Xmax AND Y GE Ymin AND Y LT Ymax,cc) IF cc EQ 0 THEN RETURN,0 IF NOT KEYWORD_SET(Npts) THEN Npts = 100 IF NOT KEYWORD_SET(Nc) THEN Nc = 247 ix = FIX((X[valid]-Xmin)*Npts/FLOAT(Xmax-Xmin)) < (Npts-1) iy = FIX((Y[valid]-Ymin)*Npts/FLOAT(Ymax-Ymin)) < (Npts-1) histoB = INTARR(Npts,Npts) FOR i=0l, cc-1 DO histoB[ix[i],iy[i]] = histoB[ix[i],iy[i]] + 1 histoB = BYTE(histoB*FLOAT(Nc)/MAX(histoB)) RETURN,histoB END ;======================================================================== PRO wpng,filename DEVICE,/DECOMPOSED WRITE_PNG,filename,TVRD(/TRUE) DEVICE,DECOMPOSED=0 END ;======================================================================== PRO mypalette,PSCRIPT=ps,NC=n loadct,13 tvlct,r,g,b,/get N = N_ELEMENTS(r) r[ 0]=255 & g[ 0]=255 & b[ 0]=255 r[N-1]= 0 & g[N-1]= 0 & b[N-1]= 0 tvlct,r,g,b !P.BACKGROUND = 0 !P.COLOR = N-1 END ;================================================================================= ; get_range : Retourne un tableau [vmin,vmax] ou vmin et vmax sont les valeurs ; telles que n % du tableau ait des valeurs inferieures a vmin, et ; n % ait des valeurs superieures a vmax. ; INPUT : data. Un tableau de dimension et type quelconque ; percent : Le pourcentage choisi (par defaut, 5) ; bad: : Permet de definit une valeur a ne pas prendre en compte ;================================================================================= FUNCTION get_range, data, PERCENT=percent, UNVALID=bad IF N_ELEMENTS(bad) EQ 0 THEN dat=data ELSE dat= data[WHERE(data NE bad)] N = N_ELEMENTS(dat) IF N_ELEMENTS(percent) EQ 0 THEN percent=5. IF percent LT 0. OR percent GT 40. THEN BEGIN PRINT, 'percent must be in range [0, 40]' RETURN,[0.,0.] ENDIF IF N EQ 0 THEN BEGIN print,' No valid element in data' RETURN,[0,0] ENDIF vmin = MIN(dat,MAX=vmax) ; ; On fabrique l'histogramme cumule ; nbins = 5000 histo = HISTOGRAM(dat,MIN=vmin,MAX=vmax,NBINS=nbins) chisto = LONG(histo) FOR i=1l, N_ELEMENTS(histo)-1 DO chisto[i] = chisto[i-1] + histo[i] ; ; recherche minimum ; vm = percent/100. * N imin = 0l & WHILE (chisto[imin] LT vm) DO imin=imin+1 vm = (1.-percent/100.) * N imax = 0l & WHILE (chisto[imax] LT vm) DO imax=imax+1 RETURN,vmin + [imin,imax]*(vmax-vmin)/FLOAT(nbins) END