Documentation page for AR_detector.py¶

AR detection functions.

Author: guangzhi XU (xugzhi1987@gmail.com) Update time: 2020-07-22 09:27:22.

AR_detector.applyCropIdx(slab, cropidx)¶

Cut out a bounding box from given 2d slab given corner indices

Parameters:

slab (NCVAR or ndarray) – 2D array to cut a box from.
cropidx (tuple) – (y, x) coordinate indices, output from cropMask().

Returns:

cropslab (NCVAR or ndarray) –

2D sub array cut from <slab> using: <cropidx> as boundary indices.

AR_detector.areaFilt(mask, area, min_area=None, max_area=None, zonal_cyclic=False)¶

Filter AR binary masks by region areas

Keyword Arguments:
Parameters:	mask (ndarray) – 2D binary mask with detected objects shown as 1s. area (ndarray) – 2D map showing grid cell areas in km^2.
	min_area (float or None) – if not None, minimum area to filter objects in <mask>. max_area (float or None) – if not None, maximum area to filter objects in <mask>. zonal_cyclic (bool) – if True, treat zonal boundary as cyclic.
Returns:	result (ndarray) – 2D binary mask with objects area-filtered.

AR_detector.cart2Spherical(x, y, z, shift_lon)¶

Convert Cartesian to spherical coordiantes :param x,y,z: x-, y- and z- coordiantes. :type x,y,z: float :param shift_lon: longitude to shift so the longitude dimension

starts from this number.

Returns:	result (ndrray) – 1x3 array, the columns are the lat-, lon- and r- coordinates. r- coordinates are all 1s (unit sphere).

AR_detector.cart2Wind(vs, lats, lons)¶

Convert winds in Cartesian coordinates to u,v, inverse to wind2Cart. :param vs: Cartesian representation of the horizontal

winds.

Parameters:	lats,lons (float or ndarray) – latitude and longitude coordinates corresponding to the wind vectors given by <u> and <v>.
Returns:	u,v (float or ndarray) – u- and v- components of horizontal winds.

AR_detector.checkCyclic(mask)¶

Check binary mask is zonally cyclic or not

Parameters:	mask (ndarray) – 2d binary array, with 1s denoting existance of a feature.
Returns:	result (bool) – True if feature in <mask> is zonally cyclic, False otherwise.

AR_detector.computeTheta(p1, p2)¶

Tangent line to the arc |p1-p2|

Parameters:	p1,p2 (float) – (lat,lon) coordinates
Returns:	theta (float) – unit vector tangent at point <p1> pointing at <p2> on unit sphere.

AR_detector.cropMask(mask, edge=4)¶

Cut out a bounding box around mask==1 areas

Parameters:

mask (ndarray) – 2D binary map showing the location of an AR with 1s.
edge (int) – number of pixels as edge at 4 sides.

Returns:

mask[y1 –

y2,x1:x2] (ndarray): a sub region cut from <mask> surrouding: regions with value=1.
(yy,xx): y-, x- indices of the box of the cut region. Can later by: used in applyCropIdx(new_slab, (yy,xx)) to crop out the same region from a new array <new_slab>.

AR_detector.crossSectionFlux(mask, quslabNV, qvslabNV, axis_rdp)¶

Compute setion-wise orientation differences and cross-section fluxes in an AR

Parameters:

mask (ndarray) – CROPPED (see cropMask and applyCropIdx) 2D binary map showing the location of an AR with 1s.
quslab (NCVAR) – CROPPED (n * m) 2D array of u-flux, in kg/m/s.
qvslab (NCVAR) – CROPPED (n * m) 2D array of v-flux, in kg/m/s.
axis_rdp (ndarray) – Nx2 array storing the (lat, lon) coordinates of rdp-simplified AR axis.

Returns:

angles (ndarray) –

2D map with the same shape as <mask>, showing: section-wise orientation differences between horizontal flux (as in <quslab>, <qvslab>) and the AR axis of that section. In degrees. Regions outside of AR (0s in <mask>) are masked.

anglesmean (float): area-weighted averaged of <angles> inside <mask>. crossflux (ndarray): 2D map with the same shape as <mask>,

the section-wise cross-section fluxes in the AR, defined as the projection of fluxes onto the AR axis, i.e. flux multiplied by the cos of <angles>.

seg_thetas (list): list of (x, y, z) Cartesian coordinates of the: tangent vectors along section boundaries.

AR_detector.cyclicLabel(mask, connectivity=1, iszonalcyclic=False)¶

Label connected region, zonally cyclic version

