;$Id: ctm_extract.pro,v 1.21 1999/03/16 22:48:46 mgs Stab $ ;------------------------------------------------------------- ;+ ; NAME: ; CTM_EXTRACT (function) ; ; PURPOSE: ; Extracts a block of data from a 3-D CTM data cube. ; ; CATEGORY: ; CTM Tools ; ; CALLING SEQUENCE: ; CTM_EXTRACT, Data, ModelInfo, GridInfo [,Keywords] ; ; INPUTS: ; DATA -> The data cube from which to extract a region. ; DATA must be 3-dimensional (e.g. lon-lat-alt, or ; lon-lat-tracer, etc). ; ; KEYWORD PARAMETERS: ; MODELINFO -> The model info structure returned from function ; CTM_TYPE. ; ; GRIDINFO -> The grid info structure returned from function ; CTM_GRID. ; ; AVERAGE -> Bit flag indicating the dimensions over which to ; average the data: ; 1 : longitudinal ; 2 : latitudinal ; 4 : vertical ; These values can be combined. E.g., to average over ; longitude and latitude use 3. A bit set in AVERAGE ; superseeds the corresponding bit in TOTAL (see below). ; ; TOTAL -> Bit flag indicating the dimensions over which ; to sum the data: ; 1 : longitudinal ; 2 : latitudinal ; 4 : vertical ; These values can be combined. E.g., to integrate over ; longitude and latitude use 3. A bit set in AVERAGE ; superseeds the corresponding bit in TOTAL (see above). ; ; /INDEX -> If set, will interpret LAT, LEV, and LON as index ; arrays. If not set, will interpret LAT, LEV, and LON as ; ranges (i.e. two-element vectors containing min and max values). ; ; LAT -> An index array of latitudes *OR* a two-element vector ; specifying the min and max latitudes to be included in ; the extracted data block. Default is [ -90, 90 ]. ; ; LEV -> An index array of sigma levels *OR* a two-element vector ; specifying the min and max sigma levels to be included ; in the extracted data block. Default is [ 1, GRIDINFO.LMX ]. ; ; LON -> An index array of longitudes *OR* a two-element vector ; specifying the min and max longitudes to be included in ; the extracted data block. Default is [ -180, 180 ]. ; ; ALTRANGE -> A vector specifying the min and max altitude ; values to be included in the extracted data block. ; ; PRANGE -> A two-element vector specifying the min and max ; pressure levels to be included in the extracted data block. ; ; WE -> Returns to the calling program the index array of longitudes ; for the extracted data region, ordered from west to east. ; ; SN -> Returns to the calling program the index array of latitudes ; for the extracted data region, ordered from South to North. ; ; UP -> Returns to the calling program the index array of vertical ; levels for the extracted data region, ordered from surface ; to top. ; ; DIM -> A named variable will return the new dimension information ; of the data block after averaging or totaling. ; ; _EXTRA=e -> Picks up extra keywords for CTM_INDEX. ; ; OUTPUTS: ; X, Y, Z -> Arrays of latitude, longitude, or altitude values ; corresponding to the the 1st, 2nd, and 3rd dimensions of ; the DATA array, respectively. ; ; SUBROUTINES: ; CTM_INDEX ; DEFAULT_RANGE (function) ; ; REQUIREMENTS: ; Uses GAMAP package subroutines. ; ; NOTES: ; (1) CTM_EXTRACT returns the extracted data region as ; the function value. ; ; (2) Assumes a 3-D data cube as input, of dimensions (lon, lat, ; alt), or for some diagnostics (lon, lat, "tracer" number). ; ; (3) In the calling program, CTM_TYPE and CTM_GRID must be ; called to compute the MODELINFO and GRIDINFO structures, ; which can then be passed to CTM_EXTRACT. ; ; (4) If any of the LAT, LON, LEV, ALTRANGE, PRANGE keywords are ; explicity specified, then CTM_EXTRACT will return to the ; calling program their original, unmodified values. If any ; of these are not explicitly specified , then CTM_EXTRACT ; will return to the calling program default values. ; ; EXAMPLE: ; (1) ; ModelInfo = CTM_TYPE( 'GEOS1' ) ; GridInfo = CTM_GRID( ModelInfo ) ; DataRegion = CTM_EXTRACT( DataCube, ModelInfo, DataInfo, $ ; Lon=[-180,0], Lat=[-30,30], $ ; Lev=[1,10] ) ; ; ; Extracts from a GEOS-1 4x5 data cube a region ; ; delimited by longitude = [-180, 0], ; ; latitude = [-30, 30], for levels 1 to 10. ; ; (2) ; Lon = indgen( 36 ) ; Lat = indgen( 16 ) + 15 ; Lev = indgen( 10 ) ; ModelInfo = CTM_TYPE( 'GEOS1' ) ; GridInfo = CTM_GRID( ModelInfo ) ; DataRegion = CTM_EXTRACT( DataCube, ModelInfo, DataInfo, $ ; /Index, Lon=Lon, Lat=Lat, Lev=Lev ) ; ; ; Extracts same data region as in Example (1) but ; ; here passes explicit index arrays instead of ranges. ; ; MODIFICATION HISTORY: ; bmy, 16 Sep 1998: VERSION 1.00 ; bmy, 17 Sep 1998: - now extracts data from data cube one ; dimension at a time to avoid errors ; bmy, 18 Sep 1998: VERSION 1.01 ; - INDEX, SN, WE, UP keywords added ; - LATRANGE, LONRANGE, LEVRANGE renamed ; to LAT, LON, LEV (since they may now ; contain arrays and not just ranges). ; mgs, 21 Sep 1998: - some more error checking ; - removed MinData and MaxData ; bmy and mgs, 22 Sep 1998: VERSION 1.02 ; - added AVERAGE and TOTAL keywords ; bmy, 24 Sep 1998: VERSION 1.03 ; - Now returns original values of LAT, LON, ; LEV, ALTRANGE, and PRANGE if those keywords ; are specified. Otherwise returns ; defaults. ; MGS, 29 SEP 1998: - Introduced new DEFAULT_RANGE function. ; bmy, 06 Oct 1998: - fixed bug: now S = size( NEWDATA ) ; bmy, 08 Oct 1998: VERSION 1.04 ; - MODELINFO and GRIDINFO are now keywords ; - added X, Y, and Z as parameters ; bmy, 11 Feb 1999: - updated comments ; bmy, 19 Feb 1999: VERSION 1.05 ; - added FIRST keyword so that the values of ; THISDATAINFO.FIRST can be passed from the ; calling routine. ; - now call ADJ_INDEX to adjust the WE, ; SN, and UP index arrays for data blocks ; that are less than global size. ; - added DEBUG keyword ; mgs, 16 Mar 1999: - cosmetic changes ; ;- ; Copyright (C) 1998, 1999, Bob Yantosca and Martin Schultz, ; Harvard University ; This software is provided as is without any warranty ; whatsoever. It may be freely used, copied or distributed ; for non-commercial purposes. This copyright notice must be ; kept with any copy of this software. If this software shall ; be used commercially or sold as part of a larger package, ; please contact the author to arrange payment. ; Bugs and comments should be directed to bmy@io.harvard.edu ; or mgs@io.harvard.edu with subject "IDL routine ctm_extract" ;------------------------------------------------------------- function CTM_Extract, Data, X, Y, Z, $ ModelInfo=ModelInfo, GridInfo=GridInfo, $ Average=Average, Total=FTotal, $ Index=Index, Lat=Lat, $ Lev=Lev, Lon=Lon, $ AltRange=AltRange, PRange=PRange, $ WE=WE, SN=SN, $ UP=UP, Dim=Dim, $ First=First, Debug=Debug, $ _EXTRA=e ;===================================================================== ; Pass external functions ;===================================================================== FORWARD_FUNCTION Adj_Index, Default_Range ;===================================================================== ; Error checking / Keyword Settings ;===================================================================== if ( N_Elements( Data ) eq 0 ) then begin Message, 'Invalid DATA array!', /Continue return, -1 endif if ( not ChkStru( ModelInfo, [ 'NAME', 'FAMILY' ] ) ) then begin Message, 'Invalid MODELINFO structure!', /Continue return, -1 endif if ( not ChkStru( GridInfo, [ 'XEDGE', 'YEDGE' ] ) ) then begin Message, 'Invalid GRIDINFO structure!', /Continue return, -1 endif ; Flag settings for AVERAGE and TOTAL keywords if ( N_Elements( Average ) eq 0 ) then Average = 0 if ( N_Elements( FTotal ) eq 0 ) then FTotal = 0 ; Very hokey way to set FIRST. if ( N_Elements( First ) eq 0 ) then First = [0, 0, 0] Flag = ( Average OR FTotal ) ; Debugging flag Debug = Keyword_Set( Debug ) ;===================================================================== ; Default setting for INDEX ;===================================================================== Index = Keyword_Set( Index ) if ( Index ) then begin ;================================================================== ; If /INDEX is set, then interpret LON, LAT, and LEV as index ; arrays instead of ranges. ; ; NOTE: Array indices are given in FORTRAN convention. ; Hence, 1 is subtracted to convert to IDL. ; ; If array indices exceed allowed range, we end up with duplicate ; entries! Warn user... ;================================================================== WE = ((Lon - 1) > 0) < GridInfo.IMX SN = ((Lat - 1) > 0) < GridInfo.JMX UP = Lev ; conversion below Test = Uniq( WE, Sort( WE ) ) if ( N_Elements( Test ) ne N_Elements( WE ) ) then $ Message, 'Invalid array indices for LON!', /Continue Test = Uniq( SN, Sort( SN ) ) if ( N_Elements( Test ) ne N_Elements( SN ) ) then $ Message, 'Invalid array indices for LAT!', /Continue endif else begin ;================================================================== ; If INDEX = 0, then interpret LAT, LON, and LEV as two-element ; vectors containing the min and max extents of latitude, ; longitude, and vertical levels. Ditto for Altitude and Pressure. ; ; NOTE: Preserve original values of LAT, LON, LEV, ALTRANGE, and ; PRANGE, since these might be needed to create labels ; for plot titles in the calling program. ;================================================================== NewLon = Default_Range( Lon, [-182.5,180.], /NOSORT ) NewLat = Default_Range( Lat, [-90.,90.] ) NewLev = Default_Range( Lev, [1, ModelInfo.NTrop] ) NewAltRange = Default_Range( AltRange, [-1.,-1.] ) NewPRange = Default_Range( PRange, [-1.,-1.],/NOSORT ) ;================================================================== ; Call CTM index to translate Lat/Lon ranges into I/J ranges ; WE = index array of longitudinal ( West -> East ) boxes ; SN = index array of latitudinal ( South -> North ) boxes ;================================================================== CTM_Index, ModelInfo, $ Edge=[ NewLat[0], NewLon[0], NewLat[1], NewLon[1] ], $ WE_Index=WE, SN_Index=SN, /Non_Interactive, _EXTRA=e ;================================================================== ; UP = index array of sigma levels for vertical direction ;================================================================== if (ChkStru(GridInfo,'LMX')) then $ MaxLev = GridInfo.LMX $ else $ MaxLev = 1 LevInd = IndGen( MaxLev ) if ( NewAltRange[0] ge 0 ) then begin ;=============================================================== ; If ALTRANGE is specified, use the min and max altitude ; values to compute the UP index array ;=============================================================== NewAltRange = NewAltRange( Sort( NewAltRange ) ) Ind = Where( GridInfo.ZEdge ge NewAltRange[0] and $ GridInfo.ZEdge le NewAltRange[1] ) if ( Ind[0] ge 0 ) $ then UP = LevInd[ Ind ] $ else UP = LevInd[ * ] endif else if ( NewPRange[0] ge 0 ) then begin ;=============================================================== ; If PRANGE is specified, use the min and max altitude values ; to compute the UP index array ;=============================================================== NewPRange = NewPRange( Sort( NewPRange ) ) Ind = Where( GridInfo.PEdge le NewPRange[0] and $ GridInfo.PEdge ge NewPRange[1] ) if ( Ind[0] ge 0 ) $ then UP = LevInd[ Ind ] $ else UP = LevInd[ * ] endif else begin ;=============================================================== ; Otherwise, use LEV to compute the UP array. ; ; NOTE: Levels are given in FORTRAN convention: ; Convert to IDL by subtracting 1 ;=============================================================== NewLev = NewLev( Sort( NewLev ) ) UP = LevInd[ ( NewLev[0]-1 ) > 0 : ( NewLev[1]-1 ) < (MaxLev-1) ] endelse ;================================================================== ; If LAT, LON, LEV, ALTRANGE, or PRANGE were not passed ; explicitly, then return the default values to the calling ; program. Otherwise, do not modify the original values. ;================================================================== if ( N_Elements( Lat ) eq 0 ) then Lat = NewLat if ( N_Elements( Lon ) eq 0 ) then Lon = NewLon if ( N_Elements( Lev ) eq 0 ) then Lev = NewLev if ( N_Elements( AltRange ) eq 0 ) then AltRange = NewAltRange if ( N_Elements( PRange ) eq 0 ) then PRange = NewPRange endelse ;===================================================================== ; Now that we have computed the WE, SN, UP index arrays, extract the ; data region from the data cube. Proceed one dimension at a time, ; in order to avoid array subscript errors. ; ; WE, SN, and UP are "global" index arrays. Thus for data blocks ; that are smaller than global size, WE[0], SN[0], and UP[0] are ; not necessarily zero. Call ADJ_INDEX to adjust these indices ; (and store them in separate index arrays) before selecting a ; subset of the DATA array. ; ; NOTE: FIRST is in FORTRAN array notation, so we need to subtract 1. ; (but make sure it really is so ...) ;===================================================================== First = First > 1 WE_Adj = Adj_Index( WE, First[0]-1, GridInfo.IMX ) SN_Adj = Adj_Index( SN, First[1]-1, GridInfo.JMX ) UP_Adj = Adj_Index( UP, First[2]-1, MaxLev ) NewData = Data [ WE_Adj, *, * ] NewData = NewData[ *, SN_Adj, * ] NewData = NewData[ *, *, UP_Adj ] ;**** Debug output (bmy, 2/19/99) if ( Debug ) then begin print, '### CTM_EXTRACT: WE : ', WE print, '### CTM_EXTRACT: WE_Adj: ', WE_Adj print, '### CTM_EXTRACT: SN : ', SN print, '### CTM_EXTRACT: SN_Adj: ', SN_Adj print, '### CTM_EXTRACT: UP : ', UP print, '### CTM_EXTRACT: UP_Adj: ', UP_Adj endif ; get dimensions of data (lon, lat, alt, time) S = Size( NewData, /N_Dimensions ) if ( S eq 2 ) then begin ; last filter removed one dim S = size(newdata, /Dimensions) NewData = Reform( NewData, S[0], S[1], 1 ) endif S = Size( Newdata, /Dimensions ) dim = [ S, 1 ] ;===================================================================== ; Compute totals and averages over selected dimensions ; Flag = 1 : longitude, 2 : latitude, 4 : altitude ; Any combination is possible ; DIMSHIFT keeps tracks of dimensions already lost ; ; NOTE: We do not have to call CTM_ADJ_INDEX here, since WE, SN, and ; UP will be used to reference GRIDINFO.XMID, GRIDINFO.YMID, ; and GRIDINFO.ZMID, which are of global size. ;===================================================================== DimShift = 0 ; longitudinal total or average ? if ((Flag AND 1) gt 0) then begin NewData = Total(NewData,1) Dim[0] = 1 DimShift = 1 if ((Average AND 1) gt 0) then NewData = NewData / N_Elements(WE) WE = -1L endif ; latitudinal total or average ? if ((Flag AND 2) gt 0) then begin NewData = Total(NewData,2-DimShift) Dim[1] = 1 DimShift = DimShift+1 if ((Average AND 2) gt 0) then NewData = NewData / N_Elements(SN) SN = -1L endif ; altitude total or average ? if ((Flag AND 4) gt 0) then begin NewData = Total(NewData,3-DimShift) Dim[2] = 1 if ((Average AND 4) gt 0) then NewData = NewData / N_Elements(UP) UP = -1L endif ;===================================================================== ; X = index array of 1st dimension of NEWDATA ; Y = index array of 2nd dimension of NEWDATA ; Z = index array of 3rd dimension of NEWDATA ;===================================================================== X = -999 Y = -999 Z = -999 if ( WE[0] ge 0 AND N_Elements( WE ) gt 1 ) then begin X = GridInfo.Xmid[ WE ] if ( SN[0] ge 0 AND N_Elements( SN ) gt 1 ) then begin Y = GridInfo.YMid[ SN ] if ( UP[0] ge 0 AND N_Elements( UP ) gt 1 ) $ then Z = GridInfo.ZMid[ UP ] endif else begin if ( UP[0] ge 0 AND N_Elements( UP ) gt 1 ) $ then Y = GridInfo.ZMid[ UP ] endelse endif else begin if ( SN[0] ge 0 AND N_Elements( SN ) gt 1 ) then begin X = GridInfo.YMid[ SN ] if ( UP[0] ge 0 AND N_Elements( UP ) gt 1 ) $ then Y = GridInfo.ZMid[ UP ] endif else begin if ( UP[0] ge 0 AND N_Elements( UP ) gt 1 ) $ then X = GridInfo.ZMid[ UP ] endelse endelse ;===================================================================== ; Call REFORM to get rid of extra dimensions and return ;===================================================================== if ( N_Elements( NewData ) gt 1 ) then NewData = Reform( NewData ) return, NewData end