;=================================================================================================== ; ; VISU_GLAS Version 1.0.0 ; ;=================================================================================================== ; ; Copyright ICARE ; Auteur : SIX Bruno, Juin 2006 ; Bruno.Six@univ-lille1.fr ; ; Ce logiciel est régi par la licence CeCILL soumise au droit français et ; respectant les principes de diffusion des logiciels libres. Vous pouvez ; utiliser, modifier et/ou redistribuer ce programme sous les conditions de la ; licence CeCILL telle que diffusée par le CEA, le CNRS et l'INRIA sur le site ; "http://www.cecill.info". ; ; En contrepartie de l'accessibilité au code source et des droits de copie, de ; modification et de redistribution accordés par cette licence, il n'est offert ; aux utilisateurs qu'une garantie limitée. Pour les mêmes raisons, seule une ; responsabilité restreinte pèse sur l'auteur du programme, le titulaire des ; droits patrimoniaux et les concédants successifs. ; ; A cet égard, l'attention de l'utilisateur est attirée sur les risques associés ; au chargement, à l'utilisation, à la modification et/ou au développement et à ; la reproduction du logiciel par l'utilisateur étant donné sa spécificité de ; logiciel libre, qui peut le rendre complexe à manipuler et qui le réserve donc ; à des développeurs et des professionnels avertis possédant des connaissances ; informatiques approfondies. Les utilisateurs sont donc invités à charger et ; tester l'adéquation du logiciel à leurs besoins dans des conditions permettant ; d'assurer la sécurité de leurs systèmes et ou de leurs données et, plus ; généralement, à l'utiliser et l'exploiter dans les mêmes conditions de ; sécurité. ; ; Le fait que vous puissiez accéder à cet en-tête signifie que vous avez pris ; connaissance de la licence CeCILL, et que vous en avez accepté les termes. ; ;=================================================================================================== PRO visu_glas FORWARD_FUNCTION Create_GDI, fileselect pstate=Create_GDI() par = (*pstate).wstrparam (*par).wxsize = 1000 & (*par).wysize = 280 (*par).SegRec2 = 2 & (*par).SegRec = 2*(*par).SegRec2 + 1 (*par).SegLthL2 = 50 & (*par).SegLthL1 = 4*(*par).SegLthL2 (*par).torig = JULDAY(1, 1, 2000, 12., 0., 0.) & (*par).LidName = 'GLAS' (*par).altitudes = '10.5' ; IF fileselect(pstate) THEN Display, pstate, /MAPFIRST END PRO setWindow, pstate, xw, yw, typ par = (*pstate).wStrParam & wid = (*pstate).wstrWidg IF (*par).numwin GE 0 THEN BEGIN DEVICE, WINDOW_STATE=winstate IF winstate[(*par).numwin] AND (*par).wxsize EQ xw AND (*par).wysize EQ yw AND typ EQ (*par).ftype THEN RETURN WDELETE, (*par).numwin END DEVICE, WINDOW_STATE=winstate FOR iwin = 1, N_ELEMENTS(winstate) DO IF winstate[iwin] EQ 0 THEN BREAK (*par).numwin = iwin & (*par).l1win0 = iwin + 7 & (*par).l1win = (*par).l1win0 WINDOW, (*par).numwin, xsize=xw, ysize=yw, TITLE='Visualisation '+typ (*par).wxsize = xw & (*par).wysize = yw & (*par).ftype = typ END ;====================================================================================== ; Initialisation paramètres et GUI en fonction du type de produit GLAS ;------------------------------------------------------------------------------------- PRO setTypeParms, pstate, GlasType par = (*pstate).wStrParam & wid = (*pstate).wstrWidg IF GlasType EQ (*par).ftype THEN RETURN (*par).type = GlasType & (*par).Nsamp = 1 DEVICE, WINDOW_STATE=winstate SWITCH GlasType OF 'GLA08' : 'GLA09' : 'GLA10' : 'GLA11' : $ BEGIN (*par).ilat = 27 & (*par).ilon = (*par).ilat + 4 (*par).LidNiv = 2 & (*par).SegLength = (*par).SegLthL2 BREAK; END ELSE : $ BEGIN IF GlasType EQ 'GLA02' THEN (*par).ilat = 3 ELSE (*par).ilat = 9 (*par).ilon = (*par).ilat + 1 (*par).LidNiv = 1 & (*par).SegLength = (*par).SegLthL1 & (*par).L1type = '' IF (*par).l1win GE 0 THEN IF winstate[(*par).l1win] THEN WDELETE, (*par).l1win END ENDSWITCH IF (*wid).wVisuBase LT 0 THEN (*wid).wVisuBase = WIDGET_BASE((*wid).wParamBase, /COLUMN, FRAME=2) $ ELSE Display, pstate, /RESET gid = WIDGET_INFO((*wid).wVisuBase, FIND_BY_UNAME='GlasBase') IF gid GT 0 THEN WIDGET_CONTROL, gid, /DESTROY GlasBase = WIDGET_BASE((*wid).wVisuBase, /COLUMN, UNAME='GlasBase') gid = WIDGET_INFO((*wid).wParamBase, FIND_BY_UNAME='Lev1Base') llev = ['Surf', 'Altitudes'] SWITCH GlasType OF 'GLA02' : 'GLA07' : $ BEGIN if gid GT 0 THEN WIDGET_CONTROL, gid, /DESTROY setMinMax, *par, GlasType MnMxBase = WIDGET_BASE(GlasBase, /ROW, UNAME='MnMxBase', YPAD=0) MinField = CW_FIELD(MnMxBase, VALUE = (*par).minv, /FLOATING, UVALUE = 'min', XSIZE = 8, TIT = 'MIN', /RETURN_EVENTS, /COLUMN) MaxField = CW_FIELD(MnMxBase, VALUE = (*par).maxv, /FLOATING, UVALUE = 'max', XSIZE = 8, TITLE='MAX', /RETURN_EVENTS, /COLUMN) (*wid).wMinF = MinField & (*wid).wMaxF = MaxField; & (*wid).wRelf = -1l lev = ['532 nm', '1064 nm'] Lev1Group=CW_BGROUP(GlasBase, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UNAME='Channel', UVALUE='chann', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) BREAK; END ELSE : $ BEGIN IF GlasType EQ 'GLA09' THEN BEGIN lev = ['4s sampling', '1s sampling'] SampGroup=CW_BGROUP(GlasBase, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UVALUE='sampl', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) ENDIF IF GlasType EQ 'GLA10' OR GlasType EQ 'GLA11' THEN BEGIN llev = ['PBL', 'Altitudes'] setMinMax, *par, GlasType MnMxBase = WIDGET_BASE(GlasBase, /ROW, UNAME='MnMxBase', YPAD=0) MinField = CW_FIELD(MnMxBase, VALUE = (*par).minv, /FLOATING, UVALUE = 'min', XSIZE = 8, TIT = 'MIN', /RETURN_EVENTS, /COLUMN) MaxField = CW_FIELD(MnMxBase, VALUE = (*par).maxv, /FLOATING, UVALUE = 'max', XSIZE = 8, TITLE='MAX', /RETURN_EVENTS, /COLUMN) (*wid).wMinF = MinField & (*wid).wMaxF = MaxField; & (*wid).wRelf = -1l lev = ['Clouds', 'Aerosols'] SampGroup=CW_BGROUP(GlasBase, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UNAME='Aerosol', UVALUE='aero', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) lev = ['Both', 'Profile', 'Layers'] Lev2Group=CW_BGROUP(GlasBase, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UNAME='TypDisp', UVALUE='tdisp', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) ENDIF StyBase = WIDGET_BASE(GlasBase, /ROW, UNAME='StyBase', YPAD=0) LabSty = WIDGET_LABEL(StyBase, UNAME='LabSty', VALUE="Style : ", UVALUE='labsty') Style = Widget_Droplist(StyBase, UNAME='Style', VALUE=['Rectangles', 'Contour'], UVALUE='style') if gid GT 0 THEN BREAK Lev1Base = WIDGET_BASE((*wid).wParamBase, /COLUMN, UNAME='Lev1Base', YPAD=0, FRAME=2) lev = ['Level 1 related segment'] RelFil = CW_BGROUP(Lev1Base, lev, BUTTON_UVALUE=lev, SET_VALUE=[0], UVALUE='relf', /NONEXCLUSIVE, /ROW) Lev1Base2 = WIDGET_BASE(Lev1Base, /COLUMN, MAP=0, UNAME='Lev1Base2', SPACE=0, XPAD=0, YPAD=0) lev = ['GLA02', 'GLA07'] Lev1Group=CW_BGROUP(Lev1Base2, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UNAME='Lev1RelSeg', UVALUE='lev1', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) ;--------------------- setMinMax, *par, 'GLA02' MnMxBase = WIDGET_BASE(Lev1Base2, /ROW, UNAME='MnMxBase', YPAD=0) MinField = CW_FIELD(MnMxBase, VALUE = (*par).minv, /FLOATING, UVALUE = 'min', XSIZE = 8, TIT = 'MIN', /RETURN_EVENTS, /COLUMN) MaxField = CW_FIELD(MnMxBase, VALUE = (*par).maxv, /FLOATING, UVALUE = 'max', XSIZE = 8, TITLE='MAX', /RETURN_EVENTS, /COLUMN) (*wid).wMinF = MinField & (*wid).wMaxF = MaxField; & (*wid).wRelf = -1l ;--------------------- lev = ['532 nm', '1064 nm'] Lev1Group2=CW_BGROUP(Lev1Base2, lev, BUTTON_UVALUE=lev, SET_VALUE=0, UNAME='Channel', UVALUE='chann', /EXCLUSIVE, /ROW, /NO_RELEASE, SPACE=0, YPAD=0) Lev1Base3 = WIDGET_BASE(Lev1Base2, /ROW, UNAME='Lev1Base3', SPACE=0, XPAD=0, YPAD=0) lev = ['Superposition'] RelSuper = CW_BGROUP(Lev1Base3, lev, BUTTON_UVALUE=lev, SET_VALUE=[0], UVALUE='super', /NONEXCLUSIVE, /ROW, SPACE=0, YPAD=0) (*wid).wRefFile = WIDGET_LABEL(Lev1Base2, VALUE=' ', UNAME='Lev1File', /DYNAMIC_RESIZE) (*wid).wsuper = RelSuper & END ENDSWITCH LineGroup=CW_BGROUP(GlasBase, llev, BUTTON_UVALUE=llev, SET_VALUE=[1b, 1b], UVALUE='lines', /NONEXCLUSIVE, /ROW, SPACE=0, YPAD=0) (*wid).wLines = LineGroup END ;================================================================================================ ; Initialisation des valeurs min et max des profils en fonction du type de produit GLAS niveau 1 ;------------------------------------------------------------------------------------------------ PRO SetMinMax, p, typ p.minv = 7. IF typ EQ 'GLA02' THEN p.maxv = 13. ELSE p.maxv = 21. END ;======================================================================== ; Param Event : Various actions depending on events ;------------------------------------------------------------------------ PRO Param_Event, event FORWARD_FUNCTION ResetPlot, fileselect WIDGET_CONTROL, event.top, GET_UVALUE=pstate WIDGET_CONTROL, event.id , GET_UVALUE=uval WIDGET_CONTROL, event.id , GET_VALUE=val par = (*pstate).wstrParam & wid = *(*pstate).wstrWidg DEVICE, WINDOW_STATE=winstate CASE uval OF 'relf' : $ BEGIN WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='Lev1Base2'), MAP=val WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='Lev1File' ), SET_VALUE=' ' IF val EQ 1 THEN BEGIN WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='Lev1RelSeg'), GET_VALUE=val2 CreateL1Window, pstate, val2 ENDIF ELSE BEGIN (*par).L1type = '' WDELETE, (*par).l1win ENDELSE END 'lev1' : CreateL1Window, pstate, val 'sampl': $ BEGIN IF val EQ 0 THEN (*par).Nsamp = 1 ELSE (*par).Nsamp = 4 Display, pstate, MAPFIRST=-1 END 'tdisp': $ BEGIN (*par).typdisp = val & Display, pstate, MAPFIRST=-1 END 'aero': $ BEGIN (*par).aerosol = val DEVICE, WINDOW_STATE=winstate IF winstate[(*par).l1win] THEN WDELETE, (*par).l1win Display, pstate, MAPFIRST=-1 END 'super': $ BEGIN DisplayProduct, pstate & DisplayL1RelatedProfile, pstate END ELSE: $ CommonParamEvent, event ENDCASE END ;=================================================================================== ; Lecture de la partie d'un profil de niveau 1 limitée à un fichier. ;----------------------------------------------------------------------------------- PRO readSingleProfl1, pstate, irec0, toread, Profile, Line, LinLeg, nrec, type, lun, recsize, high par = (*pstate).wstrparam Gla02PixOffset=2 Nshot = (*par).SegLthL1 yoffset = 1500. Nlow = 148 & Nmed = 132 & Nhigh0 = 268 Nsamp = 5 IF high THEN BEGIN yhight0 = 41000 & Nhigh = Nhigh0 ENDIF ELSE BEGIN yhight0 = 20500 & Nhigh = 0 ENDELSE yhight = yhight0 + yoffset Ntot = Nlow + Nmed + Nhigh yscale = Ntot/yhight yhgmax = (*par).channel EQ 0 ? 41000 : 20500 Nmax = Nlow + Nmed + ((*par).channel EQ 0 ? Nhigh0 : 0) ON_IOERROR, EOF CASE type OF 'GLA02' : BEGIN Profile2 = LONARR(Nmax, Nsamp) resamp = 8 dem_min = 0 & dem_max = 0 i40_lid = LONARR(Nlow, resamp*Nsamp) i5_lid = LONARR(Nmed, Nsamp) IF high THEN i1_lid = LONARR(Nhigh) offset = (*par).channel EQ 0 ? 36L : 30224L samp = REFORM(INDGEN(resamp*Nsamp), resamp, Nsamp) offset2 = (*par).channel EQ 0 ? 28476L : 29500L i_rng = 0L & i_Hsat = 0L & bin_hgt = 0L & bin_low = 10000000L FOR irec=0, toread-1 DO BEGIN rrec = irec + Nshot - toread point_lun, lun, recsize*(irec0 + irec) + 20L READU, lun, dem_min, dem_max point_lun, lun, recsize*(irec0 + irec) + offset READU, lun, i40_lid, i5_lid, i1_lid point_lun, lun, recsize*(irec0 + irec) + offset2 READU, lun, i_rng point_lun, lun, recsize*(irec0 + irec) + 56636L READU, lun, i_Hsat Gla02PixOffset = ROUND((yhgmax - 1.e-2*(i_Hsat - i_rng))*yscale) ddd = (i_Hsat - i_rng)/100 dde = i_rng/100 bin_hgt = bin_hgt > (i_Hsat - i_rng) ; pour la correction du décalage entre GLA02 et GLA07 bin_low = bin_low < (i_Hsat - i_rng) ; lié à la prise en compte de l'altitude du satellite ; Pour le niveau bas, on fait la moyenne des mesures FOR i=0, Nsamp-1 DO Profile2[Nmax-Nlow:Nmax-1, i] = TOTAL(i40_lid[*, samp[*, i]], 2)/resamp ; pour le niveau moyen, on prend la mesure directement Profile2[Nmax-Nlow-Nmed:Nmax-Nlow-1, *] = i5_lid ; pour le niveau haut, si présent, on rééchantillonne IF high THEN FOR i=0, Nsamp-1 DO Profile2[0:Nhigh0-1, i] = i1_lid decal = ABS(Gla02PixOffset) + 2 IF Ntot EQ Nmax THEN BEGIN IF Gla02PixOffset LT 0 $ THEN Profile[0:Ntot - decal - 1, *, rrec] = Profile2[Nmax - Ntot + decal:Nmax - 1, *] $ ELSE Profile[decal:Ntot - 1, *, rrec] = Profile2[Nmax - Ntot:Nmax - decal - 1, *] ENDIF ELSE BEGIN IF Gla02PixOffset LT 0 $ THEN Profile[0:Ntot - decal - 1, *, rrec] = Profile2[Nmax - Ntot + decal:Nmax - 1, *] $ ELSE Profile[0:Ntot - 1, *, rrec] = Profile2[Nmax - Ntot - decal:Nmax - decal - 1, *] ENDELSE Line[0, *, rrec] = FIX(ROUND((dem_Min + yoffset)*yscale)) Line[1, *, rrec] = FIX(ROUND((dem_Max + yoffset)*yscale)) LinLeg = [['DEM min', '-3'], ['DEM max', '3']] ENDFOR END 'GLA07' : BEGIN Profoffset = 250. ; Pour GLA07 : range -1.75km à 22.25km au lieu de -1.5km à 22.5km Gla07PixOffset = ROUND(yscale*Profoffset) elev = 0l & i_atm_dem = 0l offset = (*par).channel EQ 0 ? 1952 : 36592 Nbscs = Nlow + Nmed + ((*par).channel EQ 0 ? 268 : 0) i5_bscs = LONARR(Nbscs, Nsamp) FOR irec=0, toread-1 DO BEGIN rrec = irec + Nshot - toread point_lun, lun, recsize*(irec0 + irec) + 72 READU, lun, elev point_lun, lun, recsize*(irec0 + irec) + 1920 READU, lun, i_atm_dem point_lun, lun, recsize*(irec0 + irec) + offset READU, lun, i5_bscs Profile[0:Ntot-Gla07PixOffset-1, *, rrec] = i5_bscs[Nbscs-Ntot+Gla07PixOffset:Nbscs-1, *] Line[0, *, rrec] = FIX(ROUND((elev + yoffset)*yscale))-Gla07PixOffset Line[1, *, rrec] = FIX(ROUND((i_atm_dem + yoffset)*yscale))-Gla07PixOffset LinLeg = [['Surf', '-3'], ['DEM', '3']] ENDFOR END ENDCASE EOF: nrec=irec END ;=================================================================================== ; Lecture d'un profil Niveau 1 correspondant au segment de niveau 2 traité. ;----------------------------------------------------------------------------------- PRO readProfl1, pstate FORWARD_FUNCTION SEARCHPATH par = (*pstate).wstrparam IF (*par).L1type EQ '' THEN RETURN Nshot = (*par).SegLthL1 & Ntot = 280 & Nsamp = 5 Profile = LONARR(Ntot, Nsamp, Nshot) & Line = INTARR(2, Nsamp, Nshot) & infos = 0L res = STRSPLIT((*par).wfilein, (*par).wsep, /EXTRACT, COUNT=cnt) namein = res[cnt - 1] orb = FIX(STRMID(namein, 19, 4)) > 1 & toread = (*par).SegLthL1 orb0 = orb LinLeg = 0 WHILE toread GT 0 DO BEGIN nbread = (*par).SegLthL1 - toread wL1file = (*par).L1type + STRMID(namein, 5, 14) + STRING(orb, FORMAT = '(I4.4)') + STRMID(namein, 23) ; Ouverture éventuelle du fichier ON_IOERROR, EOF path = "@" IF wL1file NE (*par).wL1file THEN BEGIN if (*par).L1lun GT 0 then FREE_LUN, (*par).L1lun path = SEARCHPATH((*par).wliddir, wL1file) openr, L1lun, path, /get_Lun, /SWAP_IF_LITTLE GetGlasFileParms, L1lun, recsize, nbrec, hdlth (*par).wL1file=wL1file & (*par).L1lun = L1lun & (*par).L1recs = recsize & (*par).L1hdl = hdlth ENDIF point_lun, (*par).L1lun, (*par).L1hdl*(*par).L1recs + 4 readu, (*par).L1lun, infos irec0 = (*(*par).datinf)[(*par).iseg] + nbread - infos + (*par).L1hdl IF irec0 LT 0 THEN BREAK; point_lun, (*par).L1lun, (irec0 + (*par).L1hdl)*(*par).L1recs + 4 readu, (*par).L1lun, infos if infos LT (*(*par).datinf)[(*par).iseg] THEN GOTO, EOF (*par).L1datinf = infos readSingleProfl1, pstate, irec0, toread, Profile, Line, LinLeg, nrec, (*par).L1type, (*par).L1lun, (*par).L1recs, 0 toread = toread - nrec EOF: IF path EQ '' THEN BREAK orb = orb + 2 ENDWHILE rNshot = 1 > ((*par).SegLthL1 - toread) IF toread GT 0 THEN BEGIN IF toread EQ (*par).SegLthL1 THEN (*par).L1lun = -1 (*par).wL1file = toread EQ (*par).SegLthL1 ? 'No level 1 data found' : 'Level 1 file ' + wL1file + ' not found' ENDIF Profile = ROTATE(REFORM(Profile[*, *, 0:rNshot-1], Ntot, Nsamp*rNshot), 3) Line = REFORM(Line[*, *, 0:rNshot-1], 2, rNshot*Nsamp) PTR_FREE, (*par).L1Prof, (*par).L1Line (*par).L1Prof = PTR_NEW(Profile) & (*par).L1Line = PTR_NEW(Line) & (*par).L1LinLeg = PTR_NEW(LinLeg) END ;=================================================================================== ; Lectude d'un profil Lidar. On va chercher le segment donné par (*par).iseg ;----------------------------------------------------------------------------------- PRO readProf, pstate par = (*pstate).wstrparam Nshot = LONG((*par).SegLength) & yhight = 22000. & yoffset = 1500. Nlow = 148 & Nmed = 132 yscale = (*par).wysize/yhight pixoffset=1 IF (*par).type EQ 'GLA02' OR (*par).type EQ 'GLA07' THEN BEGIN Ntot = (*par).channel EQ 0 ? 548 : 280 & Nsamp = 5 Profile = LONARR(Ntot, Nsamp, Nshot) & Line = INTARR(2, Nsamp, Nshot) & Layer = 0 irec = (*par).hdlth + (*par).iseg*Nshot readSingleProfl1, pstate, irec, Nshot, Profile, Line, LinLeg, nrec, (*par).type, (*par).lun, (*par).recsize, (*par).channel EQ 0 Profile = ROTATE(REFORM(Profile, Ntot, Nsamp*Nshot), 3) & Line = REFORM(Line, 2, Nshot*Nsamp) (*par).height = Ntot & (*par).width = Nsamp*Nshot ENDIF ELSE BEGIN CASE (*par).type OF 'GLA08' : $ BEGIN Ntot = Nlow + Nmed & Nsamp = 5 badval = 32767 & Nlayer = 5 & LineNsamp = 4 & Nlines = 2 offset = 176 i_bot = INTARR(Nlayer) & i_top = INTARR(Nlayer) Layer = INTARR(Nlayer, 2, Nshot) & Line = INTARR(Nlines, LineNsamp, Nshot) & Profile = 0 i_atm_dem = LONARR(LineNsamp) i_LRpbl_ht = 0 ON_IOERROR, EOF8 FOR irec=0, Nshot-1 DO BEGIN point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + 160 READU, (*par).lun, i_atm_dem point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + 208 READU, (*par).lun, i_LRpbl_ht point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + offset READU, (*par).lun, i_bot, i_top err = WHERE(i_bot - badval EQ 0, cter) & IF cter GT 0 THEN i_bot[err] = -yoffset err = WHERE(i_top - badval EQ 0, cter) & IF cter GT 0 THEN i_top[err] = -yoffset Layer[*, 0, irec] = (i_bot*1.E01 + yoffset)*yscale & Layer[*, 1, irec] = (i_top*1.E01 + yoffset)*yscale Line [0, *, irec] = FIx((i_atm_dem + yoffset)*yscale) Line [1, *, irec] = (i_LRpbl_ht*1.E01 + yoffset)*yscale ENDFOR EOF8: FOR i = 0, Nlayer - 1 DO IF MAX(Layer[i, *, *]) LT 0 THEN BREAK IF i EQ 0 THEN Layer = 0L ELSE Layer = REFORM(Layer[0:(i-1)>0, *, *], i>1, 2, Nshot) Line = REFORM(Line, Nlines, Nshot*LineNsamp) LinLeg = [['DEM', '-3'], ['PBL', '3']] END 'GLA09' : $ BEGIN Ntot = Nlow + Nmed & Nsamp = 5 badval = 32767 & Nlayer = 10 & LineNsamp = 4 & Nlines = 2 LayerNsamp = (*par).Nsamp IF LayerNsamp EQ 1 THEN offset = 200 else offset = 244 i_bot = INTARR(Nlayer, LayerNsamp) & i_top = INTARR(Nlayer, LayerNsamp) Layer = INTARR(Nlayer, 2, LayerNsamp, Nshot) & Line = INTARR(Nlines, LineNsamp, Nshot) & Profile = 0 elev = LONARR(LineNsamp) & i_atm_dem = LONARR(LineNsamp) ON_IOERROR, EOF9 FOR irec=0, Nshot-1 DO BEGIN point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + 168 READU, (*par).lun, elev point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + 184 READU, (*par).lun, i_atm_dem point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + offset READU, (*par).lun, i_bot, i_top err = WHERE(i_bot - badval EQ 0, cter) & IF cter GT 0 THEN i_bot[err] = -yoffset err = WHERE(i_top - badval EQ 0, cter) & IF cter GT 0 THEN i_top[err] = -yoffset Layer[*, 0, *, irec] = (i_bot*1.E01 + yoffset)*yscale + pixoffset Layer[*, 1, *, irec] = (i_top*1.E01 + yoffset)*yscale + pixoffset Line[0, *, irec] = FIX((elev + yoffset)*yscale) & Line[1, *, irec] = FIX((i_atm_dem + yoffset)*yscale) ENDFOR EOF9: FOR i=0, Nlayer - 1 DO IF MAX(Layer[i, *, *, *]) LT 0 THEN BREAK IF i EQ 0 THEN Layer = 0L ELSE Layer = REFORM(Layer[0:(i-1)>0, *, *, *], i>1, 2, Nshot*LayerNsamp) Line = REFORM(Line, Nlines, Nshot*LineNsamp) LinLeg = [['Surf', '-3'], ['DEM', '3']] END 'GLA10' : $ BEGIN IF (*par).aerosol EQ 0 THEN BEGIN yhight = 22000. Nhigh = 0 pNsamp = 4 Nlayer = 10 offsetp = 4640 offsetl = 13700 Profoffset = 250. ; Pour GLA10 : bas du profil à 1km au-dessus du geoide ENDIF ELSE BEGIN yhight = 42500. Nhigh = 268 pNsamp = 1 Nlayer = 9 offsetp = 11312 offsetl = 13868 Profoffset = 500. ; Pour GLA10 : bas du profil à 1km au-dessus du geoide ENDELSE NLines = 1 offseti = 13884 Ntot = Nlow + Nmed + Nhigh & Nsamp = 5 yscale = Ntot/yhight Gla10PixOffset = ROUND(yscale*Profoffset) pbadval = 2147483647L & badval = 32767 LayerNsamp = pNsamp & LineNsamp = 1 Profile = LONARR(Ntot, pNsamp, Nshot) & Layer = INTARR(Nlayer, 2, LayerNsamp, Nshot) & Line = INTARR(NLines, LineNsamp, Nshot) i_ext_prof = LONARR(Ntot, pNsamp) i_bot = INTARR(Nlayer, LayerNsamp) & i_top = INTARR(Nlayer, LayerNsamp) i_pbl = INTARR(LineNsamp) ON_IOERROR, EOF10 FOR irec=0, Nshot-1 DO BEGIN point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + offsetp READU, (*par).lun, i_ext_prof bad = WHERE(i_ext_prof EQ pbadval, ct) & IF ct GT 0 THEN i_ext_prof[bad] = 0 Profile[0:Ntot-Gla10PixOffset-1, *, irec] = i_ext_prof[Gla10PixOffset:Ntot-1, *] point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + offsetl READU, (*par).lun, i_bot, i_top err = WHERE(i_bot - badval EQ 0, cter) & IF cter GT 0 THEN i_bot[err] = -yoffset err = WHERE(i_top - badval EQ 0, cter) & IF cter GT 0 THEN i_top[err] = -yoffset Layer[*, 0, *, irec] = (i_bot*1.E01 + yoffset)*yscale Layer[*, 1, *, irec] = (i_top*1.E01 + yoffset)*yscale point_lun, (*par).lun, (*par).recsize*((*par).hdlth + (*par).iseg*Nshot + irec) + offseti READU, (*par).lun, i_pbl Line[*, *, irec] = (i_pbl*1.E01 + yoffset)*yscale - Gla10PixOffset ENDFOR EOF10: newsamp=Nsamp*4 np = LONARR(Ntot, newsamp, Nshot) FOR i = 0, newsamp - 1 DO np[*, i, *] = Profile[*,pNsamp*i/newsamp,*] Profile = np Profile = ROTATE(REFORM(Profile, Ntot, newsamp*Nshot), 3) MaxLay = INTARR(Nlayer)*0 FOR i=0, Nlayer - 2 DO IF MAX(Layer[i, *, *, *]) GT 0 THEN MaxLay[i] = 1 LayerOK = WHERE(MaxLay GT 0, ctok) IF ctok EQ 0 THEN Layer = 0L ELSE Layer = REFORM(Layer[LayerOK, *, *, *], ctok, 2, Nshot*LayerNsamp) Line = REFORM(Line, Nlines, Nshot*LineNsamp) LinLeg = ['PBL'] END 'GLA11' : $ BEGIN ON_IOERROR, EOF11 EOF11: END ENDCASE (*par).height = Ntot & (*par).width = Nsamp*(*par).SegLthL1 ENDELSE PTR_FREE, (*par).Profile, (*par).Layer, (*par).Line (*par).Profile = PTR_NEW(Profile) & (*par).Layer = PTR_NEW(Layer) (*par).Line = PTR_NEW(Line) & (*par).LinLeg = PTR_NEW(LinLeg) END ;**************************************************************************************************** ;************************************** INTERFACE GRAPHIQUE ************************************** ;**************************************************************************************************** FUNCTION Create_GDI FORWARD_FUNCTION ResetPlot !QUIET = 1 sizeunit = 40 & gap = 5 ; Creates TopLevelBase ; =========================== TopLevelBase = WIDGET_BASE(TITLE='LIDAR Data Vizualisation', MBAR=MenuBar, /ROW) ; Creates menus ; ================== FileMenu = WIDGET_BUTTON(MenuBar , VALUE='File', /MENU, EVENT_PRO='FileMenu_Event') HelpMenu = WIDGET_BUTTON(MenuBar , VALUE='Help', /MENU, EVENT_PRO='HelpMenu_Event') OpenOption = WIDGET_BUTTON(FileMenu, VALUE='Open', UVALUE='file') ExitOption = WIDGET_BUTTON(FileMenu, VALUE='Exit', UVALUE='exit') HelpOption = WIDGET_BUTTON(HelpMenu, VALUE='Help', UVALUE='help') ; Creates sub-base ParamBase: Column with all buttons ; ====================================== ParamBase = WIDGET_BASE(TopLevelBase, /COLUMN, EVENT_PRO='Param_Event') ; Creates sub-base FileBase Selection du fichier ; ====================================== FileBase = WIDGET_BASE(ParamBase, /COLUMN, FRAME=2) TypeLabel = WIDGET_LABEL(FileBase, VALUE='Selected Product', UVALUE='label', /DYNAMIC_RESIZE) FileButt = WIDGET_BUTTON(FileBase, XSIZE=sizeunit*6, VALUE='FILE', UVALUE='file') FileLabel = WIDGET_LABEL(FileBase, VALUE='Selected File', UVALUE='label', /DYNAMIC_RESIZE) DateLabel = WIDGET_LABEL(FileBase, VALUE='Date', UVALUE='date', /DYNAMIC_RESIZE) ; Creates subbase Directorybase: All buttons to select a directory: ; ======================================== ButtBase1 = WIDGET_BASE(ParamBase, COLUMN=2, FRAME=2, UNAME='ButtBase1', SENSITIVE=0) PrevButt = WIDGET_BUTTON(ButtBase1, XSIZE=sizeunit*3, VALUE='PREV', UVALUE='prev') NextButt = WIDGET_BUTTON(ButtBase1, XSIZE=sizeunit*3, VALUE='NEXT', UVALUE='next') ButtBase2 = WIDGET_BASE(ParamBase, COLUMN=2, FRAME=2, UNAME='ButtBase2', SENSITIVE=0) PNGSave = WIDGET_BUTTON(ButtBase2, XSIZE=sizeunit*3, VALUE='PNG' , UVALUE='png' ) epsSave = WIDGET_BUTTON(ButtBase2, XSIZE=sizeunit*3, VALUE='EPS' , UVALUE='eps' ) ; Creates the visualisation window (map of orbits) ; ==================================================== draw_xsize=540 & draw_ysize=270 Draw = WIDGET_DRAW(TopLevelBase, XSIZE=draw_xsize, YSIZE=draw_ysize, $ EVENT_PRO='Fenetre_Event', /BUTTON_EVENTS, UVALUE='draw') ; Realize TopLevelBase WIDGET_CONTROL, TopLevelBase, /REALIZE WIDGET_CONTROL, Draw, GET_VALUE=draw_id CASE !VERSION.OS_FAMILY OF 'MacOS' : graph='MAC' 'unix' : graph='X' 'Windows' : graph='WIN' ENDCASE sep=PATH_SEP() ; ; Initialisation graphique ; NC = ResetPlot(graph) win_xsize = 810 & win_ysize = 583 & SegRec2 = 2 & SegLthL2 = 50 ptrtab = PTRARR(4) struct_widgets={ $ wFileBase:FileBase, wTypeLabel:TypeLabel, wFileButt:FileButt, wFileLabel:FileLabel, wDateLabel:DateLabel, $ wParamBase:ParamBase, wVisuBase:-1l, wMinF:-1l, wMaxF:-1l, wLines:-1l, wsuper:-1l, wunic:-1l, $ wPrev:PrevButt, wNext:NextButt, wpng:PNGsave, weps:epsSave, wRefFile:-1L} struct_parametres={$ wpathHelp:'Help_VisuLidar.txt', wliddir:'', wpathin:'', wfilein:'', wpathout:'', wL1file:'', wfileout:'', $ l1type:'', ftype:'', wsep:sep, graph:graph, xsize:draw_xsize, ysize:draw_ysize, wxsize:0, wysize:0, wid:Draw_id, $ numwin:-1, l1win0:-1, l1win:-1, disp:0l, NC:NC, L1datinf:0L, minv:0., maxv:0., lun:-1l, L1lun:-1l, $ latc:PTR_NEW(), lonc:PTR_NEW(), lats:PTR_NEW(), lons:PTR_NEW(), rlat:PTR_NEW(), rlon:PTR_NEW(), SolElev:PTR_NEW(), $ datinf:PTR_NEW(), Profile:PTR_NEW(), Layer:PTR_NEW(), Line:PTR_NEW(), L1Prof:PTR_NEW(), L1Line:PTR_NEW(), $ SegLthL2:0, SegLthL1:0, SegLength:0, Nseg:0, Nsamp:0, channel:0, typdisp:0, aerosol:0, iseg:0, altitudes:'', $ recsize:0l, hdlth:0, nbrec:0, L1recs:0l, L1hdl:0l, ilat:0, ilon:0, SegRec2:0, SegRec:0, prof_idx:PTR_NEW(), $ maplim:[-89., -180., 89., 180.], width:0, height:0, type:'', torig:0.D, LidName:'', LidNiv:0, $ LinLeg:PTR_NEW(), L1LinLeg:PTR_NEW() } state = { $ wstrWidg :PTR_NEW(struct_widgets , /NO_COPY), $ ; Pointer on a structure that defines widgets wstrParam:PTR_NEW(struct_parametres, /NO_COPY) $ ; Pointer on a structure that defines parameters } pstate=PTR_NEW(state, /NO_COPY) WIDGET_CONTROL, TopLevelBase, SET_UVALUE=pstate XMANAGER, 'LIDAR Analysis', TopLevelBase, /NO_BLOCK, /JUST_REG RETURN, pstate END ; ***** END OF Create_GDI ***** ;**************************************************************************************************** ;************************************** GESTION DES EVENEMENTS ************************************** ;**************************************************************************************************** ;========================================================================================= ; ***** CommonParamEvent : PRO ***** ; Gestion des évènements communs à tous les lidars ; * event (I) évenement déclencheur ;----------------------------------------------------------------------------------------- PRO CommonParamEvent, event FORWARD_FUNCTION ResetPlot, fileselect, DateTimeString WIDGET_CONTROL, event.top, GET_UVALUE=pstate WIDGET_CONTROL, event.id , GET_UVALUE=uval WIDGET_CONTROL, event.id , GET_VALUE=val par = (*pstate).wstrParam & wid = *(*pstate).wstrWidg DEVICE, WINDOW_STATE=winstate CASE uval OF 'file' : $ IF fileselect(pstate) THEN BEGIN WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='ButtBase1'), /SENSITIVE WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='ButtBase2'), /SENSITIVE Display, pstate ENDIF 'prev' : $ BEGIN (*par).iseg = ((*par).iseg - 1) > 0 & Display, pstate END 'next' : $ BEGIN (*par).iseg = ((*par).iseg + 1) < ((*par).Nseg - 1) & Display, pstate END 'chann': $ BEGIN (*par).channel = val & Display, pstate, MAPFIRST=-1 END 'style': $ BEGIN DisplayProduct, pstate & DisplayL1RelatedProfile, pstate END 'lines': $ BEGIN DisplayProduct, pstate & DisplayL1RelatedProfile, pstate END 'min' : IF (*par).L1type EQ '' THEN DisplayProduct, pstate ELSE DisplayL1RelatedProfile, pstate 'max' : IF (*par).L1type EQ '' THEN DisplayProduct, pstate ELSE DisplayL1RelatedProfile, pstate 'png' : $ BEGIN ; Creation d'une image PNG IF (*par).LidNiv EQ 1 THEN BEGIN WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='Channel'), GET_VALUE=chnl fo = (*par).wfileout+'_'+STRING((*par).iseg + 1, FORMAT='(I2.2)')+'_'+(chnl EQ 0 ? '0532' : '1064')+'.PNG' ENDIF ELSE fo = (*par).wfileout+'_'+STRING((*par).iseg + 1, FORMAT='(I2.2)')+'.PNG' fileout = DIALOG_PICKFILE(/write, title='Select Filename for PNG image file', file=fo, GET_PATH=pathout, path=(*par).wpathout) IF STRLEN(fileout) EQ 0 THEN RETURN (*par).wpathout = pathout DEVICE, DECOMPOSE=1 WSET, (*par).numwin WRITE_PNG, fileout, TVRD(/TRUE) DEVICE, DECOMPOSE=0 END 'eps' : $ BEGIN ; Creation d'un fichier poscript encapsulé IF (*par).LidNiv EQ 1 THEN BEGIN WIDGET_CONTROL, WIDGET_INFO(event.top, FIND_BY_UNAME='Channel'), GET_VALUE=chnl fo = (*par).wfileout+'_'+STRING((*par).iseg + 1, FORMAT='(I2.2)')+'_'+(chnl EQ 0 ? '0532' : '1064')+'.EPS' ENDIF ELSE fo = (*par).wfileout+'_'+STRING((*par).iseg + 1, FORMAT='(I2.2)')+'.EPS' fileout = DIALOG_PICKFILE(/write, title='Select Filename for EPS image file', file=fo, GET_PATH=pathout, path=(*par).wpathout) IF STRLEN(fileout) EQ 0 THEN RETURN (*par).wpathout = pathout SET_PLOT, 'PS' DEVICE, FILENAME=fileout, /COLOR, BITS=8, /PORTRAIT, xsize=17., ysize=17., xoff=2., yoff=2., /ENCAPSULATED DisplayProduct, pstate, /PS DEVICE, /CLOSE nc = ResetPlot((*par).graph) END ELSE: ENDCASE END ;========================================================================================= ; ***** HelpMenu_Event : PRO ***** ; Read and display Help File ; * event (I) évenement déclencheur ;----------------------------------------------------------------------------------------- PRO HelpMenu_Event, event WIDGET_CONTROL, event.top, GET_UVALUE=pstate ; Tests whether file is easily located. If not, asks for the location, and keeps information IF NOT FILE_TEST((*(*pstate).wstrParam).wpathHelp) THEN $ (*(*pstate).wstrParam).wpathHelp = DIALOG_PICKFILE(FILTER='Help_visu_lidar.txt', /MUST_EXIST, TITLE='Locate Help File') XDISPLAYFILE, (*(*pstate).wstrParam).wpathHelp, TITLE='Short Help', DONE_BUTTON='Quit', WIDTH=100 END ; ***** END OF HelpMenu_Event ***** ;========================================================================================= ; ***** FileMenu_Event : PRO ***** ; File menu only contains "exit" in this program ; * event (I) évenement déclencheur ;----------------------------------------------------------------------------------------- PRO FileMenu_Event, event WIDGET_CONTROL, event.id, GET_UVALUE=uval DEVICE, WINDOW_STATE=winstate CASE uval OF 'exit': $ BEGIN FOR id=0, N_ELEMENTS(winstate) - 1 DO IF winstate[id] THEN WDELETE, id ; Erase all opened windows WIDGET_CONTROL, event.top, /DESTROY ; Erase main window HEAP_GC ; Clean Heap Memory END ELSE: CommonParamEvent, event ENDCASE END ; ***** END OF FileMenu_Event ***** ;=================================================================================================== ; ***** Fenetre_Event : PRO ***** ; A click on the window. Look for the closest targets and displays it. ; * event (I) évenement déclencheur ;---------------------------------------------------------------------------------------------------- PRO Fenetre_Event, event IF event.release THEN RETURN WIDGET_CONTROL, event.top, GET_UVALUE=pstate par = (*pstate).wstrParam lon0 = ((*par).maplim[1] + event.x*((*par).maplim[3] - (*par).maplim[1])/(*par).xsize + 360D) MOD 360.D lat0 = (*par).maplim[0] + event.y*((*par).maplim[2] - (*par).maplim[0])/(*par).ysize tempo = MIN((*(*par).latc - lat0)^2 + (*(*par).lonc - lon0)^2, imin) (*par).iseg = imin Display, pstate END ; ***** END OF Fenetre_Event ***** ;**************************************************************************************************** ;************************************** GESTION DES AFFICHAGES ************************************** ;**************************************************************************************************** ;=========================================================================================== ; ***** ResetPlot : FUNCTION ***** ; Réinitialisation de l'environnement graphique. Utile après une génération PS, par exemple. ; * graph (I) device à initialiser ;------------------------------------------------------------------------------------------- FUNCTION ResetPlot, graph SET_PLOT, graph DEVICE, DECOMPOSE=0 LOADCT, 13 TVLCT, r, g, b, /GET nc = N_ELEMENTS(r) r[0] = 0 & g[0] = 0 & b[0] = 0 & r[nc-1] = 255 & g[nc-1] = 255 & b[nc-1] = 255 TVLCT, r, g, b !P.COLOR = 0 & !P.BACKGROUND = nc-1 return, nc END ; ***** END OF ResetPlot ***** ;========================================================================================= ; ***** Display : PRO ***** ; Visualisation des données et paramètres d'un segment du produit principal ; * pstate (I) pointeur sur la zone contenant les données à visualiser ; * /MAPFIRST (I) affichage de la carte d'abord (si MAPFIRST=-1, pas d'affichage de carte) ; * /RESET (I) si simple réinitialisation des différentes fenêtres d'affichage ; * /NOLEVEL1 (I) si pas d'affichage du segment de niveau 1 correspondant ; * /LEVEL1 (I) si affichage uniquement du segment de niveau 1 correspondant ;----------------------------------------------------------------------------------------- PRO Display, pstate, MAPFIRST=mapfirst, RESET=reset, NOLEVEL1=nolevel1, LEVEL1=level1 IF KEYWORD_SET(RESET) THEN BEGIN DisplayMap, pstate, /RESET & DisplayProduct, pstate, /RESET RETURN ENDIF IF KEYWORD_SET(MAPFIRST) THEN IF mapfirst NE -1 THEN DisplayMap, pstate IF NOT KEYWORD_SET(LEVEL1) THEN BEGIN readProf, pstate DisplayProduct, pstate ENDIF IF NOT KEYWORD_SET(NOLEVEL1) THEN BEGIN readProfl1, pstate DisplayL1RelatedProfile, pstate ENDIF IF NOT KEYWORD_SET(MAPFIRST) THEN DisplayMap, pstate END ; ***** END OF Display ***** ;======================================================================================== ; ***** DisplayMap : PRO ***** ; Visualisation d'une carte gobale. On met les orbites dessus. Le code couleur est ; fonction de l'angle solaire. Les segments en noir indiquent une mesure de nuit. ; * pstate (I) pointeur sur la zone contenant les paramètres à visualiser sur la carte ; * /RESET (I) si simple réaffichage de la carte sans les orbites et segments ;---------------------------------------------------------------------------------------- PRO DisplayMap, pstate, RESET=reset FORWARD_FUNCTION DateTimeRange par = *(*pstate).wstrparam WSET, par.wid MAP_SET, 0., 0., LIMIT=par.maplim, /CONTINENTS, /NOBORDER, POS=[0., 0., 1., 1.] IF KEYWORD_SET(RESET) THEN RETURN color = BYTE(SIN(*par.SolElev > 0.)*(par.NC - 3)) + 1 FOR i = 0, par.Nseg - 1 DO BEGIN FOR j = i*par.SegRec, (i+1)*par.SegRec - 1 DO BEGIN OPLOT, (*par.rlon)[j:j+1], (*par.rlat)[j:j+1], COL=color[i], THICK=2 ENDFOR ENDFOR OPLOT, (*par.lonc), (*par.latc), PSYM=5 iseg = par.iseg FOR j = iseg*par.SegRec, (iseg+1)*par.SegRec - 1 DO BEGIN OPLOT, (*par.rlon)[j:j+1], (*par.rlat)[j:j+1], COL=par.NC-3, THICK=3 ENDFOR WIDGET_CONTROL, (*(*pstate).wstrWidg).wDateLabel, $ SET_VALUE=DateTimeRange(par.torig, (*(par).datinf)[iseg], par.SegLthL1) END ; ***** END OF DisplayMap ***** ;========================================================================================= ; ***** DisplayProduct : PRO ***** ; Visualisation des données d'un segment du produit principal dans une fenêtre ; * pstate (I) pointeur sur la zone contenant le profil à visualiser ; * /PS (I) si génération Postscript ; * /RESET (I) si simple effacement de la fenêtre ;----------------------------------------------------------------------------------------- PRO DisplayProduct, pstate, PS=ps, RESET=reset wid = *(*pstate).wStrWidg & par = *(*pstate).wStrParam IF SIZE(*(par).Profile, /N_DIMENSIONS) GT 0 THEN BEGIN WIDGET_CONTROL, wid.wMinF, GET_VALUE=vmin & WIDGET_CONTROL, wid.wMaxF, GET_VALUE=vmax END DisplayWindow, pstate, vmin, vmax, ps, reset END ; ***** END OF DisplayProduct ***** ;========================================================================================= ; ***** DisplayL1RelatedProfile : PRO ***** ; Visualisation du profil du segment de niveau 1 correspondant dans une fenêtre ; * pstate (I) pointeur sur la zone contenant le profil Lidar à visualiser ; * /RESET (I) si simple effacement de la fenêtre ;----------------------------------------------------------------------------------------- PRO DisplayL1RelatedProfile, pstate, RESET=reset FORWARD_FUNCTION GetPathFile wid = *(*pstate).wStrWidg & par = *(*pstate).wStrParam IF par.L1type EQ '' THEN RETURN WIDGET_CONTROL, wid.wMinF, GET_VALUE=vmin & WIDGET_CONTROL, wid.wMaxF, GET_VALUE=vmax DisplayWindow, pstate, vmin, vmax, ps, reset, /LEVEL1 ; DisplayWindow, pstate, par.minv, par.maxv, ps, reset, /LEVEL1 WIDGET_CONTROL, wid.wRefFile, SET_VALUE=GetPathFile(par.wL1file,SUFFIX=".hdf") END ; ***** END OF DisplayL1RelatedProfile ***** ;========================================================================================= ; ***** FillLayer : PRO ***** ; Affichage d'une couche ; * col (I) couleur d'affichage de la couche ; * xsc (I) échelle en x ; * bot (I) vecteur des altitudes du bas de la couche ; * top (I) vecteur des altitudes du bas de la couche ;----------------------------------------------------------------------------------------- PRO FillLayer, col, xsc, style, bot, top, SUPER=super IF style EQ 0 THEN FillLayer0, col, xsc, bot, top, super $ ELSE FillLayer1, col, xsc, bot, top, super END ; ***** END OF FillLayer ***** PRO FillLayer0, col, xsc, bot, top, super FOR k = 0, N_ELEMENTS(bot) - 1 DO BEGIN IF top[k] LT 0 THEN CONTINUE x = INTARR(4) & y = INTARR(4) x[0:1] = xsc*(indgen(2) + k) & y[0:1] = bot[k] x[2:3] = REVERSE(x[0:1]) & y[2:3] = top[k] IF NOT KEYWORD_SET(SUPER) THEN POLYFILL, x, y, COLOR=col, /DEVICE $ ELSE BEGIN POLYFILL, x, y, COLOR=col, /LINE_FILL, SPACING=0.15, /DEVICE, ORIENTATION=+45 POLYFILL, x, y, COLOR=col, /LINE_FILL, SPACING=0.15, /DEVICE, ORIENTATION=-45 ENDELSE ENDFOR IF KEYWORD_SET(SUPER) THEN BEGIN deb = 1 FOR k = 0, N_ELEMENTS(bot) - 2 DO BEGIN IF top[k] LT 0 THEN CONTINUE IF deb THEN PLOTS, [xsc*k, xsc*k], [bot[k], top[k]], COLOR=col, /DEVICE PLOTS, [xsc*k, xsc*(k+1)], [bot[k], bot[k]], COLOR=col, /DEVICE PLOTS, [xsc*k, xsc*(k+1)], [top[k], top[k]], COLOR=col, /DEVICE IF top[k] LT bot[k+1] OR bot[k] GT top[k+1] THEN BEGIN deb = 1 PLOTS, [xsc*(k+1), xsc*(k+1)], [bot[k], top[k]], COLOR=col, /DEVICE ENDIF ELSE BEGIN deb = 0 PLOTS, [xsc*(k+1), xsc*(k+1)], [bot[k], bot[k+1]], COLOR=col, /DEVICE PLOTS, [xsc*(k+1), xsc*(k+1)], [top[k], top[k+1]], COLOR=col, /DEVICE ENDELSE ENDFOR IF top[k] GE 0 THEN BEGIN PLOTS, [xsc*k, xsc*(k+1)], [bot[k], bot[k]], COLOR=col, /DEVICE PLOTS, [xsc*k, xsc*(k+1)], [top[k], top[k]], COLOR=col, /DEVICE PLOTS, [xsc*(k+1), xsc*(k+1)], [bot[k], top[k]], COLOR=col, /DEVICE ENDIF ENDIF END ; ***** END OF FillLayer0 ***** PRO FillLayer1, col, xsc, bot, top, super FOR k = 0, N_ELEMENTS(bot) - 2 DO BEGIN IF bot[k] LT 0 THEN CONTINUE FOR j = k + 1, N_ELEMENTS(bot) - 1 DO IF bot[j] LT 0 THEN BREAK IF bot[k] LE 0 THEN CONTINUE FOR j = k + 1, N_ELEMENTS(bot) - 1 DO IF bot[j] LE 0 THEN BREAK IF bot[k] LT 0 THEN CONTINUE FOR j = k + 1, N_ELEMENTS(bot) - 1 DO IF bot[j] GT top[j-1] OR top[j] LT bot[j-1] THEN BREAK x = INTARR(2*(j - k)) & y = INTARR(2*(j - k)) x[0:j-k-1] = xsc*(0.5 + indgen(j - k) + k) & y[0:j-k-1] = bot[k:j-1] x[j-k:2*(j-k)-1] = REVERSE(x[0:j-k-1]) & y[j-k:2*(j-k)-1] = REVERSE(top[k:j-1]) IF NOT KEYWORD_SET(SUPER) THEN POLYFILL, x, y, COLOR=col, /DEVICE $ ELSE BEGIN POLYFILL, x, y, COLOR=col, /LINE_FILL, SPACING=0.15, /DEVICE, ORIENTATION=+45 POLYFILL, x, y, COLOR=col, /LINE_FILL, SPACING=0.15, /DEVICE, ORIENTATION=-45 FOR i = 0, N_ELEMENTS(x) - 2 DO PLOTS, [x[i], x[i+1]], [y[i], y[i+1]], COLOR=col, /DEVICE PLOTS, [x[i], x[0]], [y[i], y[0]], COLOR=col, /DEVICE ENDELSE k = j ENDFOR END ; ***** END OF FillLayer1 ***** ;========================================================================================= ; ***** DrawLine : PRO ***** ; Affichage d'une ligne ; * col (I) couleur d'affichage de la ligne ; * xsc (I) échelle en x ; * line (I) vecteur des altitudes ;----------------------------------------------------------------------------------------- PRO DrawLine, col, xsc, line, LEG=leg, NUM=num FOR j = 0, N_ELEMENTS(line) - 2 DO PLOTS, [xsc*j, xsc*(j+1)-1], line[j]*[1, 1], COLOR=col, /DEVICE IF KEYWORD_SET(LEG) THEN BEGIN jj = N_ELEMENTS(line) & xx = xsc*j IF KEYWORD_SET(NUM) THEN l1 = FIX(leg[0]) ELSE l1 = 0 py = l1 GE 0 ? 3 : -11 l1 = ABS(l1) & jjj = l1 EQ 1 ? 0 : l1 EQ 3 ? jj - 1 : jj/2 & px = jjj*xsc + 10*((2-l1)/2) algn = l1 EQ 1 ? 0. : l1 EQ 3 ? 1. : 0.5 XYOUTS, px, line[jjj] + py, leg[0], COLOR=col, /DEVICE, ALIGNMENT=algn ENDIF END ; ***** END OF DrawLine ***** ;**************************************************************************************************** ;****************************************** UTILITAIRES ********************************************* ;**************************************************************************************************** ;============================================================================================ ; ***** ElevSol : FUNCTION ***** ; Calcul de l'élévation solaire en fonction du jour julien, de l'heure TU et de la position. ; Tous les angles, y compris latitutde et longitude, sont en radians. ; p (I) structure contenant les paramètres spatio-temporels (date, heure, lat, lon) ;-------------------------------------------------------------------------------------------- FUNCTION ElevSol, p FORWARD_FUNCTION ElevationSolaire RETURN, ElevationSolaire(p.torig + *(p.datinf)/86400.d, *(p.lats), *(p.lons)) END ; ***** END OF ElevSol ***** ;=========================================================================================== ; ***** DateTimeRange : FUNCTION ***** ; Formate les dates et heures de début et fin d'un intervalle de temps ; orig (I) origine des temps, sous la forme d'une date julienne ; t (I) temps, en s., depuis l'origne des temps ; delay (I) durée, en s., de l'intervalle de temps ;------------------------------------------------------------------------------------------- FUNCTION DateTimeRange, orig, t, delay FORWARD_FUNCTION DateTimeString RETURN, 'du ' + DateTimeString(orig, t) + ' au ' + DateTimeString(orig, t + delay) END ; ***** END OF DateTimeRange ***** ;=========================================================================================== ; ***** DateTimeString : FUNCTION ***** ; Formatage date/heure ; orig (I) origine des temps, sous la forme d'une date julienne ; t (I) temps, en s., depuis l'origne des temps ;------------------------------------------------------------------------------------------- FUNCTION DateTimeString, orig, t jul = orig + t/86400.d CALDAT, jul, month, day, year, hour, minu, sec RETURN, STRING(year, month, day, hour, minu, sec, FORMAT='(I4.4, 2("/", I2.2), ",", I2.2, 2(":", I2.2))') END ; ***** END OF DateTimeRange ***** ;**************************************************************************************************** ;************************************** GESTION DES AFFICHAGES ************************************** ;**************************************************************************************************** ;========================================================================================= ; ***** CreateL1Window : PRO ***** ; Création et affichage de la fenêtre du segment de niveau 1 associé ; * pstate (I) pointeur sur la zone contenant les données à visualiser ; * type (I) GLA02 si type = 0, GLA07 sinon ;----------------------------------------------------------------------------------------- PRO CreateL1Window, pstate, type wid = *(*pstate).wStrWidg & par = (*pstate).wstrparam IF type EQ 0 THEN (*par).L1type = 'GLA02' ELSE (*par).L1type = 'GLA07' SetMinMax, *par, (*par).L1type WIDGET_CONTROL, wid.wMinF, SET_VALUE=(*par).minv & WIDGET_CONTROL, wid.wMaxF, SET_VALUE=(*par).maxv resol = GET_SCREEN_SIZE() WINDOW, (*par).l1win, XSIZE=(*par).wxsize, YSIZE=(*par).wysize, YPOS=(resol[1]-2*(28+(*par).wysize) > 0), TITLE='Related '+(*par).L1type+' segment' Display, pstate, /LEVEL1 END ; ***** END OF CreateL1Window ***** ;========================================================================================= ; ***** DisplayWindow : PRO ***** ; Visualisation des données d'un segment dans une fenêtre (profils, couches et lignes) ; * pstate (I) pointeur sur la zone contenant les données à visualiser ; * vmin, vmax (I) valeurs min et max pour les profils Lidar ; * ps (I) défini si génération Postscript ; * reset (I) défini si simple effacement de la fenêtre ; * /LEVEL1 (I) s'il s'agit du profil d'un segment de niveau 1 associé ;----------------------------------------------------------------------------------------- PRO DisplayWindow, pstate, vmin, vmax, ps, reset, LEVEL1=level1 wid = *(*pstate).wStrWidg & par = *(*pstate).wStrParam Layer = *(par).Layer IF KEYWORD_SET(LEVEL1) THEN BEGIN Profile = *(par).L1Prof & Line = *(par).L1Line & LinLeg = *(par).L1LinLeg IF NOT KEYWORD_SET(PS) THEN BEGIN DEVICE, WINDOW_STATE=winstate IF NOT winstate[par.l1win] THEN CreateL1Window, pstate, par.L1type ENDIF numwin = par.l1win ENDIF ELSE BEGIN Profile = *(par).Profile & Line = *(par).Line & LinLeg = *(par).LinLeg IF NOT KEYWORD_SET(PS) THEN setWindow, pstate, par.width, par.height, par.type par = *(*pstate).wStrParam ; Car possible modif de wysize numwin = par.numwin ENDELSE IF NOT KEYWORD_SET(PS) THEN BEGIN WSET, numwin & WSHOW, numwin ENDIF colbg = 0 & colfg = !P.BACKGROUND xsz = par.wxsize & ysz = par.wysize ymid = 12000.d & yhg1 = 22000.d & yhg2 = 42500.d ymax = par.height GT 300 ? yhg2 : yhg1 yscale = ysz/ymax & yoffset = 1500 s_alts = STRSPLIT(par.altitudes, '@', /EXTRACT) & alts = (1000*DOUBLE(s_alts) + yoffset)*yscale IF KEYWORD_SET(RESET) OR KEYWORD_SET(LEVEL1) THEN TV, MAKE_ARRAY(xsz, ysz, /BYTE, VALUE=colbg) IF KEYWORD_SET(RESET) OR (KEYWORD_SET(LEVEL1) AND par.L1lun LT 0) THEN RETURN IF SIZE(Profile, /N_DIMENSIONS) GT 0 AND (par.typdisp LE 1 OR KEYWORD_SET(LEVEL1)) THEN BEGIN scale = (par.NC - 3) /(vmax - vmin) visu = BYTE((((ALOG(Profile > 1) - vmin) > 0)*scale) < (par.NC - 3)) + 1 TV, visu super = 1 ENDIF ELSE BEGIN TV, MAKE_ARRAY(xsz, ysz, /BYTE, VALUE=colbg) super = 0 ENDELSE IF SIZE(Layer, /N_DIMENSIONS) GT 0 AND NOT KEYWORD_SET(LEVEL1) AND par.typdisp NE 1 THEN BEGIN style = WIDGET_INFO(WIDGET_INFO(wid.wVisuBase, FIND_BY_UNAME='Style'), /DROPLIST_SELECT) Nsamp = (SIZE(Layer, /DIMENSIONS))[2] NLayer = (SIZE(Layer, /DIMENSIONS))[0] xscale = xsz/Nsamp FOR i = 0, NLayer - 1 DO FillLayer, ((10 - i)*(par.NC - 3))/10, xscale, style, Layer[i, 0, *], Layer[i, 1, *], SUPER=super ENDIF IF KEYWORD_SET(LEVEL1) THEN BEGIN WIDGET_CONTROL, wid.wsuper, GET_VALUE=choix IF choix THEN BEGIN IF SIZE(Layer, /N_DIMENSIONS) GT 0 THEN BEGIN style = WIDGET_INFO(WIDGET_INFO(wid.wVisuBase, FIND_BY_UNAME='Style'), /DROPLIST_SELECT) Nsamp = (SIZE(Layer, /DIMENSIONS))[2] NLayer = (SIZE(Layer, /DIMENSIONS))[0] xscale = xsz/Nsamp FOR i = 0, NLayer - 1 DO FillLayer, ((10 - i)*(par.NC-3))/10, xscale, style, Layer[i, 0, *], Layer[i, 1, *], /SUPER ENDIF ENDIF ENDIF WIDGET_CONTROL, wid.wLines, GET_VALUE=choix IF choix[1] THEN FOR i = 0, N_ELEMENTS(alts) - 1 DO DrawLine, colfg, xsz, alts[i]*[1, 1], LEG=[s_alts[i], '1'] IF choix[0] AND SIZE(Line, /N_DIMENSIONS) GT 0 THEN BEGIN NLines = (SIZE(Line, /DIMENSIONS))[0] & Nsamp = (SIZE(Line, /DIMENSIONS))[1] WIDGET_CONTROL, wid.wLines, GET_VALUE=choix FOR i = 0, NLines - 1 DO DrawLine, ((10 - i)*(par.NC - 3))/10, FLOAT(xsz)/Nsamp, line[i, *], LEG=[LinLeg[*,i]] ENDIF END ; ***** END OF DisplayWindow ***** ;**************************************************************************************************** ;********************************** SELECTION DES FICHIERS GLAS ************************************* ;**************************************************************************************************** ;========================================================================================= ; ***** FileSelect : FUNCTION ***** ; Choix d'un fichier GLAS; lecture des informations générales. On découpe en segments. ; On lit lat, lon et heure et on en déduit l'élévation solaire. ; * pstate (I) pointeur sur la zone contenant la structure où renseigner ces informations ;----------------------------------------------------------------------------------------- FUNCTION FileSelect, pstate FORWARD_FUNCTION GetGlasFile, ElevSol IF NOT GetGlasFile(pstate) THEN RETURN, 0; Attention : Pb.: nbrec non renseigné pour GLA07 par = (*pstate).wstrParam & wid = (*pstate).wstrWidg & lun = (*par).lun idat = 1 & ilat = (*par).ilat & ilon = (*par).ilon & ninf = ilon + 1 infos = LONARR(ninf) & Sinfos = LONARR(ninf, 10000) iseg = 0 & SegLength = LONG((*par).SegLength) ON_IOERROR, EOF1 WHILE NOT EOF(lun) DO BEGIN POINT_LUN, lun, (iseg*SegLength/2 + (*par).hdlth)*(*par).recsize READU, lun, infos Sinfos[*, iseg] = infos iseg = iseg + 1 ENDWHILE EOF1: Nseg = iseg/2 Sinfos = Sinfos[*, 2*INDGEN(Nseg)] & nbrec = Nseg*(*par).SegRec & lth = SegLength/(*par).SegRec llinfos = LONARR(ninf) & latloninf = LONARR(ninf - ilat, nbrec + 1) ON_IOERROR, EOF2 FOR irec = 0, nbrec - 1 DO BEGIN POINT_LUN, lun, (irec*lth + (*par).hdlth)*(*par).recsize READU, lun, llinfos latloninf[*, irec] = llinfos[ilat:ilon] ENDFOR EOF2: latloninf[*, nbrec] = llinfos[ilat:ilon] rlat = REFORM(latloninf[0, *]*1.E-6) & rlon = REFORM(latloninf[ilon - ilat, *]*1.E-6) ; Lats et lons donnent les position de l'extremite des segments, latc et lonc celle du milieu du segment lats = REFORM(Sinfos[ilat, *]*1.E-6) & lons = REFORM(Sinfos[ilon, *]*1.E-6) latc = rlat[INDGEN(Nseg)*(*par).SegRec + (*par).SegRec2] & lonc = rlon[INDGEN(Nseg)*(*par).SegRec + (*par).SegRec2] PTR_FREE, (*par).rlat, (*par).rlon, (*par).lats, (*par).lons, (*par).latc, (*par).lonc, (*par).SolElev, (*par).datinf (*par).rlat = PTR_NEW(rlat) & (*par).lats = PTR_NEW(lats) & (*par).latc = PTR_NEW(latc) (*par).rlon = PTR_NEW(rlon) & (*par).lons = PTR_NEW(lons) & (*par).lonc = PTR_NEW(lonc) (*par).datinf = PTR_NEW(REFORM(Sinfos[idat, *])) & (*par).SolElev = PTR_NEW(ElevSol(*par)) (*par).Nseg = Nseg & (*par).iseg = 0 RETURN, 1 END ; ***** END OF FileSelect ***** ;======================================================================================= ; ***** GetGlasFile : FUNCTION ***** ; Choix d'un fichier GLAS et récupération des informations générales ; pstate (I) pointeur sur la zone contenant la structure où renseigner ces informations ;--------------------------------------------------------------------------------------- FUNCTION GetGlasDir, fpath l0 = 0 & p = fpath WHILE STRLEN(p) GE 5 DO BEGIN print,p plast = STRPOS(STRMID(p, 0, STRLEN(p) - l0 - 1), PATH_SEP(), /REVERSE_SEARCH) part = STRMID(p, plast + 1) & p = STRMID(p, 0, plast) IF STRLEN(part) EQ 5 AND STRMID(part, 0, 3) EQ "GLA" THEN BREAK ENDWHILE print,p RETURN, p + PATH_SEP() END FUNCTION GetGlasFile, pstate FORWARD_FUNCTION GetGlasDir par = (*pstate).wstrParam filein = DIALOG_PICKFILE(FILTER='GLA*.DAT', /MUST_EXIST, TITLE='GLAS File Selection', GET_PATH=pathin, PATH=(*par).wpathin) IF filein EQ '' THEN RETURN, 0 (*par).wpathin = pathin & (*par).wfilein = filein & (*par).wliddir = GetGlasDir(filein) (*par).wfileout = GetPathFile(filein, SUFFIX='.dat') ; On cherche le label et on affiche du nom de l'orbite dans la widget OrbitLabel sl = STRLEN(filein) & sdeb = STRPOS(filein, (*par).wsep, /REVERSE_SEARCH) orbit_label = STRMID(filein, sdeb + 1, sl - sdeb - 5) ;***** Pour l'instant, tant que les données niveau 2 ne sont pas toutes valides snumprd = STRMID(orbit_label, 3, 2) & numprd = FIX(snumprd) IF numprd NE 2 AND (numprd GE 11 OR numprd LT 7) THEN BEGIN Result = DIALOG_MESSAGE('Données GLA' + snumprd + ' non traitées actuellement...', /ERROR) RETURN, 0 ENDIF ;***** setTypeParms, pstate, STRMID(orbit_label, 0, 5) WIDGET_CONTROL, (*(*pstate).wstrWidg).wTypeLabel, SET_VALUE='Product : '+(*par).ftype WIDGET_CONTROL, (*(*pstate).wstrWidg).wFileLabel, SET_VALUE=orbit_label ; Ouverture du fichier et identification du Header IF (*par).lun GT 0 THEN FREE_LUN, (*par).lun OPENR, lun, (*par).wfilein, /GET_LUN, /SWAP_IF_LITTLE GetGlasFileParms, lun, recsize, nbrec, hdlth (*par).lun = lun & (*par).recsize = recsize & (*par).hdlth = hdlth & (*par).nbrec = nbrec RETURN, 1 END ; ***** END OF GetGlasFile ***** ;======================================================================================== ; ***** GetGlasFileParms : PRO ***** ; lecture des informations générales d'un fichier GLAS: longueur d'enregistrement, nombre ; d'enregistrements, nombre d'enregistrements d'entête ; * lun (I) unité de lecture du fichier ; * recsize (O) taille d'enregistrement ; * nbrec (O) nombre d'enregistrement (problème pour GLA07) ; * hdlth (O) nombre d'enregistrements de l'entête ;---------------------------------------------------------------------------------------- PRO GetGlasFileParms, lun, recsize, nbrec, hdlth meta = 'RECL= 9999999999;|NUMHEAD= 9999999999;|size_mb_ecs_data_granule=9999.99999999999999' ; Header type READU, lun, meta strs = STRSPLIT(meta, '[^0123456789.]', /EXTRACT, /REGEX) recsize = LONG(strs[0]) & hdlth = LONG(strs[1]) nbrec = LONG(DOUBLE(strs[2])*1024l*1024l + 0.5)/recsize END ; ***** END OF GetGlasFileParms ***** ;==================================================================================================== ; DEBUT DE LA FONCTION DateTimeStr ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : DateTimeStr ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne une chaîne de caractères représentant une date et une heure ; sous la forme "YYYY-MM-JJ HH:MN:SS.CCCCCC" ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; String : date et heure formattées ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; dts = DateTimeStr(jd) ; double, jd : (IN) date julienne, incluant l'heure TU à formatter ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonctions CALDAT, STRING, FIX et ROUND (IDL) ;==================================================================================================== ; VARIABLES LOCALES: aucune ;==================================================================================================== ; PRINCIPE: s.o. ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION DateTimeStr, jd CALDAT, jd, mm, dd, yy, hh, mn, ss RETURN, STRING(yy, mm, dd, hh, mn, FIX(ss), ROUND(1.E6*(ss MOD 1)), $ FORMAT='(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2,".",i6.6)') END ;==================================================================================================== ; FIN DE LA FONCTION DateTimeStr ;==================================================================================================== ;============================================================================================== ; Cette fonction recherche le chemin d'un fichier dans une arborescence ;---------------------------------------------------------------------------------------------- FUNCTION SEARCHPATH, basedir, fname FORWARD_FUNCTION SearchMacPath, SearchUnxPath, SearchWinPath CASE !VERSION.OS_FAMILY OF 'MacOS' : RETURN, SearchMacPath(basedir, fname) 'unix' : RETURN, SearchUnxPath(basedir, fname) 'Windows': RETURN, SearchWinPath(basedir, fname) ENDCASE END ;---------------------------------------------------------------------------------------------- FUNCTION SearchMacPath, basedir, fname Result = FINDFILE(basedir+fname, COUNT=pnb) IF pnb EQ 1 THEN return, Result[0] Result = FINDFILE(basedir+'*', COUNT=pnb) FOR i = 2, pnb-1 DO BEGIN IF NOT FILE_TEST(Result[i], /DIRECTORY) THEN CONTINUE p = SearchMacPath(Result[i], fname) IF p NE '' THEN return, p ENDFOR return, '' END ;---------------------------------------------------------------------------------------------- FUNCTION SearchUnxPath, basedir, fname Result = FINDFILE(basedir+fname, COUNT=pnb) IF pnb EQ 1 THEN return, Result[0] Result = FINDFILE(basedir, COUNT=pnb) FOR i = 0, pnb-1 DO BEGIN IF NOT (FILE_TEST(basedir+Result[i], /DIRECTORY) OR $ FILE_TEST(basedir+Result[i], /SYMLINK)) THEN CONTINUE p = SearchUnxPath(basedir+Result[i]+'/', fname) IF p NE '' THEN return, p ENDFOR return, '' END ;---------------------------------------------------------------------------------------------- FUNCTION SearchWinPath, basedir, fname Result = FINDFILE(basedir+fname, COUNT=pnb) IF pnb EQ 1 THEN return, Result[0] Result = FINDFILE(basedir+'*', COUNT=pnb) FOR i = 2, pnb-1 DO BEGIN IF NOT FILE_TEST(Result[i], /DIRECTORY) THEN CONTINUE p = SearchWinPath(Result[i], fname) IF p NE '' THEN return, p ENDFOR return, '' END ;======================================================================================= ; ***** GetPathDir : FUNCTION ***** ; Récupère la partie répertoire dans un chemin ; fpath (S) chemin du fichier ;--------------------------------------------------------------------------------------- FUNCTION GetPathDir, fpath plast = STRPOS(STRMID(fpath, 0, STRLEN(fpath) - 1), PATH_SEP(), /REVERSE_SEARCH) RETURN, STRMID(fpath, 0, plast) + PATH_SEP() END ;======================================================================================= ; ***** GetPathFile : FUNCTION ***** ; Récupère la partie fichier dans un chemin ; fpath (S) chemin du fichier ;--------------------------------------------------------------------------------------- FUNCTION GetPathFile, fpath, SUFFIX=suffix pfirst = STRPOS(fpath, PATH_SEP(), /REVERSE_SEARCH) p = STRMID(fpath, pfirst + 1) & l = STRLEN(p) IF KEYWORD_SET(SUFFIX) AND STRUPCASE(STRMID(p, l - STRLEN(suffix), STRLEN(suffix))) $ EQ STRUPCASE(suffix) THEN p = STRMID(p, 0, l - STRLEN(suffix)) RETURN, p END FUNCTION RayonTerrestreEquateur RETURN, 6378.2064D END FUNCTION RayonTerrestrePole RETURN, 6356.5838D END ;==================================================================================================== ; DEBUT DE LA FONCTION RayonTerrestre ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : RayonTerrestre ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne le rayon terrestre pour une latitude donnée. ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; double: rayon terrestre pour la latitude fournie ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; r = RayonTerrestre(dlat) ; double, dlat : (IN) latitude à laquelle le rayon terrestre est demandé ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonctions RayonTerrestrePole et RayonTerrestrePole (dans ce module) ; Fonctions SIN, COS et SQRT (IDL) ;==================================================================================================== ; VARIABLES LOCALES: aucune ;==================================================================================================== ; PRINCIPE: le rayon terrestre est calculé en assimilant la terre à un ellipsoïde de révolution ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION RayonTerrestre, dlat FORWARD_FUNCTION RayonTerrestrePole, RayonTerrestreEquateur RETURN, SQRT((RayonTerrestrePole( )*SIN(!DTOR*dlat))^2 + $ (RayonTerrestreEquateur()*COS(!DTOR*dlat))^2) END ;==================================================================================================== ; FIN DE LA FONCTION RayonTerrestre ;==================================================================================================== ;==================================================================================================== ; DEBUT DE LA FONCTION ElevationSolaire ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : ElevationSolaire ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne l'élévation solaire en fonction du jour julien, de l'heure TU ; et de la position en latitude et longitude (en radians). ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; double: rayon terrestre pour la latitude fournie ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; es = ElevationSolaire(jdat, lats, lons) ; double, jdat : (IN) date julienne, incluant l'heure TU, de détermination de l'élévation solaire ; double, lats : (IN) latitude de détermination de l'élévation solaire ; double, lons : (IN) longitude de détermination de l'élévation solaire ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonctions CALDAT, INTARR, SIN, COS et ASIN (IDL) ;==================================================================================================== ; VARIABLES LOCALES: ;==================================================================================================== ; PRINCIPE: ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION ElevationSolaire, jdat, lats, lons CALDAT, jdat, month, day, year, hour, minu, sec Ydays = 365.25 Mdays = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] Nday = INTARR(12)*0 FOR i = 1, 11 DO Nday[i] = Nday[i-1] + Mdays[i-1] ; On determine l'elevation solaire xlat = lats*!DTOR & xlon = lons*!DTOR tu = hour + minu/60. + sec/3600. & tet = 2.*!PI*(day + Nday[month - 1])/Ydays et = 0.000075 + 0.001868*COS(tet) - 0.032077*SIN(tet) - 0.014615*COS(2*tet) - 0.040849*SIN(2*tet) delta = 0.006918 - 0.399912*COS(tet) + 0.070257*SIN(tet) - 0.006758*COS(2*tet) + 0.000907*SIN(2*tet) - 0.002697*COS(3*tet) + 0.001480*SIN(3*tet) ah = (tu/12. - 1.)*!PI + xlon + et RETURN, ASIN(SIN(xlat)*SIN(delta) + COS(xlat)*COS(delta)*COS(ah)) END ;==================================================================================================== ; FIN DE LA FONCTION ElevationSolaire ;==================================================================================================== FUNCTION Jour, jdat, lats, lons FORWARD_FUNCTION ElevationSolaire RETURN, ElevationSolaire(jdat, lats, lons) GE 0 END FUNCTION PiedSolaire, jdat CALDAT, jdat, month, day, year, hour, minu, sec Ydays = 365.25 Mdays = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] Nday = INTARR(12)*0 FOR i = 1, 11 DO Nday[i] = Nday[i-1] + Mdays[i-1] tu = hour + minu/60. + sec/3600. & tet = 2.*!PI*(day + Nday[month - 1])/Ydays et = 0.000075 + 0.001868*COS(tet) - 0.032077*SIN(tet) - 0.014615*COS(2*tet) - 0.040849*SIN(2*tet) lat = 0.006918 - 0.399912*COS(tet) + 0.070257*SIN(tet) - 0.006758*COS(2*tet) + 0.000907*SIN(2*tet) - 0.002697*COS(3*tet) + 0.001480*SIN(3*tet) lon = - et - (tu/12. - 1.)*!PI RETURN, [lat, lon]*!RADEG END FUNCTION CycOrb, orb, orbcyc RETURN, [orb/orbcyc + 1, orb MOD orbcyc + 1] END ;==================================================================================================== ; DEBUT DE LA FONCTION rgeodang ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : rgeodang ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne l'angle géodésique entre 2 points sur la terre ; en radians, dans l'intervalle [0,Pi] ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; double: angle géodésique entre les deux points ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; dif = rgeodang(rlat1, rlon1, rlat2, rlon2) ; double, rlat1, rlat2 : (IN) latitudes des points, en radians dans l'intervalle [-Pi/2,Pi/2] ; double, rlon1, rlon2 : (IN) longitudes des points, en radians dans l'intervalle [0,2*Pi] ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonctions SIN, COS et ACOS (IDL) ;==================================================================================================== ; VARIABLES LOCALES: aucune ;==================================================================================================== ; PRINCIPE: s.o. ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION rgeodang, rlat1, rlon1, rlat2, rlon2 RETURN, ACOS(SIN(rlat1)*SIN(rlat2) + COS(rlat1)*COS(rlat2)*COS(rlon1 - rlon2)) END ;==================================================================================================== ; FIN DE LA FONCTION rgeodang ;==================================================================================================== ;==================================================================================================== ; DEBUT DE LA FONCTION dgeodang ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : dgeodang ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne l'angle géodésique entre 2 points sur la terre ; en degrés, dans l'intervalle [0°,180°] ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; double: angle géodésique entre les deux points ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; dif = dgeodang(dlat1, dlon1, dlat2, dlon2) ; double, dlat1, dlat2 : (IN) latitudes des points, en degrés dans l'intervalle [-90°,+90°] ; double, dlon1, dlon2 : (IN) longitudes des points, en degrés dans l'intervalle [0°,360°] ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonction rgeodang (dans ce module) ;==================================================================================================== ; VARIABLES LOCALES: aucune ;==================================================================================================== ; PRINCIPE: s.o. ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION dgeodang, dlat1, dlon1, dlat2, dlon2 FORWARD_FUNCTION rgeodang RETURN, !RADEG*rgeodang(!DTOR*dlat1, !DTOR*dlon1, !DTOR*dlat2, !DTOR*dlon2) END ;==================================================================================================== ; FIN DE LA FONCTION dgeodang ;==================================================================================================== ;==================================================================================================== ; DEBUT DE LA FONCTION dgeoddist ;==================================================================================================== ; DEBUT D'ENTETE DE FONCTION ; --------------------------------------------------------------------------------------------------- ; NOM : dgeoddist ; VERSION : 1.0 ; AUTEUR : Bruno SIX, CGTD ICARE ; CREATION : Novembre 2004 ; DESCRIPTION: Cette fonction retourne une estimation de la distance géodésique entre 2 points sur la ; terre, en km ;==================================================================================================== ; VALEUR DE RETOUR ; --------------------------------------------------------------------------------------------------- ; double: distance géodésique entre les deux points ;==================================================================================================== ; MODE D'APPEL ET PARAMETRES ; --------------------------------------------------------------------------------------------------- ; dif = dgeoddist(dlat1, dlon1, dlat2, dlon2) ; double, dlat1, dlat2 : (IN) latitudes des points, en degrés dans l'intervalle [-90°,+90°] ; double, dlon1, dlon2 : (IN) longitudes des points, en degrés dans l'intervalle [0°,360°] ;==================================================================================================== ; COMMON UTILISES: aucun ;==================================================================================================== ; ROUTINES UTILISES: ; --------------------------------------------------------------------------------------------------- ; Fonctions RayonTerrestreEquateur, RayonTerrestrePole, rgeodang (dans ce module) ;==================================================================================================== ; VARIABLES LOCALES: ; --------------------------------------------------------------------------------------------------- ; array[2] of double, l_d_RT_t : valeurs du rayon terrestre pour les latitudes des points ;==================================================================================================== ; PRINCIPE: ; --------------------------------------------------------------------------------------------------- ; On calcule l'angle géodésique entre les points, ainsi que les rayons terrestres correspondant à ; leur latitude respective. La distance géodésique est alors approchée par celle de même angle ; géodésique sur une sphère dont le rayon est la moyenne des rayons terrestres trouvés. ;==================================================================================================== ; FIN D'ENTETE DE FONCTION ;==================================================================================================== FUNCTION dgeoddist, dlat1, dlon1, dlat2, dlon2 FORWARD_FUNCTION RayonTerrestre, rgeodang ; Calcul des rayons terrestres pour les latitudes l_d_RT_t = RayonTerrestre([dlat1, dlat2]) ; On retourne une distance géodésique approchée par celle correspondant à l'angle géodésique sur ; une sphère dont le rayon est la moyenne des rayons terrestres trouvés RETURN, rgeodang(!DTOR*dlat1, !DTOR*dlon1, !DTOR*dlat2, !DTOR*dlon2)*(l_d_RT_t[0] + l_d_RT_t[1])/2. END ;==================================================================================================== ; FIN DE LA FONCTION dgeoddist ;====================================================================================================