Parameters:

mask (ndarray) – 2d binary mask with 1s denoting feature regions.

Keyword Arguments:

connectivity (int) – 1 or 2 connectivity. 2 probaby won’t work for zonally cyclic labelling.
iszonalcyclic (bool) – doing zonally cyclic labelling or not. If False, call skimage.measure.label(). If True, call skimage.measure.label() for initial labelling, then shift the map zonally by half the zonal length, do another measure.label() to find zonally linked regions. Then use this to update the original labels.

Returns:

result (ndarray) –

2d array with same shape as <mask>. Each connected: region in <mask> is given an int label.

AR_detector.determineThresLow(anoslab, sill=0.8)¶

Determine the threshold for anomalous IVT field, experimental

Keyword Arguments:
Parameters:	anoslab (ndarray) – (n * m) 2D anomalous IVT slab, in kg/m/s.
	sill (float) – float in (0, 1). Fraction of max score to define as the 1st time the score is regarded as reaching stable level.
Returns:	result (float) – determined lower threshold.

Method of determining the threshold:

1. make a loop through an array of thresholds from 1 to the 99th percentile of <ivtano>. 2. at each level, record the number of pixels > threshold. 3. after looping, pixel number counts will have a curve with a

lower-left elbow. Compute the product of number counts and thresholds P. P will have a peak value around the elbow.

choose the 1st time P reaches the 80% of max(P), and pick the corresponding threshold as result.

Largely empirical, but seems to work good on MERRA2 IVT results.

AR_detector.findARAxis(quslab, qvslab, armask_list, costhetas, sinthetas, param_dict, verbose=True)¶

Find AR axis

Parameters:

quslab (ndarray) – (n * m) 2D u-flux slab, in kg/m/s.
qvslab (ndarray) – (n * m) 2D v-flux slab, in kg/m/s.
armask_list (list) – list of 2D binary masks, each with the same shape as <quslab> etc., and with 1s denoting the location of a found AR.
costhetas (ndarray) – (n * m) 2D slab of grid cell shape: cos=dx/sqrt(dx^2+dy^2).
sinthetas (ndarray) – (n * m) 2D slab of grid cell shape: sin=dy/sqrt(dx^2+dy^2).
param_dict (dict) – a dict containing parameters controlling the detection process. Keys of the dict: ‘thres_low’, ‘min_area’, ‘max_area’, ‘max_isoq’, ‘max_isoq_hard’, ‘min_lat’, ‘max_lat’, ‘min_length’, ‘min_length_hard’, ‘rdp_thres’, ‘fill_radius’, ‘single_dome’, ‘max_ph_ratio’, ‘edge_eps’. See the doc string of findARs() for more.

Keyword Arguments:

verbose (bool) – print some messages or not.

Returns:

axes (list) –

list of AR axis coordinates. Each coordinate is defined: as a Nx2 ndarray storing (y, x) indices of the axis (indices defined in the matrix of corresponding mask in <armask_list>.)
axismask (ndarray): 2D binary mask showing all axes in <axes> merged: into one map.

New in v2.0.

AR_detector.findARs(ivt, ivtrec, ivtano, qu, qv, lats, lons, times=None, ref_time='days since 1900-01-01', thres_low=1, min_area=500000.0, max_area=18000000.0, min_LW=2, min_lat=20, max_lat=80, min_length=2000, min_length_hard=1500, rdp_thres=2, fill_radius=None, single_dome=False, max_ph_ratio=0.6, edge_eps=0.4, zonal_cyclic=False, verbose=True)¶

Find ARs from THR results, get all results in one go.

Parameters:

ivt (ndarray) – 3D or 4D input IVT data, with dimensions (time, lat, lon) or (time, level, lat, lon).
ivtrec (ndarray) – 3D or 4D array, the reconstruction component from the THR process.
ivtano (ndarray) – 3D or 4D array, the difference between input <ivt> and <ivtrec>.
qu (ndarray) – 3D or 4D array, zonal component of integrated moisture flux.
qv (ndarray) – 3D or 4D array, meridional component of integrated moisture flux.
lats (ndarray) – 1D, latitude coordinates, the length needs to be the same as the lat dimension of <ivt>.
lons (ndarray) – 1D, longitude coordinates, the length needs to be the same as the lon dimension of <ivt>.

Keyword Arguments:

