;+ ; NAME: ; SWARM_PCP_FIND_PATCHES ; ; PURPOSE: ; This function identifies polar cap patches within a given set of plasma density data; typically one day's worth. It begins by calculating the background plasma density as a moving percentile over a certain length, followed by the calculation of the foreground plasma density, again as a moving percentile over a certain length. It then loops through all values above 77 degrees MLAT to find those regions, where the foreground is twice as big (or larger) than the background. From this regions, the polar cap patch proper, it finds the edges of each patch. ; ; CALLING SEQUENCE: ; result = swarm_pcp_find_patches( mlat, ndens, backfilterlen, backfilterperc, forefilterlen, forefilterperc ) ; ; INPUTS: ; mlat: The geomagnetic latitudes of the satellite. When run at UiO, we use the AACGM coordinates. For n data point inside the file this is a vector of length n. ; ndens: The plasma density provided by the Langmuir probe of the EFI, from the EFI data files. For n data point inside the file this is a vector of length n. ; backfilterlen: The length (in number of data points) of the percentile filter over which to calculate the background plasma density. At UiO, we use 551 points. ; backfilterperc: The percentile within the window of length backfilterlen which is defined a the background plasma density. At UiO, we use the 35th percentile. ; forefilterlen: The length (in number of data points) of the percentile filter over which to calculate the foreground plasma density. At UiO, we use 7 points. ; forefilterperc: The percentile within the window of length forefilterlen which is defined a the foreground plasma density. At UiO, we use the 50th percentile. ; ; OUTPUTS: ; pcp_flag: Each plasma density measurement is associated with this value, indicating whether the particular measurement is part of a polar cap patch, or not. The function sets the value of this flag, if inside a patch, to 2; it is 0 otherwise. ; ; OPTIONAL INPUTS: ; None. ; ; OPTIONAL OUTPUTS: ; None. ; ; KEYWORD PARAMETERS: ; pcount: Set this to a named variable that will upon completion contain the number of patches found within the data. ; ps1idx: Set this to a named variable that will upon completion contain an array of indices which mark the beginning of the first edge of all patches found. It ; has a length of pcount. ; ps2idx: Set this to a named variable that will upon completion contain an array of indices which mark the beginning of the patch proper of all patches found. It ; has a length of pcount. ; pe2idx: Set this to a named variable that will upon completion contain an array of indices which mark the end of the patch proper of all patches found. It ; has a length of pcount. ; pe1idx: Set this to a named variable that will upon completion contain an array of indices which mark the end of the second edge of all patches found. It ; has a length of pcount. ; background: Set this to a named variable that will upon completion contain the background plasma density calculated through percentile filtering. ; foreground: Set this to a named variable that will upon completion contain the foreground plasma density calculated through percentile filtering. ; PMININC: Set this to a named variable that will upon completion contain the minimum relative increase of the foreground above the background to qualify as a patch - this will always be returned as 2. ; PMAXDIFF: Set this to a named variable that will upon completion contain the percent density increase that marks the beginning and end of the patch edges - this will always be returned as 30. ; PMINLEN: Set this to a named variable that will upon completion contain the minimum number of points which make up the entire patch (including edges) - this will always be returned as 15. ; PMINMLAT: Set this to a named variable that will upon completion contain the minimum geomagnetic latitude below which we do not expect patches - this will always be returned as 77. ; ; EXAMPLE: ; ; read the electric field data, i.e., the data from both the Thermal Ion Imager and the ; ; Langmuir Probe ; swarm_efi_read, adate, probe, $ ; juls=juls, glats=glats, glons=glons, alts=alts, $ ; ndens=ndens, v_sc=v_sc, v_ion=v_ion, te=te, $ ; preliminary=preliminary ; ; ; convert goegraphic latitude to geomagnetic latitude ; mlats = swarm_pcp_convert_mlat( glats, glons, alts ) ; ; ; Find patches ; pcp_flag = swarm_pcp_find_patches( mlats, ndens, 551, 35, 7, 50, $ ; pcount=pcount, ps1idx=ps1idx, ps2idx=ps2idx, pe1idx=pe1idx, pe2idx=pe2idx, $ ; background=background, foreground=foreground, $ ; PMININC=PMININC, PMAXDIFF=PMAXDIFF, PMINLEN=PMINLEN, PMINMLAT=PMINMLAT ) ; ; COPYRIGHT: ; This file was developed as part of the Swarm + Innovation: Polar Cap Products project. ; Please visit http://www.mn.uio.no/fysikk/english/research/projects/swarm/ for more information. ; PCP (Polar Cap Products) is a project funded by the European Space Agency ; (Contract No. 4000114121/15/NL/MP) in the framework of the STSE (Support To Science Element) ; Swarm + Innovation Program. ; If you have any questions, please contact Lasse Clausen (lasse.clausen@fys.uio.no) or Wojciech ; Miloch (w.j.miloch@fys.uio.no). ;- function swarm_pcp_find_patches, mlat, ndens, $ backfilterlen, backfilterperc, forefilterlen, forefilterperc, $ pcount=pcount, ps1idx=ps1idx, ps2idx=ps2idx, pe1idx=pe1idx, pe2idx=pe2idx, $ background=background, foreground=foreground, $ PMININC=PMININC, PMAXDIFF=PMAXDIFF, PMINLEN=PMINLEN, PMINMLAT=PMINMLAT ; backfilterlen: number of points over which to percentile filter for the background ; backfilterperc: percentile that defines the background ; forefilterlen: number of points over which to percentile filter for the foreground ; forefilterperc: percentile that defines the foreground ; minimum relative increase of the foreground above the background to ; qualify as a patch PMININC = 2. ; percent density increase that marks the beginning and end of the ; patch edges PMAXDIFF = 30. ; minimum number of points during which the relative increase has to be ; higher than PINCMIN PMINLEN = 15. ; minimum geomagnetic latitude below which we expect no patches PMINMLAT = 77. ; minimum geomagnetic latitude below which we expect no patche edges PEDGEMINMLAT = PMINMLAT - 2. ; number of density measurements np = n_elements(ndens) ; percentile filter density for background background = fltarr(np) for p=backfilterlen/2, np-backfilterlen/2-1 do begin tmp = ndens[p-backfilterlen/2:p+backfilterlen/2-1] & si = sort(tmp) background[p] = tmp[si[round(backfilterperc/100.*backfilterlen)]] endfor ; percentile filter density for foreground foreground = fltarr(np) for p=forefilterlen/2, np-forefilterlen/2-1 do begin tmp = ndens[p-forefilterlen/2:p+forefilterlen/2-1] & si = sort(tmp) foreground[p] = tmp[si[round(forefilterperc/100.*forefilterlen)]] endfor ; calculate the increase between the background and the ; foreground, i.e., foreground = inc*background inc = foreground/background ; array to be returned, holding the patch flag ; pcp_flag = 1: first edge ; pcp_flag = 2: patch proper ; pcp_flag = 3: second edge pcp_flag = intarr(n_elements(inc)) ; initialize some temporary variables ; plen is the current patch length in points plen = 0 ; pstarted is a boolean value equal to TRUE if inc is currently over ; PMININC pstarted = 0 ; initial array size for pidx and plength - this will increase as ; patches are found, if necessary arrlen = 5 ; pcount is number of patches found pcount = 0 ; ps1idx holds the index of the start of each identified patch, i.e., ; the point where the foreground rises 5% above the background ps1idx = lonarr(arrlen) ; ps2idx holds the index of the point where the density rises above the ; PMININC threshold, the beginning of the patch proper ps2idx = lonarr(arrlen) ; pe2idx holds the index of the last point before the density drops ; below the PMININC threhold, the end of the patch proper pe2idx = lonarr(arrlen) ; pe1idx holds the index of the end of each identified patch, i.e., the ; last point before the foreground drops below the 5% above the ; background pe1idx = lonarr(arrlen) ; loop over all relative increase values for i=backfilterlen/2L, np-backfilterlen/2L-1L do begin ; we are only interested in the polar cap, hence ignore all values ; below PMINMLAT MLAT if abs(mlat[i]) lt PMINMLAT then begin ; if we crossed the 77 degree boundary while in a potential ; patch, do not count that one in. pstarted = 0 continue endif ; if the relative increase is above PMININC if inc[i] gt PMININC then begin ; and if we are already in a patch, increase plen if pstarted then begin plen = plen + 1 ; else start the patch by setting pstarted to TRUE and plen to 1 endif else begin pstarted = 1 plen = 1 endelse ; if the relative increase is below PMININC endif else begin ; and we were inside a patch if pstarted then begin ; ; check if the length of the patch exceeds PMINLEN ; if plen gt PMINLEN then begin ; Average backgrounc density during patch proper mi = mean(background[i-plen:i-1]) ; we define ps1idx where the relative density increase rises above ; PMAXDIF ; walk backward in time to find that point bcnt = 1 while 1 do begin ; check if index in within bounds and whether edge is not too long if i-plen-bcnt eq backfilterlen/2 $ or abs(mlat[i-plen-bcnt]) lt PEDGEMINMLAT then begin break endif ; the foreground drops below PMAXDIF percent of the average ; background during the patch proper if foreground[i-plen-bcnt] lt (1.+PMAXDIFF/100.)*mi then begin break endif bcnt = bcnt + 1 endwhile ; we define the pe1idx where the relative density increase drops ; below PMAXDIF ; walk forward in time to find that point fcnt = 1 while 1 do begin ; check if index is within bounds if i+fcnt eq np-backfilterlen/2-1 $ ; or fcnt gt PMAXEDGELEN then begin or abs(mlat[i+fcnt]) lt PEDGEMINMLAT then begin break endif ; the foreground drops below PMAXDIF percent of the average ; background during the patch proper if foreground[i+fcnt] lt (1.+PMAXDIFF/100.)*mi then begin break endif fcnt = fcnt + 1 endwhile ; if the entire length of the patch exceeds PMINLEN ; and it is located above PEDGEMINMLAT ; and it is not longer than the averaging length for the background, ; and the patch does not contain NaNs ; save the patch position if plen + bcnt + fcnt gt PMINLEN $ and min(abs(mlat[i-plen-bcnt:i+fcnt-1])) ge PEDGEMINMLAT $ and plen + bcnt + fcnt lt backfilterlen $ and total(~finite(foreground[i-plen-bcnt:i+fcnt-1])) eq 0 then begin ; check whether previous patch has been eaten while ( i - plen - bcnt ) le pe1idx[pcount-1] do begin pcount -= 1 endwhile ; increase patch count pcount = pcount + 1 ; make sure the arrays holding the patch information ; are big enough if pcount gt arrlen then begin ; if not, increase their size arrlen = 2*arrlen t_ps1idx = lonarr( arrlen ) t_ps2idx = lonarr( arrlen ) t_pe1idx = lonarr( arrlen ) t_pe2idx = lonarr( arrlen ) t_ps1idx[0:pcount-2] = ps1idx t_ps2idx[0:pcount-2] = ps2idx t_pe1idx[0:pcount-2] = pe1idx t_pe2idx[0:pcount-2] = pe2idx ps1idx = t_ps1idx ps2idx = t_ps2idx pe1idx = t_pe1idx pe2idx = t_pe2idx endif ; save index of left edge where foreground rises above PMAXDIF ; percent of the average background during the patch proper ps1idx[pcount-1] = i - plen - bcnt ; save index of right edge where foreground drops above PMAXDIF ; percent of the average background during the patch proper pe1idx[pcount-1] = i + fcnt - 1 ; save indeces of the patch proper, i.e., where the foreground is ; above PMININC times the background ps2idx[pcount-1] = i - plen pe2idx[pcount-1] = i - 1 ; we cannot decide which is the trailing/leading edge, so set PCP_FLAG to 2 here ; the proper flagging of leading/trailing edge is done in SWARM_PCP_CALC_GDI pcp_flag[ ps1idx[pcount-1]:pe1idx[pcount-1] ] = 2 ; make sure we start the search for the new patch ; after the end of the last one i = i + fcnt - 1 endif endif ; we are now not inside a patch anymore pstarted = 0 endelse endfor ; check if any patches where found if pcount eq 0 then begin ps1idx = [] ps2idx = [] pe1idx = [] pe2idx = [] return, pcp_flag endif ; strip empty elements from the arrays ps1idx = ps1idx[0:pcount-1] ps2idx = ps2idx[0:pcount-1] pe1idx = pe1idx[0:pcount-1] pe2idx = pe2idx[0:pcount-1] return, pcp_flag end