;---------------------------------------------------------------------------- ; SCCS/s.cyclone_tracks.pro 1.16 01/07/09 11:15:50 ; ; function cyclone_tracks ; ; Purpose: ; Given a timeseries of lat.-lon. grids of sea level pressure, function ; returns a structure that describes where every cyclone is located in ; space, at every moment in time, using the detection and tracking ; method of Serreze (1995) and Serreze et al. (1997). ; ; Because of the assumptions made in CYCLONE_LOC, this function supports ; "pseudo-arbitrary" spacing where the grid is a 2-D grid where each ; internal point is surrounded by 8 points, and the boundaries of the 2-D ; array are also assumed to be the boundaries of the domain. No other ; such assumptions, including in terms of grid size and spacing, are ; made. ; ; NB: Default values of all keywords except Count_Cyc are set assuming ; a timestep of 6 hrs. The values have been translated accordingly from ; Serreze (1995, 1997), which used a timestep of 12 hrs. Chart 1 lows ; are at timestep t-1; Chart 2 lows are at timestep t. ; ; File I/O: ; None. ; ; Entry and Exit States: ; N/A. ; ; Input Parameter: ; in_slp Sea level pressure (in hPa) at locations defined ; by in_lon and in_lat, at NT timesteps. 3-D ; floating or double array dimensioned (NX,NY,NT). ; If keyword Assoc_Slp is used (for cases where the ; SLP timeseries is very long), value of in_slp is ; ignored, though it still must be included. ; Unchanged by procedure. ; in_lon Longitude of points given by in_slp (in decimal ; deg). Floating or double array. Dimensioned ; (NX,NY), as in_slp is dimensioned in space. ; Unchanged by procedure. ; in_lat Latitude of points given by in_slp (in decimal ; deg). Floating or double array. Dimensioned ; (NX,NY), as in_slp is dimensioned in space. ; Unchanged by procedure. ; ; Output Parameter: ; out_tracks Created structure where each tag is a structure of ; info. about each cyclone that is ever present in ; the dataset for the NT timesteps. Variable ; out_tracks is a named structure of the form: ; {tracks, C0:LONARR(NT), C1:LONARR(NT) ... } ; where tag C0 means "cyclone number 0". The tag is ; set to a long integer array of dimension NT that ; gives the 1-D subscript index that locates where ; the center of the cyclone is on a spatial grid of ; dimension (NX,NY). At timesteps where the cyclone ; doesn't exist, the value of C0 at that point is -1L. ; If there are no cyclones in in_slp, out_tracks is ; a scalar set to -1L. ; ; Keywords (optional unless otherwise noted): ; ASSOC_SLP If the SLP timeseries is very long (e.g. more than ; one month of 6 hrly data), it is not reasonable to ; store the SLP data in memory. Instead, the SLP ; variable should be associated with a file that ; contains the data, which is what this keyword is ; for. This keyword is set to a named structure with ; tags (and variable types): ; { assoc_slp, fnassoc:STRING ; , NX:LONG, NY:LONG, NT:LONG ; , type:LONG } ; where fnassoc is the filename the SLP timeseries is ; in, NX and NY are the dimensions same as the longit. ; and latit. arrays, NT is the number of timesteps, and ; type is the IDL data type code (e.g. as from SIZE) ; of the SLP data that was written out. The SLP data ; is assumed to have been written to fnassoc in a way ; similar to the following (note in the example below, ; variable orig_data is the SLP data dimensioned as ; (NX,NY,NT)): ; OPENW, iu, fnassoc, /Get_Lun ; for it=0,NT-1 do begin ; tmpdata = REFORM(orig_data[*,*,it]) ; WRITEU, iu, tmpdata ; endfor ; FREE_LUN, iu ; The Assoc_Slp keyword structure is unchanged by the ; procedure. All values in the structure are scalars. ; COUNT_CYC A named long integer scalar variable that returns ; the number of cyclones in out_tracks. If there ; are no cyclones in the entire timeseries, then ; out_tracks=-1, Count_Cyc=0. Created and/or ; overwritten by procedure. ; DIST_MOVE_1DELT Search radius in which it is assumed if a cyclone ; at time t-1 is seen within this distance at time t, ; the cyclone is the same one. In meters. Scalar. ; Unchanged by proced. Default is 200e3 meters. ; SAVE_EACH_DELT If this keyword is set to true, at the end of each ; timestep, some variables are written out to the ; default file specified by SAVE. ; TEST2_DEG_SO_VALUE In the secondary test (i.e. the test that is ; conducted after using Dist_Move_1Delt, in Serreze ; (1995), for pairing Chart 1 and 2 cyclones) one of ; the conditions for deciding whether a Chart 2 ; cyclone is a continuation of a Chart 1 cyclone is ; based on whether it has moved Test2_Deg_So_Value ; degrees south (north) in the NH (SH). In deg. ; Scalar. Unchanged by proced. Default is 5 deg. ; TEST2_DIST_LIMIT The maximum distance the secondary test (described ; above) for pairing Chart 1 and 2 cyclones will ; search before deciding whether the Chart 1 cyclone ; died over the timestep. In meters. Scalar. ; Unchanged by proced. Default is 600e3 meters. ; TEST2_SLPT_VALUE One of the tests applied to help decide whether a ; Chart 2 cyclone is new or a continuation of a Chart 1 ; cyclone is through looking at SLP tendency between ; the two. Slp_Tend_Test is the test value used for ; this test, in units of hPa from timestep t-1 to t. ; Scalar. Default is set to 20 hPa. ; ; Revision History: ; - 8 May 2001: Orig. ver. by Johnny Lin, CIRES/Univ. of Colo. Email: ; air_jlin@yahoo.com. Passed passably reasonable tests. ; - 10 May 2001: Added in latitude correction (see Serreze 1995, p. 8) for ; sea-level pressure tendency, which was accidentally left out. ; - 17 May 2001: Assoc_Slp keyword added. Passed passably reasonable ; tests. ; ; Notes: ; - Written for IDL 5.2.1. ; - Copyright (c) 2001 Johnny Lin. For licensing and contact information ; see http://www.johnny-lin.com/lib.html. ; - Here I use "low" and "cyclone" interchangably, even though technically ; they have different meanings. In both cases, I'm trying to describe ; synoptic-scale low pressure centers. ; - Procedure should work with lows in both the N. and S. hemispheres. ; - Procedures called with _Extra keyword invoked: CYCLONE_LOC. Most of ; the keywords in that procedure are likely candidates for invoking via ; the calling line for this current proced. ; - User-written procedures called: ALL_VARS_TO_STRUCT, CYCLONE_COMPARE, ; CYCLONE_LOC. ; - No common blocks are used in this program. ; - References: ; * Serreze, M. C. (1995), Climatological aspects of cyclone development ; and decay in the Arctic, Atmos.-Oc., v. 33, pp. 1-23; ; * Serreze, M. C., F. Carse, R. G. Barry, and J. C. Rogers (1997), ; Icelandic low cyclone activity: Climatological features, linkages ; with the NAO, and relationships with recent changes in the Northern ; Hemisphere circulation, J. Clim., v. 10, pp. 453-464. ;---------------------------------------------------------------------------- FUNCTION CYCLONE_TRACKS, in_slp, in_lon, in_lat $ , ASSOC_SLP = in_assoc_slp $ , COUNT_CYC = out_count $ , DIST_MOVE_1DELT = in_dist_move_1delt $ , SAVE_EACH_DELT = save_each_delt $ , TEST2_DEG_SO_VALUE = in_test2_deg_so_value $ , TEST2_DIST_LIMIT = in_test2_dist_limit $ , TEST2_SLPT_VALUE = in_test2_slpt_value $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ON_ERROR, 0 FORWARD_FUNCTION ALL_VARS_TO_STRUCT, CYCLONE_COMPARE, CYCLONE_LOC if (N_PARAMS() ne 3) then MESSAGE, 'error-bad param list' lon = in_lon ;- Protect lon/lat input parameters lat = in_lat ; --- Set SLP array and affiliated parameters, depending on whether SLP ; is in a binary, unformatted file or not: if (N_ELEMENTS(in_assoc_slp) gt 0) then begin ;+ SLP is in file assoc_slp = in_assoc_slp OPENR, iu, assoc_slp.fnassoc, /Get_Lun NX = assoc_slp.NX NY = assoc_slp.NY NT = assoc_slp.NT case assoc_slp.type of 4: slp = ASSOC(iu, FLTARR(NX,NY)) 5: slp = ASSOC(iu, DBLARR(NX,NY)) else: MESSAGE, 'error--unrecognized type for slp' endcase endif else begin ;+ SLP is not in slp = in_slp ; file if (SIZE(slp, /N_Dimensions) ne 3) then $ MESSAGE, 'error--bad slp dims' sd_slp = SIZE(slp) ;- Dimensions of slp array NX = sd_slp[1] NY = sd_slp[2] NT = sd_slp[3] endelse if (N_ELEMENTS(lon) ne NX*NY) then MESSAGE, 'error--bad dim match' if (N_ELEMENTS(lat) ne NX*NY) then MESSAGE, 'error--bad dim match' ; --- Set other keywords: if (N_ELEMENTS(in_dist_move_1delt) eq 1) then $ ;- Set max. dist. [m] a dist_move_1delt = in_dist_move_1delt $ ; continuing cyclone else $ ; would move in 1 time- dist_move_1delt = 200e3 ; step: Default is 200 km ; --- Set secondary test (Test 2) keywords (see header comments): if (N_ELEMENTS(in_test2_deg_so_value) eq 1) then $ ;- Units deg: Default test2_deg_so_value = in_test2_deg_so_value $ ; is 5.0 deg else $ test2_deg_so_value = 5.0 if (N_ELEMENTS(in_test2_dist_limit) eq 1) then $ ;- Units m: Default test2_dist_limit = in_test2_dist_limit $ ; is 600 km else $ test2_dist_limit = 600e3 if (N_ELEMENTS(in_test2_slpt_value) eq 1) then $ ;- Units hPa: Default test2_slpt_value = in_test2_slpt_value $ ; is 20 hPa else $ test2_slpt_value = 20.0 ; ------------ Initialize a Few Special Variables and Counter(s) ------------ i_cyclone = -1L ;- Initialize cyclone number variable null_loc = REPLICATE(-1L, NT) ;- Vector used as initial value for ; timeseries of cyclone location evol. lat_ref = 60.0 ;- Ref. latitude for adjusting SLPT [deg]. This is ; arbitrary, so I've set the value to the one followed by ; Serreze (1995). NB: This is in the NH, so remember ; that fact when you use it. ; ---------------- Set Values Before Starting First Timestep ---------------- slp_tm1 = 'NA' ;- For initial timestep, set the "previous" names_tm1 = 'NA' ; timestep (i.e. Chart 1) as having no cyclones cyc_loc_tm1 = -1L count_tm1 = 0 dead_cyc = 'NA' ; -------------------- Start Loop Through All Timesteps --------------------- for it=0L,NT-1L do begin PRINT, 'Processing cyclone tracks for timestep ' $ + STRTRIM(STRING(it+1),2) ; ------------- Cyclone Locations at Timestep it (i.e. Chart 2) ------------- if (N_ELEMENTS(assoc_slp) gt 0) then $ ;+ SLP is in file slp_t = slp[it] $ else $ ;+ SLP is not in file slp_t = REFORM(slp[*,*,it]) cyc_loc = CYCLONE_LOC(slp_t, lon, lat, Count=count_t, _Extra=extra) ; ------------- Compare Chart 2 Cyclones With Chart 1 Cyclones -------------- ; ; Variable notes: ; - Each cyclone detected is appended to an anonymous structure all_tracks. ; Later, this is copied to named structure out_tracks for output. ; - dead_cyc is a vector of tag names of all cyclones that have died. At ; the beginning of the loop, it will not contain any 'NA' elements, unless ; the vector would otherwise be empty, in which case it will be a scalar ; of value 'NA'. At the end of the loop, if names_tm1 had no cyclones, ; dead_cyc will contain some 'NA' elements, and is purged of those ; elements before the loop restarts (with the exception that if the ; vector would otherwise be empty, it will contain the scalar 'NA'. ; - names_t is a string vector of the cyclone names that exist for Chart 2, ; corresponding to element-by-element the cyclones in cyc_loc. Variable ; is a scalar of value 'NA' if there are no cyclones at that timestep. ; - names_tm1 is a vector of the cyclone names that exist for Chart 1, and ; follows the same format as names_t. ; ; Method: Comparison follows this decision tree: ; - Case A: No Chart 2 cyclones exist, and ; + Case Aa: No Chart 1 cyclones exist: Thus, do nothing. ; + Case Ab: Chart 1 cyclones do exist: Thus, all Chart 1 cyclones ; are assumed to have died. ; - Case B: Chart 2 cyclones exist, and ; + Case Ba: No Chart 1 cyclones exist: Thus, all Chart 2 cyclones are ; new. ; + Case Bb: Chart 1 cyclones do exist: Thus, for each Chart 2 cyclone, ; one analyzes whether the Chart 2 cyclone is a continuation ; of a Chart 1 cyclone or is new. if (count_t eq 0) then begin ;- Case A: No Chart 2 cyclones (begin) names_t = 'NA' ;+ no cyclone names for Chart 2 if (count_tm1 eq 0) then begin ;+ case Aa: no Chart 1 cyclones ;Do nothing, since both Chart 1 ;and Chart 2 have no cyclones endif else begin ;+ case Ab: Chart 1 cyclones exist dead_cyc = [dead_cyc, names_tm1] ;> add Chart 1 cyclones to list ; of all that have died endelse endif $ ;- Case A: No Chart 2 cyclones (end) else begin ;- Case B: Chart 2 cyclones exist (begin) names_t = STRARR(count_t) ;+ initialize Chart 2 names vector if (count_tm1 eq 0) then begin ;+ case Ba: no Chart 1 cyclones for icl=0L,count_t-1L do begin ;> look at every Chart 2 cyclone: i_cyclone = i_cyclone + 1L ;* name of Chart 2 name = 'C' + STRTRIM(STRING(i_cyclone),2) ; cyclone names_t[icl] = name tmp_loc = null_loc ;* add cyclone loc tmp_loc[it] = cyc_loc[icl] ; to track vector if (N_ELEMENTS(all_tracks) eq 0) then begin ;* fill structure all_tracks = CREATE_STRUCT(name, tmp_loc) endif else begin all_tracks = CREATE_STRUCT(all_tracks, name, tmp_loc) endelse endfor endif else begin ;+ case Bb: Chart 1 cyclones exist all_vars = ALL_VARS_TO_STRUCT(Exclude_Var=['slp']) ;* all curr. ; defined vars ; (sans slp) if (SIZE(all_vars, /Type) ne 8) then begin MESSAGE, 'error--bad vars dump' endif else begin result = CYCLONE_COMPARE(all_vars) ;* result from comparing all_vars = 0 ; Charts 1 and 2 endelse all_tracks = result.all_tracks ;* replace this level's dead_cyc = result.dead_cyc ; values names_t = result.names_t i_cyclone = result.i_cyclone endelse endelse ;- Case B: Chart 2 cyclones exist (end) slp_tm1 = TEMPORARY(slp_t) ;- Set Chart 2 SLP and cyclone names to names_tm1 = TEMPORARY(names_t) ; Chart 1 for next cycle cyc_loc_tm1 = TEMPORARY(cyc_loc) count_tm1 = TEMPORARY(count_t) ; - Remove 'NA' from dead_cyc vector, if elements other than 'NA' exist; ; if all are 'NA', then set to scalar value 'NA': not_na_loc = WHERE(dead_cyc ne 'NA', count_nnl) if (count_nnl gt 0) then $ dead_cyc = dead_cyc[not_na_loc] $ else $ dead_cyc = 'NA' if (KEYWORD_SET(save_each_delt) eq 1) then begin ; - Save some variab. SAVE, it, i_cyclone, all_tracks, slp_tm1 $ ; at each timestep , names_tm1, cyc_loc_tm1 $ ; processed , count_tm1, dead_cyc, /System_Variables endif ; --------------------- End Loop Through All Timesteps ---------------------- endfor ; ---------------------------- Clean-Up and Output -------------------------- if (N_ELEMENTS(assoc_slp) gt 0) then begin ;- Close SLP file, if FREE_LUN, iu ; applicable endif if (N_ELEMENTS(all_tracks) ne 0) then begin ;- Output all_tracks and out_tracks = all_tracks ; the number of cyclones out_count = N_TAGS(out_tracks) ; in all_tracks endif else begin out_tracks = -1L out_count = 0L endelse RETURN, out_tracks END ; ===== end of function ===== ; ========== end file ==========