times (list or array) – time stamps of the input data as a list of strings, e.g. [‘2007-01-01 06:00:00’, ‘2007-01-01 12:00’, …]. Needs to have the same length as the time dimension of <ivt>. If None, default to create a dummy 6-hourly time axis, using <ref_time> as start, with a length as the time dimension of <ivt>.
ref_time (str) – reference time point to create dummy time axis, if no time stamps are given in <times>.
thres_low (float or None) –
kg/m/s, define AR candidates as regions >= this anomalous ivt level. If None is given, compute a threshold based on anomalous ivt data in <ivtano> using an empirical method:
1. make a loop through an array of thresholds from 1 to the 99th percentile of <ivtano>. 2. at each level, record the number of pixels > threshold. 3. after looping, pixel number counts will have a curve with a

lower-left elbow. Compute the product of number counts and thresholds P. P will have a peak value around the elbow.
1. choose the 1st time P reaches the 80% of max(P), and pick the corresponding threshold as result.
min_area (float) – km^2, drop AR candidates smaller than this area.
max_area (float) – km^2, drop AR candidates larger than this area.
min_LW (float) – minimum length/width ratio.
min_lat (float) – degree, exclude systems whose centroids are lower than this latitude. NOTE this is the absolute latitude for both NH and SH. For SH, systems with centroid latitude north of - min_lat will be excluded.
max_lat (float) – degree, exclude systems whose centroids are higher than this latitude. NOTE this is the absolute latitude for both NH and SH. For SH, systems with centroid latitude south of - max_lat will be excluded.
min_length (float) – km, ARs shorter than this length is treated as relaxed.
min_length_hard (float) – km, ARs shorter than this length is discarded.
rdp_thres (float) – degree lat/lon, error when simplifying axis using rdp algorithm.
fill_radius (int or None) –
number of grids as radius to fill small holes in AR contour. If None, computed as

max(1,int(4*0.75/RESO))

where RESO is the approximate resolution in degrees of lat/lon, estimated from <lat>, <lon>.
single_dome (bool) – do peak partition or not, used to separate systems that are merged together with an outer contour.
max_ph_ratio (float) – max prominence/height ratio of a local peak. Only used when single_dome=True
edge_eps (float) – minimal proportion of flux component in a direction to total flux to allow edge building in that direction.
zonal_cyclic (bool) – if True, treat the data as zonally cyclic (e.g. entire hemisphere or global). ARs covering regions across the longitude bounds will be correctly treated as one. If your data is not zonally cyclic, or a zonal shift of the data can put the domain of interest to the center, consider doing the shift and setting this to False, as it will save computations.
verbose (bool) – print some messages or not.

Returns:

time_idx (list) – indices of the time dimension when any AR is found. labels_allNV (NCVAR): 3D array, with dimension

(time, lat, lon). At each time slice, a unique int label is assign to each detected AR at that time, and the AR region is filled out with the label value in the (lat, lon) map.

angles_allNV (NCVAR): 3D array showing orientation: differences between AR axes and fluxes, for all ARs. In degrees.
crossfluxes_allNV (NCVAR): 3D array showing cross-: sectional fluxes in all ARs. In kg/m/s.
result_df (DataFrame): AR record table. Each row is an AR, see code: in getARData() for columns.

See also

findARs(): collect and return all results in one go.

New in v2.0.

AR_detector.getARAxis(g, quslab, qvslab, mask)¶

Find AR axis from AR region mask

Parameters:

g (networkx.DiGraph) – directed planar graph constructed from AR mask and flows. See maskToGraph().
quslab (ndarray) – 2D map of u-flux.
qvslab (ndarray) – 2D map of v-flux.
mask (ndarray) – 2D binary map showing the location of an AR with 1s.

Returns:

path (ndarray) –

Nx2 array storing the AR axis coordinate indices in: (y, x) format.
axismask (ndarray): 2D binary map with same shape as <mask>, with: grid cells corresponding to coordinates in <path> set to 1s.

AR_detector.getARData(slab, quslab, qvslab, anoslab, quano, qvano, areas, lats, lons, mask_list, axis_list, timestr, param_dict)¶

Fetch AR related data

Parameters:

