;---------------------------------------------------------------------------- ; SCCS/s.cyclone_compare.pro 1.9 01/05/10 17:07:03 ; ; function cyclone_compare ; ; Purpose: ; Given cyclones at two timesteps (Chart 1 and Chart 2, to use notation ; from Serreze [1995]), compares to find whether the Chart 2 cyclones are ; new or continuations of a Chart 1 cyclone, and what Chart 1 cyclones ; die. ; ; Although this procedure technically "stands alone," it's designed for ; use with the function CYCLONE_TRACKS. For instance, all input variable ; names are from CYCLONE_TRACKS. Thus, much of the discussion in ; CYCLONE_TRACKS about the algorithm and variable names used is also not ; repeated here. ; ; NB: 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_struct Anon. structure of input that is unchanged by procedure. ; Structure is created by ALL_VAR_STRUCT, which creates a ; structure containing all the variables that are defined ; at the time ALL_VAR_STRUCT is called. In other words, ; in_struct is a copy of all variables defined in ; CYCLONE_TRACKS at the time CYCLONE_COMPARE is called. ; ; Output Parameter: ; out_struct Anon. structure of output that is created by the function. ; The tags have the same names as the counterparts in ; CYCLONE_TRACKS, since this output structure is designed to ; update these specific variables that were read in as ; in_struct: ; out_struct = { all_tracks:all_tracks $ ; , dead_cyc:dead_cyc $ ; , names_t:names_t $ ; , i_cyclone:i_cyclone } ; ; Keywords (optional unless otherwise noted): ; None. ; ; 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. ; ; 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. ; - No procedures called with _Extra keyword invoked. ; - User-written procedures called: CMREPLICATE (by Craig B. Markwardt), ; CONVERT_INDEX (by Bob Yantosca and Martin Schultz). ; - 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_COMPARE, in_struct $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ON_ERROR, 0 FORWARD_FUNCTION CMREPLICATE, CONVERT_INDEX if (N_PARAMS() ne 1) then MESSAGE, 'error-bad param list' input = in_struct ;- Protect input parameters tags = TAG_NAMES(input) ;- Extract all variables from for itag=0L,N_ELEMENTS(tags)-1L do begin ; input structure command_str = tags[itag] + '=' $ + 'input.' + tags[itag] exec_flag = EXECUTE(command_str) if (exec_flag ne 1) then $ MESSAGE, 'error--bad execute' endfor if (N_ELEMENTS(maxlong) ne 0) then $ ;- Est. maxlong, a temp. variable MESSAGE, 'error--duplicate maxlong' $ ; used for testing mins. in a else $ ; location vector maxlong = 9999999L if ( (count_t ge maxlong) or (count_tm1 ge maxlong) or $ (count_t*count_tm1 ge maxlong) )then $ MESSAGE, 'error--too many cyclones' ; ---------------------------- Initialize Flags ----------------------------- t_paired = LONARR(count_t) ;- Initialize flag as to whether Chart 2 ; cyclones in cyc_loc have been paired to ; a Chart 1 cyclone to "no" (i.e. 0) tm1_paired = LONARR(count_tm1) ;- Initialize flag as to whether Chart 1 ; cyclones in cyc_loc_tm1 have been paired ; to a Chart 2 cyclone to "no" (i.e. 0) ; ------ (1) Check This Chart 2 Cyclone's Dist. to Each Chart 1 Cyclone ----- ; ; Method: Is the distance from the Chart 2 cyclone close enough (i.e. under ; dist_move_1delt away) to a Chart 1 cyclone that we can consider the Chart 2 ; cyclone a continuation of that Chart 1 cyclone? All Chart 1 cyclones ; meeting that condition is specified in cont_loc. If there are more than ; two such Chart 1 cyclones, the number is narrowed to 1 by choosing the ; closest one (arbitarily pick the lowest index one, if there's a tie for ; "closest"). for icl=0L,count_t-1L do begin ;- Loop Through Every Chart 2 Cyclone (begin) tmplon = REPLICATE( lon[ cyc_loc[icl] ], count_tm1 ) ;- Lon/lat for tmplat = REPLICATE( lat[ cyc_loc[icl] ], count_tm1 ) ; dist. test dist_from_icl = $ ;- Calc. dist. DISTANCE( tmplon, tmplat $ , lon[cyc_loc_tm1], lat[cyc_loc_tm1] ) cont_loc = WHERE( dist_from_icl le dist_move_1delt, count_cl ) if (count_cl ge 2) then begin ;- Case of 2+ PRINT, '*** warning--2+ cyclones within ' $ ; Chart 1 lows + 'dist_move_1delt so choose any of the closest' ; near the icl tmp = MIN(dist_from_icl, min_loc) ; Chart 2 low cont_loc = [min_loc[0]] endif if (count_cl gt 0) then begin if (N_ELEMENTS(cont_loc) ne 1) then $ ;- Small error check MESSAGE, 'error--cont_loc bad' name = names_tm1[cont_loc] ;- Extract track for this name = name[0] ; continued Chart 1 low tags_at = TAG_NAMES(all_tracks) ; (name must be scalar) tm1_at_pt = WHERE(tags_at eq name) all_tracks.(tm1_at_pt[0])[it] = cyc_loc[icl] ;- Put Chart 2 track ; location into track ; structure names_t[icl] = name ;- Add name to names_t t_paired[icl] = 1 ;- Change flag to "Yes, this Chart 2 cyclone has ; been paired with a Chart 1 cyclone" (i.e. 1) tm1_pt = WHERE( names_tm1 eq name $ ;- Change flag to "Yes, this , count_tp ) ; Chart 1 cyclone has been if (count_tp ne 1) then $ ; paired with a Chart 2 MESSAGE, 'error--bad name' ; cyclone" (i.e. 1), and ; check also whether if (MAX(tm1_paired[tm1_pt]) gt 0) then $ ; the Chart 1 low was MESSAGE, 'error--previously paired' ; previously paired tm1_paired[tm1_pt] = 1 endif endfor ;- Loop Through Every Chart 2 Cyclone (end) ; ------------------ (2) Check Remaining Unpaired Systems ------------------- ; --> (2a) Create a matrix dist_2t dimensioned (count_tm1, count_t), which ; gives the distance between all pairs of unpaired Chart 1 and ; Chart 2 cyclones. Pairs that would involve either a previously ; paired Chart 1 or Chart 2 cyclone have their pair distance set ; to !VALUES.F_INFINITY. ; *** Fill matrices of dist_2d dim. w/ t_paired and tm1_paired: tm1_paired_2t = CMREPLICATE(tm1_paired, count_t) t_paired_2t = TRANSPOSE( CMREPLICATE(t_paired, count_tm1) ) tm1_paired_2t = REFORM(tm1_paired_2t, count_tm1, count_t) t_paired_2t = REFORM(t_paired_2t, count_tm1, count_t) ; *** Lon and lat for these t_paired vs. tm1_paired couples: lon_tm1_2t = CMREPLICATE(lon[cyc_loc_tm1], count_t) lat_tm1_2t = CMREPLICATE(lat[cyc_loc_tm1], count_t) lon_t_2t = TRANSPOSE( CMREPLICATE(lon[cyc_loc], count_tm1) ) lat_t_2t = TRANSPOSE( CMREPLICATE(lat[cyc_loc], count_tm1) ) lon_tm1_2t = REFORM(lon_tm1_2t, count_tm1, count_t) lat_tm1_2t = REFORM(lat_tm1_2t, count_tm1, count_t) lon_t_2t = REFORM(lon_t_2t, count_tm1, count_t) lat_t_2t = REFORM(lat_t_2t, count_tm1, count_t) ; *** SLP at t, t-1 and ABS(SLP tendency) for these couples: slp_tm1_2t = CMREPLICATE(slp_tm1[cyc_loc_tm1], count_t) slp_t_2t = TRANSPOSE( CMREPLICATE(slp_t[cyc_loc], count_tm1) ) slp_tm1_2t = REFORM(slp_tm1_2t, count_tm1, count_t) slp_t_2t = REFORM(slp_t_2t, count_tm1, count_t) tmp = WHERE(ABS(lat_t_2t) le 5.0, count_lz) ;- Adjust ABS(SLPT) if (count_lz eq 0) then begin ; using latit. slptend_2t = ABS(slp_t_2t - slp_tm1_2t) $ * ABS( SIN(lat_ref * !DPI / 180.d0) $ / SIN(lat_t_2t * !DPI / 180.d0) ) endif else begin PRINT, '*** warning--SLPT not adjusted for latit. because ' $ , 'cyclones are too near the equator (within 5 deg)' slptend_2t = ABS(slp_t_2t - slp_tm1_2t) endelse ; *** Fill dist_2d array, and set already paired dist. to Inf: dist_2t = DISTANCE( lon_tm1_2t, lat_tm1_2t, lon_t_2t, lat_t_2t ) paired_loc = WHERE( (tm1_paired_2t eq 1) or (t_paired_2t eq 1), count_bl ) if (count_bl gt 0) then $ dist_2t[paired_loc] = !VALUES.F_INFINITY ; --> (2b) Variable rank_dist_tm1_2t is a long matrix dimensioned same as ; dist_2t that shows the order of increasing distance of the Chart 1/ ; Chart 2 couples, for and wrt each cyc_loc_tm1 (i.e. Chart 1) loc. ; Values at previously paired Chart 1 or Chart 2 lows have ; rank_dist_tm1_2t set to maxlong. Variable rank_dist_t_2t uses the ; same idea as rank_dist_tm1_2t, but for each cyc_loc (i.e. Chart 2) ; location. rank_dist_tm1_2t = LONARR(count_tm1, count_t) rank_dist_t_2t = LONARR(count_tm1, count_t) for icltm1=0L,count_tm1-1L do begin ascending_loc = SORT( dist_2t[icltm1,*] ) rank_dist_tm1_2t[icltm1,ascending_loc] = LINDGEN(count_t) if (count_bl gt 0) then $ rank_dist_tm1_2t[paired_loc] = maxlong endfor for icl=0L,count_t-1L do begin ascending_loc = SORT( dist_2t[*,icl] ) rank_dist_t_2t[ascending_loc,icl] = LINDGEN(count_tm1) if (count_bl gt 0) then $ rank_dist_t_2t[paired_loc] = maxlong endfor ; --> (2c) Look at each Chart 1 low. Which Chart 2 lows claim this Chart 1 ; low as the Chart 2 low's min. distance? For these "candidate" ; Chart 2 lows that claim this Chart 1 as its min. dist., does the ; closest of these Chart 2 lows (to the Chart 1 low) pass the tests ; (describ. below) wrt the Chart 1 low? If so, this Chart 2 is a ; continuation. If not, this Chart 2 is new, and look at next ; closest Chart 2. Repeat with other min. dist. pairs until search ; limit, and if no match, then the Chart 1 low has died. for icltm1=0L,count_tm1-1L do begin ;- Loop Every Chart 1 Cyclone (begin) other_min = WHERE( rank_dist_t_2t[icltm1,*] eq 0, count_om) if (count_om gt 0) then begin still_loop = 1B while still_loop do begin tests_passed = 1B ;- Flag default set: "Yes, candidate Chart 2 low ; passed all tests for pairing to Chart 1 low ; *** Loc. of closest candidate Chart 2 low to this Chart 1 low: tmp_2t = REPLICATE(maxlong, count_tm1, count_t) tmp_2t[icltm1, other_min] = rank_dist_tm1_2t[icltm1, other_min] tmp = MIN(tmp_2t, closest_loc) ; *** Test 1: Candidate Chart 2 low moves < test2_dist_limit ? ; If so, passes the test: if (dist_2t[closest_loc] ge test2_dist_limit) then $ tests_passed = 0B ; *** Test 2: Candidate NH (SH) Chart 2 low moves south (north) less ; than test2_deg_so_value ? If so, passes the test: case 1 of (lat_tm1_2t[closest_loc] gt 0.0): begin ;+ N. hemis. tmp = lat_tm1_2t[closest_loc] $ - lat_t_2t[closest_loc] if (tmp ge test2_deg_so_value) then $ tests_passed = 0B end (lat_tm1_2t[closest_loc] lt 0.0): begin ;+ S. hemis. tmp = lat_t_2t[closest_loc] $ - lat_tm1_2t[closest_loc] if (tmp ge test2_deg_so_value) then $ tests_passed = 0B end (lat_tm1_2t[closest_loc] eq 0.0): begin ;+ Equator tmp = ABS(lat_t_2t[closest_loc] - 0.0) if (tmp ge test2_deg_so_value) then $ tests_passed = 0B end else: MESSAGE, 'error--bad lat' endcase ; *** Test 3: Candidate Chart 2 low SLPT (SLP tendency) has magnit. ; less than test2_slpt_value ? If so, passes the test: if (slptend_2t[closest_loc] ge test2_slpt_value) then $ tests_passed = 0B ; *** The tm1 and t indices for this Chart 1 and 2 low pair: tmploc = CONVERT_INDEX(closest_loc , [count_tm1, count_t]) if (tmploc[0] ne icltm1) then MESSAGE, 'error--bad index' icl = tmploc[1] ; *** If candidate Chart 2 passed Tests 1-3: if tests_passed then begin ;+ pair Chart 2 low to Chart 1 low ; --- Check t and tm1 not previously paired: if (t_paired[icl] ne 0) then MESSAGE, 'error-prev. paired' if (tm1_paired[icltm1] ne 0) then MESSAGE, 'error-prev. paired' ; --- Extract track for this continued Chart 1 low: name = names_tm1[icltm1] name = name[0] ;+ name must be scalar tags_at = TAG_NAMES(all_tracks) tm1_at_pt = WHERE(tags_at eq name) ; --- Fill names_t and put Chart 2 track loc. into track struct.: names_t[icl] = name all_tracks.(tm1_at_pt[0])[it] = cyc_loc[icl] ; --- Change flag to "Yes, this Chart 2 cyclone has been paired ; with a Chart 1 cyclone" (i.e. 1) and vice versa: t_paired[icl] = 1 tm1_paired[icltm1] = 1 endif else begin ;+ else, Chart 2 low is new 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 to track struct. tmp_loc[it] = cyc_loc[icl] all_tracks = CREATE_STRUCT(all_tracks, name, tmp_loc) t_paired[icl] = 1 ;* change flag to "Yes, this Chart 2 low ; has been "paired" (since it's new) endelse ; *** Finished considering candidate Chart 2 low: rank_dist_tm1_2t[closest_loc] = maxlong ; *** Test for any other candidate Chart 2 lows: tmp = MIN( rank_dist_tm1_2t[icltm1, other_min] ) if (tmp ge maxlong) then still_loop = 0B endwhile endif endfor ;- Loop Through Every Chart 1 Cyclone (end) ; --> (2d) Any remaining unpaired Chart 2 lows are assumed to be new. t_unpaired_loc = WHERE(t_paired eq 0, count_tup) if (count_tup gt 0) then begin for itup=0L,count_tup-1L do begin icl = t_unpaired_loc[itup] 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 to track struct. tmp_loc[it] = cyc_loc[icl] all_tracks = CREATE_STRUCT(all_tracks, name, tmp_loc) t_paired[icl] = 1 ;* change flag to "Yes, this Chart 2 low ; has been "paired" (since it's new) endfor endif ; --> (2e) Any remaining unpaired Chart 1 lows are assumed to have died. dead_loc = WHERE(tm1_paired eq 0, count_dl) if (count_dl gt 0) then begin dead_cyc = [ dead_cyc, names_tm1[dead_loc] ] tm1_paired[dead_loc] = 1 endif ; --------------- Test All Chart 1 and 2 Lows Have Been Paired -------------- if ( (MIN(tm1_paired) ne 1) or (MIN(t_paired) ne 1) ) then $ MESSAGE, 'error--procedure failed somehow' ; ---------------------------- Clean-Up and Output -------------------------- out_struct = { all_tracks:all_tracks $ , dead_cyc:dead_cyc $ , names_t:names_t $ , i_cyclone:i_cyclone } RETURN, out_struct END ; ===== end of function ===== ; ========== end file ==========