;---------------------------------------------------------------------------- ; SCCS/s.cyclone_loc.pro 1.15 01/05/08 12:49:12 ; ; function cyclone_loc ; ; Purpose: ; Given a lat.-lon. grid of sea level pressure, fctn. finds where the ; centers of cyclones are (using a form of the Serreze (1995) and Serreze ; et al. (1997) algorithms) and returns a vector of the locations of those ; centers, in 1-D array index form. This function supports "pseudo- ; arbitrary" spacing: For the purposes of calculating local maximum, it ; is assumed that the grid is a 2-D grid where each internal point is ; surrounded by 8 points). The boundaries of the 2-D array are also ; assumed to be the boundaries of the domain. However, no other ; assumptions, including in terms of grid size and spacing, are made. ; ; If eiher your top/bottom or left/right boundaries are periodic, see the ; keyword list discussion of Lr_Periodic and Tb_Periodic below. Note that ; although these keywords are included, I have not tested whether ; specifying those keywords will properly detect cyclones at periodic ; boundaries; I have only tested whether the specification of those ; keywords will turn on or off the edge effect filter. ; ; File I/O: ; None. ; ; Entry and Exit States: ; N/A. ; ; Input Parameter: ; in_press Sea level pressure (in hPa) at locations defined by ; in_lon and in_lat. 2-D floating or double array. ; Unchanged by procedure. ; in_lon Longitude of points given by in_press (in decimal deg). ; Floating or double array. Same dimensions as in_press. ; Unchanged by procedure. ; in_lat Latitude of points given by in_press (in decimal deg). ; Floating or double array. Same dimensions as in_press. ; Unchanged by procedure. ; ; Output Parameter: ; out_loc Vector of the locations of cyclone centers on the ; input grid defined by in_press. Created. 1-D integer ; vector of array indices, in terms of the input array ; in_press. If there are no cyclone centers, out_loc is ; set to a scalar value of -1L. ; ; Keywords (optional unless otherwise noted): ; COUNT A named variable that returns the number of elements ; in out_loc. Created/overwrites by procedure. ; EDGE_DISTANCE Distance defining how from the edge is the "good" domain ; to be considered; if you're within this distance of the ; edge (and you're boundary is non-periodic), it's assumed ; that cyclone centers cannot be detected there. In ; meters. Not changed by function. Scalar. Default is ; 800e3 meters. ; LR/TB_PERIODIC If LR_PERIODIC is true, the left and right (i.e. col. at ; rows 0 and NX_P-1) bound. are assumed periodic, and the ; edge effect for IDing cyclones (i.e. that cyclones found ; near the edge are not valid) is assumed not to apply. ; If TB_PERIODIC is true, the top and bottom bound. (i.e. ; rows at col. 0 and NY_P-1) are assumed periodic. If ; neither are true (default), none of the bound. are ; assumed periodic. ; SEARCH_RAD_MAX Max. radius defining the region from a point the ; procedure searches to try and determine whether a given ; location is a low pressure center. In meters. Not ; changed by function. This can either be a scalar ; (which is applied to all locations) or a vector of ; the same size as in_press of Search_Rad_Max to use for ; each location. Default is 1200e3 meters. ; SEARCH_RAD_MIN Min. radius defining the region from a point the ; procedure searches to try and determine whether a given ; location is a low pressure center. In meters. Not ; changed by function. This can either be a scalar ; (which is applied to all locations) or a vector of the ; same size as in_press that gives Search_Rad_Min to use ; for each location. Default is 400e3 meters. This ; value is also used to test for multiple lows (see ; commenting below). ; SEARCH_RAD_NDIV Integer number of shells between Search_Rad_Min and ; Search_Rad_Max to search. Scalar. Default is 3. ; SLP_DIFF_TEST A low pressure center is identified if it is entirely ; surrounded by grid points in the region between ; Search_Rad_Min and Search_Rad_Max that are all higher ; in SLP than the point in question by a min. of ; Slp_Diff_Test. In hPa. Not changed by function. This ; can either be a scalar (which is applied to all ; locations) or a vector of the same size as in_press of ; slp_diff_test to use for each location. Default is 2 ; hPa. ; ; Revision History: ; - 25 Apr 2001: Orig. ver. by Johnny Lin, CIRES/Univ. of Colo. Email: ; air_jlin@yahoo.com. 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. ; - This routine is suitable for entire Earth's range of parameters. ; - No procedures called with _Extra keyword invoked. ; - User-written procedures called: DISTANCE, ELIM_MULT_CENTERS, ; EXTREMES_2D (by Ray Sterner, Johns Hopkins University/Applied Physics ; Laboratory). ; - 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_LOC, in_press, in_lon, in_lat $ , COUNT = out_count $ , EDGE_DISTANCE = in_edge_distance $ , LR_PERIODIC = lr_periodic $ , TB_PERIODIC = tb_periodic $ , SEARCH_RAD_MAX = in_search_rad_max $ , SEARCH_RAD_MIN = in_search_rad_min $ , SEARCH_RAD_NDIV = in_search_rad_ndiv $ , SLP_DIFF_TEST = in_slp_diff_test $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ON_ERROR, 0 FORWARD_FUNCTION DISTANCE, ELIM_MULT_CENTERS, EXTREMES_2D if (N_PARAMS() ne 3) then MESSAGE, 'error-bad param list' press = in_press ;- Protect input lon = in_lon lat = in_lat if (SIZE(press, /N_Dimensions) ne 2) then MESSAGE, 'error--bad press dims' npress = N_ELEMENTS(press) ;- Total num. elements in press if (npress ne N_ELEMENTS(lon)) then MESSAGE, 'error--bad dimension match' if (npress ne N_ELEMENTS(lat)) then MESSAGE, 'error--bad dimension match' if (N_ELEMENTS(in_search_rad_max) ge 1) then $ ;- Set max. search rad. (in search_rad_max = in_search_rad_max $ ; m): Default is 1200 km else $ search_rad_max = 1200e3 if (N_ELEMENTS(in_search_rad_min) ge 1) then $ ;- Set min. search rad. (in search_rad_min = in_search_rad_min $ ; m): Default is 400 km else $ search_rad_min = 400e3 if (N_ELEMENTS(in_search_rad_ndiv) eq 1) then $ ;- Set num search rad. search_rad_ndiv = in_search_rad_ndiv $ ; incr.: Default is 3 else $ search_rad_ndiv = 3 if (N_ELEMENTS(in_edge_distance) eq 1) then $ ;- Set dist. to edge (in edge_distance = in_edge_distance $ ; m): Default is 800 km else $ edge_distance = 800e3 if (N_ELEMENTS(in_slp_diff_test) ge 1) then $ ;- Set slp test value (in slp_diff_test = in_slp_diff_test $ ; hPa): Default is 2 hPa else $ slp_diff_test = 2.0 ; ------------ Start Cycling Through Each Point in Entire Domain ------------ for i=0L,npress-1L do begin ; ------ What Array Indices Surround Each Index for a Shell of Points ------- ; ; Variable shell_loc_for_i is a vector of the subscripts of the points that ; are within the region defined by search_rad_min and search_rad_top of ; the element i, and are not i itself. ; ; For each point in the spatial domain, we search through a number of ; shells (where search_rad_top expands outwards by search_rad_ndiv steps ; until it reaches search_rad_max). This enables more flexibility in ; finding centers of various sizes. idx = LINDGEN(npress) ;- Array indices dist_from_i = DISTANCE( REPLICATE(lon[i], npress) $ ;- Dist. of each , REPLICATE(lat[i], npress) $ ; pt. from i , lon, lat ) incr = (search_rad_max - search_rad_min) $ ;- Make array of / FLOAT(search_rad_ndiv) ; the lower lim. search_rad_top = ( (FINDGEN(search_rad_ndiv) + 1.0) $ ; of the search * incr ) + search_rad_min ; shell for ndiv=0L,search_rad_ndiv-1L do begin ;- Cycle through each search_rad ; division (begin) shell_loc_for_i = $ WHERE( (dist_from_i le search_rad_top[ndiv]) and $ (dist_from_i ge search_rad_min) and $ (idx ne i) $ , count ) if (count eq 0) then $ PRINT, '*** warning--domain may be too spread out ***' ; --------------- Find Locations That Pass the Low Pressure Test ------------ ; ; Method: For each location, check that the pressure of all the points in ; the shell around i, defined by search_rad_top and search_rad_min, is ; slp_diff_test higher. If so, and the shell of points around that location ; is ge 4 (which is a test to help make sure the location isn't being ; examined on the basis of just a few points), then that location is labeled ; as passing the low pressure test. ; ; Note that since the shell is based upon distance which is based on lat/lon, ; this low pressure test automatically accomodates for periodic bound., if ; the bounds are periodic. For non-periodic bounds, some edge points may ; pass this test, and thus must be removed later on in the edge effects ; removal section. ndim = SIZE( shell_loc_for_i $ ;+ test if shell empty: ndim=0 , /N_Dimensions ) ; means shell is empty if (ndim ne 0) then begin slp_diff = press[shell_loc_for_i] - press[i] ;+ diff. in slp npts_shell = N_ELEMENTS(shell_loc_for_i) ;+ num pts. surrounding i tmp = WHERE(slp_diff ge slp_diff_test, count) ;+ identify loc. if ( (count eq npts_shell) and $ ; that pass the low (npts_shell ge 4) ) then begin ; pressure test if (N_ELEMENTS(tmp_loc) eq 0) then $ tmp_loc = [i] $ else $ tmp_loc = [tmp_loc, i] endif endif endfor ;- Cycle through each search_rad ; division (end) ; ------------- End Cycling Through Each Point in Entire Domain ------------- endfor ; ----------------- Identify Low Pressure Centers Candidates ---------------- ; ; Method: From the locations that pass the SLP difference test, we find ; which ones could be low pressure centers by finding the locations that ; are local minimums in SLP. Note low_loc values are in units of indices ; of the orig. pressure array. if (N_ELEMENTS(tmp_loc) eq 0) then begin ;- Case no pts. at all meet the out_count = 0L ; low pressure test out_loc = -1L RETURN, out_loc endif tmp_loc = tmp_loc[UNIQ( tmp_loc, SORT(tmp_loc) )] ;- Keep only unique values ; of pts. that meet the ; low pressure test test_slp = press ;- Set all points in array that are test_slp[*,*] = 100000.0 ; not possible low press. candidates test_slp[tmp_loc] = press[tmp_loc] ; to 100,000 hPa. low_loc = EXTREMES_2D(test_slp, -1, /Edges) ;- Local min.: all boundaries ; assumed periodic in search ; ----- Test For Multiple Systems In a Region Defined By Search_Rad_Min ----- ; ; Method: If two low centers identified in low_loc are less than ; Search_Rad_Min apart, it is assumed that they are not actually ; separate systems, and the value with the lowest SLP value is ; retained as describing the true low center. if (SIZE(low_loc, /N_Dimensions) eq 1) then begin ;- If low_loc has ID'd ; some low centers test_slp_ll = test_slp[low_loc] lon_ll = lon[low_loc] lat_ll = lat[low_loc] emc_loc = ELIM_MULT_CENTERS( test_slp_ll, lon_ll, lat_ll $ , Type = -1 $ , Search_Rad = search_rad_min ) out_loc = low_loc[emc_loc] endif else begin ;- Case where there are no ; low_loc ID's clusters: out_count = 0L ; output function -1 out_loc = -1L RETURN, out_loc endelse ; --------------------------- Eliminate Edge Points ------------------------- ; ; Method: Eliminate all out_loc candidate points that are a distance ; Edge_Distance away from the edge, for the boundaries that are non- ; periodic. sd_p = SIZE(press) ;- Dimensions of 2-D press array nx_p = sd_p[1] ny_p = sd_p[2] ; --> Examine only non-periodic edges lat and lon: ielim_flag = 1 ;- Flag to elim. edge: default is on (=1) case 1 of ((KEYWORD_SET(lr_periodic) eq 0) and (KEYWORD_SET(tb_periodic) eq 0)): $ begin edge_lon = [ REFORM(lon[*,0]), REFORM(lon[*,ny_p-1L]) $ , REFORM(lon[0,*]), REFORM(lon[nx_p-1L,*]) ] edge_lat = [ REFORM(lat[*,0]), REFORM(lat[*,ny_p-1L]) $ , REFORM(lat[0,*]), REFORM(lat[nx_p-1L,*]) ] end ((KEYWORD_SET(lr_periodic) eq 1) and (KEYWORD_SET(tb_periodic) eq 0)): $ begin edge_lon = [ REFORM(lon[*,0]), REFORM(lon[*,ny_p-1L]) ] edge_lat = [ REFORM(lat[*,0]), REFORM(lat[*,ny_p-1L]) ] end ((KEYWORD_SET(lr_periodic) eq 0) and (KEYWORD_SET(tb_periodic) eq 1)): $ begin edge_lon = [ REFORM(lon[0,*]), REFORM(lon[nx_p-1L,*]) ] edge_lat = [ REFORM(lat[0,*]), REFORM(lat[nx_p-1L,*]) ] end ((KEYWORD_SET(lr_periodic) eq 1) and (KEYWORD_SET(tb_periodic) eq 1)): $ begin PRINT, '*** note--all bound. are per. so no near-edge pts. are elim.' ielim_flag = 0 ;+ set flag to elim. edge to off end else: MESSAGE, 'error--bad periodic keywords' endcase if (ielim_flag eq 1) then begin ;- Case elim. at least some edges ; --> Compare each low press. candidate to each edge location: n_elon = N_ELEMENTS(edge_lon) n_elat = N_ELEMENTS(edge_lat) if (n_elon ne n_elat) then MESSAGE, 'error--bad dims elon, elat' for i=0L,N_ELEMENTS(out_loc)-1L do begin dist_from_ol_i = DISTANCE( REPLICATE(lon[out_loc[i]], n_elon) $ , REPLICATE(lat[out_loc[i]], n_elat) $ , edge_lon, edge_lat ) tmp = WHERE( dist_from_ol_i le edge_distance, count ) if (count gt 0) then out_loc[i] = -1L endfor ; --> Keep only those points not near edge: good_pts = WHERE(out_loc ge 0L, count) if (count gt 0) then begin out_loc = out_loc[good_pts] endif else begin out_count = 0L out_loc = -1L RETURN, out_loc endelse endif ; ---------------------------- Clean-Up and Output -------------------------- out_loc = out_loc[SORT(out_loc)] out_count = N_ELEMENTS(out_loc) RETURN, out_loc END ; ===== end of function ===== ; ========== end file ==========