slab (ndarray) – (n * m) 2D array of IVT, in kg/m/s.
quslab (ndarray) – (n * m) 2D array of u-flux, in kg/m/s.
qvslab (ndarray) – (n * m) 2D array of v-flux, in kg/m/s.
anoslab (ndarray) – (n * m) 2D array of IVT anomalies, in kg/m/s.
quano (ndarray) – (n * m) 2D array of u-flux anomalies, in kg/m/s.
qvano (ndarray) – (n * m) 2D array of v-flux anomalies, in kg/m/s.
areas (ndarray) – (n * m) 2D grid cell area slab, in km^2.
lats (ndarray) – 1d array, latitude coordinates.
lons (ndarray) – 1d array, longitude coordinates.
mask_list (list) – list of 2D binary masks, each with the same shape as <anoslab> etc., and with 1s denoting the location of a found AR.
axis_list (list) – list of AR axis coordinates. Each coordinate is defined as a Nx2 ndarray storing (y, x) indices of the axis (indices defined in the matrix of corresponding mask in <masks>.)
timestr (str) – string of time snap.
param_dict (dict) – a dict containing parameters controlling the detection process. Keys of the dict: ‘thres_low’, ‘min_area’, ‘max_area’, ‘max_isoq’, ‘max_isoq_hard’, ‘min_lat’, ‘max_lat’, ‘min_length’, ‘min_length_hard’, ‘rdp_thres’, ‘fill_radius’, ‘single_dome’, ‘max_ph_ratio’, ‘edge_eps’. See the doc string of findARs() for more.

Returns:

labelsNV (NCVAR) –

(n * m) 2D int map showing all ARs at current time.: Each AR is labeled by an int label, starting from 1. Background is filled with 0s.
anglesNV (NCVAR): (n * m) 2D map showing orientation differences: between AR axes and fluxes, for all ARs. In degrees.
crossfluxesNV (NCVAR): (n * m) 2D map showing cross- sectional fluxes: in all ARs. In kg/m/s.
df (pandas.DataFrame): AR record table. Each row is an AR, see code: below for columns.

AR_detector.getNormalVectors(point_list, idx)¶

Get the normal vector and the tagent vector to the plane dividing 2 sections along the AR axis.

Parameters:

point_list (list) – list of (lat, lon) coordinates.
idx (int) – index of the point in <point_list>, denoting the point in question.

Returns:

normi (tuple) –

the (x, y, z) Cartesian coordinate of the unit: normal vector, at the point denoted by <idx>, on the Earth surface. This is the normal vector to the plane spanned by the vector Theta and P. Where P is the vector pointing to the point in question (point_list[idx]), and Theta is the tangent vector evenly dividing the angle formed by <P,P1>, and <P,P2>. Where P1, P2 are 2 points on both side of P.
thetai (tuple): the (x, y, z) Cartesian coordinate of the tangent: vector Theta above.

AR_detector.insertCropSlab(shape, cropslab, cropidx)¶

Insert the cropped sub-array back to a larger empty slab

Parameters:

shape (tuple) – (n, m) size of the larger slab.
cropslab (ndarray) – 2D array to insert.
cropidx (tuple) – (y, x) coordinate indices, output from cropMask(), defines where <cropslab> will be inserted into.

Returns:

result (ndarray) –

2D slab with shape (n, m), an empty array with a: box at <cropidx> replaced with data from <cropslab>.

AR_detector.maskToGraph(mask, quslab, qvslab, costhetas, sinthetas, edge_eps, connectivity=2)¶

Create graph from AR mask

Parameters:

mask (ndarray) – 2D binary map showing the location of an AR with 1s.
quslab (ndarray) – 2D map of u-flux.
qvslab (ndarray) – 2D map of v-flux.
costhetas (ndarray) – (n * m) 2D slab of grid cell shape: cos=dx/sqrt(dx^2+dy^2).
sinthetas (ndarray) – (n * m) 2D slab of grid cell shape: sin=dy/sqrt(dx^2+dy^2).
edge_eps (float) – float in (0,1), minimal proportion of flux component in a direction to total flux to allow edge building in that direction. Defined in Global preamble.
connectivity (int) – 1 or 2. 4- or 8- connectivity in defining neighbor- hood relationship in a 2D square grid.

Returns:

g (networkx.DiGraph) –

directed planar graph constructed from AR mask: and flows.

AR_detector.partPeaks(cropmask, cropidx, orislab, area, min_area, max_ph_ratio, fill_radius)¶

Separate local maxima by topographical prominence, watershed version

Parameters:

cropmask (ndarray) – 2D binary array, defines regions of local maxima.
cropidx (tuple) – (y, x) coordinate indices, output from cropMask().
orislab (ndarray) – 2D array, giving magnitude/height/intensity values defining the topography.
area (ndarray) – (n * m) 2D grid cell area slab, in km^2.
min_area (float) – km^2, drop AR candidates smaller than this area.
max_ph_ratio (float) – maximum peak/height ratio. Local peaks with a peak/height ratio larger than this value is treated as an independent peak.
fill_radius (int) – number of grids as radius to further separate peaks.

Returns:

result (ndarray) –

2D binary array, similar as the input <cropmask>: but with connected peaks (if any) separated so that each connected region (with 1s) denotes an independent local maximum.

AR_detector.partPeaksOld(cropmask, cropidx, orislab, max_ph_ratio)¶

Separate local maxima by topographical prominence

Parameters:

cropmask (ndarray) – 2D binary array, defines regions of local maxima.
cropidx (tuple) – (y, x) coordinate indices, output from cropMask().
orislab (ndarray) – 2D array, giving magnitude/height/intensity values defining the topography.
max_ph_ratio (float) – maximum peak/height ratio. Local peaks with a peak/height ratio larger than this value is treated as an independent peak.

Returns:

result (ndarray) –

2D binary array, similar as the input <cropmask>: but with connected peaks (if any) separated so that each connected region (with 1s) denotes an independent local maximum.

AR_detector.plotGraph(graph, ax=None, show=True)¶

Helper func to plot the graph of an AR coordinates. For debugging. :param graph: networkx Graph obj to visualize. :type graph: networkx.Graph

Keyword Arguments:
	ax (matplotlib axis obj) – axis to plot on. If None, create a new one. show (bool) – whether to show the plot or not.

AR_detector.prepareMeta(lats, lons, times, ntime, nlat, nlon, ref_time='days since 1900-01-01', verbose=True)¶

Prepare metadata for AR detection function calls

Parameters:

lats (ndarray) – 1D, latitude coordinates, the length needs to equal <nlat>.
lons (ndarray) – 1D, longitude coordinates, the length needs to equal <nlon>.
times (list or array) – time stamps of the input data as a list of strings, e.g. [‘2007-01-01 06:00:00’, ‘2007-01-01 12:00’, …]. Needs to have the a length of <ntime>.
ntime (int) – length of the time axis, should equal the length of <times>.
nlat (int) – length of the latitude axis, should equal the length of <lats>.
nlon (int) – length of the longitude axis, should equal the length of <lons>.

Keyword Arguments:

ref_time (str) – reference time point to create time axis.
verbose (bool) – print some messages or not.

Returns:

timeax (list) – a list of datetime objs. areamap (ndarray): grid cell areas in km^2, with shape (<nlat> x <nlon>). costhetas (ndarray): ratios of dx/sqrt(dx^2 + dy^2) for all grid cells.

with shape (<nlat> x <nlon>).

sinthetas (ndarray): ratios of dy/sqrt(dx^2 + dy^2) for all grid cells.: with shape (<nlat> x <nlon>).
reso (float): (approximate) horizontal resolution in degrees of lat/lon: estimate from <lats> and <lons>.

New in v2.0.

AR_detector.save2DF(result_dict)¶

Save AR records to a pandas DataFrame

Parameters:	result_dict (dict) – key: time str in ‘yyyy-mm-dd hh:00’ value: pandas dataframe. See getARData().
Returns:	result_df (pandas.DataFrame) – AR record table containing records from multiple time steps sorted by time.

AR_detector.spherical2Cart(lat, lon)¶

Convert spherical to Cartesian coordiantes :param lat,lon: latitude and longitude coordinates. :type lat,lon: float or ndarray

Returns:	result (ndarray) – Nx3 array, the columns are the x-, y- and z- coordinates.

AR_detector.uvDecomp(u0, v0, i1, i2)¶

Decompose background-transient components of u-, v- fluxes

Parameters:

u0 (ndarray) – nd array of total u-flux.
v0 (ndarray) – nd array of total v-flux.
i1 (ndarray) – nd array of the reconstruction component of IVT.
i2 (ndarray) – nd array of the anomalous component of IVT (i2 = IVT - i1).

Returns:

u1 (ndarray) –

nd array of the u-flux component corresponding to <i1>,: i.e. the background component.
v1 (ndarray): nd array of the v-flux component corresponding to <i1>,: i.e. the background component.
u2 (ndarray): nd array of the u-flux component corresponding to <i2>,: i.e. the transient component.
v2 (ndarray): nd array of the v-flux component corresponding to <i2>,: i.e. the transient component.

AR_detector.wind2Cart(u, v, lats, lons)¶

Convert u,v winds to Cartesian, consistent with spherical2Cart.

Parameters:

u,v (float or ndarray) – u- and v- components of horizontal winds.
lats,lons (float or ndarray) – latitude and longitude coordinates corresponding to the wind vectors given by <u> and <v>.

Returns:

vs (float or ndarray) –

Cartesian representation of the horizontal: winds.