+"""
+Centrographic measures for point patterns
+
+TODO
+
+- testing
+- documentation
+
+"""
+
+__author__ = "Serge Rey sjsrey@gmail.com"
+
+__all__ = [
+ "mbr",
+ "hull",
+ "mean_center",
+ "weighted_mean_center",
+ "manhattan_median",
+ "std_distance",
+ "euclidean_median",
+ "ellipse",
+ "minimum_rotated_rectangle",
+ "minimum_bounding_rectangle",
+ "skyum",
+ "dtot",
+ "_circle",
+]
+
+
+import sys
+import numpy as np
+import warnings
+import copy
+import math
+
+from math import pi as PI
+from scipy.spatial import ConvexHull
+from libpysal.cg import get_angle_between, Ray, is_clockwise
+from scipy.spatial import distance as dist
+from scipy.optimize import minimize
+import shapely
+
+not_clockwise = lambda x: not is_clockwise(x)
+
+MAXD = sys.float_info.max
+MIND = sys.float_info.min
+
+
+
+[docs]
+def minimum_bounding_rectangle(points):
+ """
+ Find minimum bounding rectangle of a point array.
+
+ Parameters
+ ----------
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ min_x : float
+ leftmost value of the vertices of minimum bounding rectangle.
+ min_y : float
+ downmost value of the vertices of minimum bounding rectangle.
+ max_x : float
+ rightmost value of the vertices of minimum bounding rectangle.
+ max_y : float
+ upmost value of the vertices of minimum bounding rectangle.
+
+ """
+ points = np.asarray(points)
+ min_x = min_y = MAXD
+ max_x = max_y = MIND
+ x, y = zip(*points)
+ min_x = min(x)
+ min_y = min(y)
+ max_x = max(x)
+ max_y = max(y)
+ return min_x, min_y, max_x, max_y
+
+
+
+def minimum_rotated_rectangle(points, return_angle=False):
+ """
+ Compute the minimum rotated rectangle for an input point set.
+
+ This is the smallest enclosing rectangle (possibly rotated)
+ for the input point set. It is computed using Shapely.
+
+ Parameters
+ ----------
+ points : numpy.ndarray
+ A numpy array of shape (n_observations, 2) containing the point
+ locations to compute the rotated rectangle
+ return_angle : bool
+ whether to return the angle (in degrees) of the angle between
+ the horizontal axis of the rectanle and the first side (i.e. length).
+
+ Returns
+ -------
+ an numpy.ndarray of shape (4, 2) containing the coordinates
+ of the minimum rotated rectangle. If return_angle is True,
+ also return the angle (in degrees) of the rotated rectangle.
+
+ """
+ points = np.asarray(points)
+ out_points = shapely.get_coordinates(
+ shapely.minimum_rotated_rectangle(shapely.multipoints(points))
+ )[:-1]
+ if return_angle:
+ angle = (
+ math.degrees(
+ math.atan2(
+ out_points[1][1] - out_points[0][1],
+ out_points[1][0] - out_points[0][0],
+ )
+ )
+ % 90
+ )
+ return (out_points[::-1], angle)
+ return out_points[::-1]
+
+
+def mbr(points):
+ warnings.warn(
+ "This function will be deprecated in the next release of pointpats.",
+ FutureWarning,
+ stacklevel=2,
+ )
+ return minimum_bounding_rectangle(points)
+
+
+mbr.__doc__ = minimum_bounding_rectangle.__doc__
+
+
+
+[docs]
+def hull(points):
+ """
+ Find convex hull of a point array.
+
+ Parameters
+ ----------
+ points: arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : array
+ (h,2), points defining the hull in counterclockwise order.
+ """
+
+ points = np.asarray(points)
+ h = ConvexHull(points)
+ return points[h.vertices]
+
+
+
+
+[docs]
+def mean_center(points):
+ """
+ Find mean center of a point array.
+
+ Parameters
+ ----------
+ points: arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : array
+ (2,), (x,y) coordinates of the mean center.
+ """
+
+ points = np.asarray(points)
+ return points.mean(axis=0)
+
+
+
+
+[docs]
+def weighted_mean_center(points, weights):
+ """
+ Find weighted mean center of a marked point pattern.
+
+ Parameters
+ ----------
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+ weights : arraylike
+ a series of attribute values of length n.
+
+ Returns
+ -------
+ _ : array
+ (2,), (x,y) coordinates of the weighted mean center.
+ """
+
+ points, weights = np.asarray(points), np.asarray(weights)
+ w = weights * 1.0 / weights.sum()
+ w.shape = (1, len(points))
+ return np.dot(w, points)[0]
+
+
+
+
+[docs]
+def manhattan_median(points):
+ """
+ Find manhattan median of a point array.
+
+ Parameters
+ ----------
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : array
+ (2,), (x,y) coordinates of the manhattan median.
+ """
+
+ points = np.asarray(points)
+ if not len(points) % 2:
+ s = "Manhattan Median is not unique for even point patterns."
+ warnings.warn(s)
+ return np.median(points, axis=0)
+
+
+
+
+[docs]
+def std_distance(points):
+ """
+ Calculate standard distance of a point array.
+
+ Parameters
+ ----------
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : float
+ standard distance.
+ """
+
+ points = np.asarray(points)
+ n, p = points.shape
+ m = points.mean(axis=0)
+ return np.sqrt(((points * points).sum(axis=0) / n - m * m).sum())
+
+
+
+
+[docs]
+def ellipse(points):
+ """
+ Calculate parameters of standard deviational ellipse for a point pattern.
+
+ Parameters
+ ----------
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : float
+ semi-major axis.
+ _ : float
+ semi-minor axis.
+ theta : float
+ clockwise rotation angle of the ellipse.
+
+ Notes
+ -----
+ Implements approach from:
+
+ https://www.icpsr.umich.edu/CrimeStat/files/CrimeStatChapter.4.pdf
+ """
+
+ points = np.asarray(points)
+ n, k = points.shape
+ x = points[:, 0]
+ y = points[:, 1]
+ xd = x - x.mean()
+ yd = y - y.mean()
+ xss = (xd * xd).sum()
+ yss = (yd * yd).sum()
+ cv = (xd * yd).sum()
+ num = (xss - yss) + np.sqrt((xss - yss) ** 2 + 4 * (cv) ** 2)
+ den = 2 * cv
+ theta = np.arctan(num / den)
+ cos_theta = np.cos(theta)
+ sin_theta = np.sin(theta)
+ n_2 = n - 2
+ sd_x = (2 * (xd * cos_theta - yd * sin_theta) ** 2).sum() / n_2
+ sd_y = (2 * (xd * sin_theta - yd * cos_theta) ** 2).sum() / n_2
+ return np.sqrt(sd_x), np.sqrt(sd_y), theta
+
+
+
+
+[docs]
+def dtot(coord, points):
+ """
+ Sum of Euclidean distances between event points and a selected point.
+
+ Parameters
+ ----------
+ coord : arraylike
+ (x,y) coordinates of a point.
+ points : arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ d : float
+ sum of Euclidean distances.
+
+ """
+ points = np.asarray(points)
+ xd = points[:, 0] - coord[0]
+ yd = points[:, 1] - coord[1]
+ d = np.sqrt(xd * xd + yd * yd).sum()
+ return d
+
+
+
+
+[docs]
+def euclidean_median(points):
+ """
+ Calculate the Euclidean median for a point pattern.
+
+ Parameters
+ ----------
+ points: arraylike
+ (n,2), (x,y) coordinates of a series of event points.
+
+ Returns
+ -------
+ _ : array
+ (2,), (x,y) coordinates of the Euclidean median.
+
+ """
+ points = np.asarray(points)
+ start = mean_center(points)
+ res = minimize(dtot, start, args=(points,))
+ return res["x"]
+
+
+
+def minimum_bounding_circle(points):
+ """
+ Implements Skyum (1990)'s algorithm for the minimum bounding circle in R^2.
+
+ Store points clockwise.
+ Find p in S that maximizes angle(prec(p), p, succ(p) THEN radius(prec(
+ p), p, succ(p)). This is also called the lexicographic maximum, and is the last
+ entry of a list of (radius, angle) in lexicographical order.
+
+ * If angle(prec(p), p, succ(p)) <= 90 degrees, then finish.
+
+ * If not, remove p from set.
+
+ Parameters
+ ----------
+ points : numpy.ndarray
+ a numpy array of shape (n_observations, 2) to compute
+ the minimum bounding circle
+
+ Returns
+ -------
+ (x,y),center for the minimum bounding circle.
+ """
+ points = hull(points)
+ if not_clockwise(points):
+ points = points[::-1]
+ if not_clockwise(points):
+ raise Exception("Points are neither clockwise nor counterclockwise")
+ POINTS = copy.deepcopy(points)
+ removed = []
+ i = 0
+ if HAS_NUMBA:
+ circ = _skyum_numba(POINTS)[0]
+ else:
+ circ = _skyum_lists(POINTS)[0]
+ return (circ[1], circ[2]), circ[0]
+
+
+
+[docs]
+def skyum(points):
+ warnings.warn(
+ "This function will be deprecated in the next release of pointpats.",
+ FutureWarning,
+ stacklevel=2,
+ )
+ return minimum_bounding_circle(points)
+
+
+
+skyum.__doc__ = (
+ "WARNING: This function is deprecated in favor of minimum_bounding_circle\n"
+ + minimum_bounding_circle.__doc__
+)
+
+
+def _skyum_lists(points):
+ points = points.tolist()
+ removed = []
+ i = 0
+ while True:
+ angles = [
+ _angle(
+ _prec(p, points),
+ p,
+ _succ(p, points),
+ )
+ for p in points
+ ]
+ circles = [
+ _circle(
+ _prec(p, points),
+ p,
+ _succ(p, points),
+ )
+ for p in points
+ ]
+ radii = [c[0] for c in circles]
+ lexord = np.lexsort((radii, angles)) # confusing as hell defaults...
+ lexmax = lexord[-1]
+ candidate = (
+ _prec(points[lexmax], points),
+ points[lexmax],
+ _succ(points[lexmax], points),
+ )
+ if angles[lexmax] <= (np.pi / 2.0):
+ # print("Constrained by points: {}".format(candidate))
+ return _circle(*candidate), points, removed, candidate
+ else:
+ try:
+ removed.append((points.pop(lexmax), i))
+ except IndexError:
+ raise Exception("Construction of Minimum Bounding Circle failed!")
+ i += 1
+
+
+try:
+ from numba import njit, boolean
+
+ HAS_NUMBA = True
+
+ @njit(fastmath=True)
+ def _skyum_numba(points):
+ i = 0
+ complete = False
+ while not complete:
+ complete, points, candidate, circle = _skyum_iteration(points)
+ if complete:
+ return circle, points, None, candidate
+
+ @njit(fastmath=True)
+ def _skyum_iteration(points):
+ points = points.reshape(-1, 2)
+ n = points.shape[0]
+ angles = np.empty((n,))
+ circles = np.empty((n, 3))
+ for i in range(n):
+ p = points[(i - 1) % n]
+ q = points[i % n]
+ r = points[(i + 1) % n]
+ angles[i] = _angle(p, q, r)
+ circles[i] = _circle(p, q, r)
+ radii = circles[:, 0]
+ # workaround for no lexsort in numba
+ angle_argmax = angles.argmax()
+ angle_max = angles[angle_argmax]
+ # the maximum radius for the largest angle
+ lexmax = (radii * (angles == angle_max)).argmax()
+
+ candidate = (lexmax - 1) % n, lexmax, (lexmax + 1) % n
+ if angles[lexmax] <= (np.pi / 2.0):
+ return True, points, lexmax, circles[lexmax]
+ else:
+ mask = np.ones((n,), dtype=boolean)
+ mask[lexmax] = False
+ new_points = points[mask, :]
+ return False, new_points, lexmax, circles[lexmax]
+
+
+except ModuleNotFoundError:
+
+ def njit(func, **kwargs):
+ return func
+
+
+@njit
+def _angle(p, q, r):
+ pq = p - q
+ rq = r - q
+ magnitudes = np.linalg.norm(pq) * np.linalg.norm(rq)
+ return np.abs(np.arccos(np.dot(pq, rq) / magnitudes))
+
+
+def _prec(p, l):
+ """
+ retrieve the predecessor of p in list l
+ """
+ pos = l.index(p)
+ if pos - 1 < 0:
+ return l[-1]
+ else:
+ return l[pos - 1]
+
+
+def _succ(p, l):
+ """
+ retrieve the successor of p in list l
+ """
+ pos = l.index(p)
+ if pos + 1 >= len(l):
+ return l[0]
+ else:
+ return l[pos + 1]
+
+
+@njit
+def _euclidean_distance(px, py, qx, qy):
+ return np.sqrt((px - qx) ** 2 + (py - qy) ** 2)
+
+
+@njit
+def _circle(p, q, r, dmetric=_euclidean_distance):
+ """
+ Returns (radius, (center_x, center_y)) of the circumscribed circle by the
+ triangle pqr.
+
+ note, this does not assume that p!=q!=r
+ """
+ px, py = p
+ qx, qy = q
+ rx, ry = r
+ angle = np.abs(_angle(p, q, r))
+ if np.abs(angle - np.pi) < 1e-5: # angle is pi
+ radius = _euclidean_distance(px, py, rx, ry) / 2.0
+ center_x = (px + rx) / 2.0
+ center_y = (py + ry) / 2.0
+ elif np.abs(angle) < 1e-5: # angle is zero
+ radius = _euclidean_distance(px, py, qx, qy) / 2.0
+ center_x = (px + qx) / 2.0
+ center_y = (py + qy) / 2.0
+ else:
+ D = 2 * (px * (qy - ry) + qx * (ry - py) + rx * (py - qy))
+ center_x = (
+ (px**2 + py**2) * (qy - ry)
+ + (qx**2 + qy**2) * (ry - py)
+ + (rx**2 + ry**2) * (py - qy)
+ ) / float(D)
+ center_y = (
+ (px**2 + py**2) * (rx - qx)
+ + (qx**2 + qy**2) * (px - rx)
+ + (rx**2 + ry**2) * (qx - px)
+ ) / float(D)
+ radius = _euclidean_distance(center_x, center_y, px, py)
+ return radius, center_x, center_y
+
+import numpy
+import warnings
+from scipy import spatial, interpolate
+from collections import namedtuple
+from .geometry import (
+ area as _area,
+ k_neighbors as _k_neighbors,
+ build_best_tree as _build_best_tree,
+ prepare_hull as _prepare_hull,
+ TREE_TYPES,
+)
+from .random import poisson
+
+
+__all__ = [
+ "f",
+ "g",
+ "k",
+ "j",
+ "l",
+ "f_test",
+ "g_test",
+ "k_test",
+ "j_test",
+ "l_test",
+]
+
+
+def _prepare(coordinates, support, distances, metric, hull, edge_correction):
+ """
+ prepare the arguments to convert into a standard format
+ 1. cast the coordinates to a numpy array
+ 2. precomputed metrics must have distances provided
+ 3. metrics must be callable or string
+ 4. warn if distances are specified and metric is not default
+ 5. make distances a numpy.ndarray
+ 6. construct the support, accepting:
+ - num_steps -> a linspace with len(support) == num_steps
+ from zero to a quarter of the bounding box's smallest side
+ - (stop, ) -> a linspace with len(support) == 20
+ from zero to stop
+ - (start, stop) -> a linspace with len(support) == 20
+ from start to stop
+ - (start, stop, num_steps) -> a linspace with len(support) == num_steps
+ from start to stop
+ - numpy.ndarray -> passed through
+ """
+ # Throw early if edge correction is requested
+ if edge_correction is not None:
+ raise NotImplementedError("Edge correction is not currently implemented.")
+
+ # cast to coordinate array
+ if isinstance(coordinates, TREE_TYPES):
+ tree = coordinates
+ coordinates = tree.data
+ else:
+ coordinates = numpy.asarray(coordinates)
+ hull = _prepare_hull(coordinates, hull)
+
+ # evaluate distances
+ if (distances is None) and metric == "precomputed":
+ raise ValueError(
+ "If metric =`precomputed` then distances must"
+ " be provided as a (n,n) numpy array."
+ )
+ if not (isinstance(metric, str) or callable(metric)):
+ raise TypeError(
+ f"`metric` argument must be callable or a string. Recieved: {metric}"
+ )
+ if distances is not None and metric != "euclidean":
+ warnings.warn(
+ "Distances were provided. The specified metric will be ignored."
+ " To use precomputed distances with a custom distance metric,"
+ " do not specify a `metric` argument.",
+ stacklevel=2,
+ )
+ metric = "euclidean"
+
+ if support is None:
+ support = 20
+
+ if isinstance(support, int): # if just n_steps, use the max nnd
+ # this is O(n log n) for kdtrees & balltrees
+ tmp_tree = _build_best_tree(coordinates, metric=metric)
+ max_dist = _k_neighbors(tmp_tree, coordinates, 1)[0].max()
+ support = numpy.linspace(0, max_dist, num=support)
+ # otherwise, we need to build it using (start, stop, step) semantics
+ elif isinstance(support, tuple):
+ if len(support) == 1: # assuming this is with zero implicit start
+ support = numpy.linspace(0, support[0], num=20) # default support n bins
+ elif len(support) == 2:
+ support = numpy.linspace(*support, num=20) # default support n bins
+ elif len(support) == 3:
+ support = numpy.linspace(support[0], support[1], num=support[2])
+ else: # try to use it as is
+ try:
+ support = numpy.asarray(support)
+ except:
+ raise TypeError(
+ "`support` must be a tuple (either (start, stop, step), (start, stop) or (stop,)),"
+ " an int describing the number of breaks to use to evalute the function,"
+ " or an iterable containing the breaks to use to evaluate the function."
+ " Recieved object of type {}: {}".format(type(support), support)
+ )
+
+ return coordinates, support, distances, metric, hull, edge_correction
+
+
+# ------------------------------------------------------------#
+# Statistical Functions #
+# ------------------------------------------------------------#
+
+
+
+[docs]
+def f(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+):
+ """
+ Ripley's F function
+
+ The so-called "empty space" function, this is the cumulative density function of
+ the distances from a random set of points to the known points in the pattern.
+
+ Parameters
+ ----------
+ coordinates : numpy.ndarray of shape (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+
+ Returns
+ -------
+ a tuple containing the support values used to evalute the function
+ and the values of the function at each distance value in the support.
+ """
+ coordinates, support, distances, metric, hull, _ = _prepare(
+ coordinates, support, distances, metric, hull, edge_correction
+ )
+ if distances is not None:
+ n = coordinates.shape[0]
+ if distances.ndim == 2:
+ k, p = distances.shape
+ if k == p == n:
+ warnings.warn(
+ f"A full distance matrix is not required for this function, and"
+ f" the intput matrix is a square {n},{n} matrix. Only the"
+ f" distances from p random points to their nearest neighbor within"
+ f" the pattern is required, as an {n},p matrix. Assuming the"
+ f" provided distance matrix has rows pertaining to input"
+ f" pattern and columns pertaining to the output points.",
+ stacklevel=2,
+ )
+ distances = distances.min(axis=0)
+ elif k == n:
+ distances = distances.min(axis=0)
+ else:
+ raise ValueError(
+ f"Distance matrix should have the same rows as the input"
+ f" coordinates with p columns, where n may be equal to p."
+ f" Recieved an {k},{p} distance matrix for {n} coordinates"
+ )
+ elif distances.ndim == 1:
+ p = len(distances)
+ else:
+ # Do 1000 empties. Users can control this by computing their own
+ # empty space distribution.
+ n_empty_points = 1000
+
+ randoms = poisson(hull=hull, size=(n_empty_points, 1))
+ try:
+ tree
+ except NameError:
+ tree = _build_best_tree(coordinates, metric)
+ finally:
+ distances, _ = tree.query(randoms, k=1)
+ distances = distances.squeeze()
+
+ counts, bins = numpy.histogram(distances, bins=support)
+ fracs = numpy.cumsum(counts) / counts.sum()
+
+ return bins, numpy.asarray([0, *fracs])
+
+
+
+
+[docs]
+def g(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ edge_correction=None,
+):
+ """
+ Ripley's G function
+
+ The G function is computed from the cumulative density function of the nearest neighbor
+ distances between points in the pattern.
+
+ Parameters
+ -----------
+ coordinates : numpy.ndarray of shape (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, n) or (n,)
+ distances from every point in the point to another point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+
+ Returns
+ -------
+ a tuple containing the support values used to evalute the function
+ and the values of the function at each distance value in the support.
+
+ """
+
+ coordinates, support, distances, metric, *_ = _prepare(
+ coordinates, support, distances, metric, None, edge_correction
+ )
+ if distances is not None:
+ if distances.ndim == 2:
+ if distances.shape[0] == distances.shape[1] == coordinates.shape[0]:
+ warnings.warn(
+ "The full distance matrix is not required for this function,"
+ " only the distance to the nearest neighbor within the pattern."
+ " Computing this and discarding the rest.",
+ stacklevel=2,
+ )
+ distances = distances.min(axis=1)
+ else:
+ k, p = distances.shape
+ n = coordinates.shape[0]
+ raise ValueError(
+ " Input distance matrix has an invalid shape: {k},{p}."
+ " Distances supplied can either be 2 dimensional"
+ " square matrices with the same number of rows"
+ " as `coordinates` ({n}) or 1 dimensional and contain"
+ " the shortest distance from each point in "
+ " `coordinates` to some other point in coordinates."
+ )
+ elif distances.ndim == 1:
+ if distances.shape[0] != coordinates.shape[0]:
+ raise ValueError(
+ f"Distances are not aligned with coordinates! Distance"
+ f" matrix must be (n_coordinates, n_coordinates), but recieved"
+ f" {distances.shape} instead of ({coordinates.shape[0]},)"
+ )
+ else:
+ raise ValueError(
+ "Distances supplied can either be 2 dimensional"
+ " square matrices with the same number of rows"
+ " as `coordinates` or 1 dimensional and contain"
+ " the shortest distance from each point in "
+ " `coordinates` to some other point in coordinates."
+ " Input matrix was {distances.ndim} dimensioanl"
+ )
+ else:
+ try:
+ tree
+ except NameError:
+ tree = _build_best_tree(coordinates, metric)
+ finally:
+ distances, indices = _k_neighbors(tree, coordinates, k=1)
+
+ counts, bins = numpy.histogram(distances.squeeze(), bins=support)
+ fracs = numpy.cumsum(counts) / counts.sum()
+
+ return bins, numpy.asarray([0, *fracs])
+
+
+
+
+[docs]
+def j(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ truncate=True,
+):
+ """
+ Ripely's J function
+
+ The so-called "spatial hazard" function, this is a function relating the F and G functions.
+
+ Parameters
+ -----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: tuple of numpy.ndarray
+ precomputed distances to use to evaluate the j function.
+ The first must be of shape (n,n) or (n,) and is used in the g function.
+ the second must be of shape (n,p) or (p,) (with p possibly equal to n)
+ used in the f function.
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern for the f function.
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ truncate: bool (default: True)
+ whether or not to truncate the results when the F function reaches one. If the
+ F function is one but the G function is less than one, this function will return
+ numpy.nan values.
+
+ Returns
+ -------
+ a tuple containing the support values used to evalute the function
+ and the values of the function at each distance value in the support.
+ """
+ if distances is not None:
+ g_distances, f_distances = distances
+ else:
+ g_distances = f_distances = None
+ fsupport, fstats = f(
+ coordinates,
+ support=support,
+ distances=f_distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ )
+
+ gsupport, gstats = g(
+ coordinates,
+ support=support,
+ distances=g_distances,
+ metric=metric,
+ edge_correction=edge_correction,
+ )
+
+ if isinstance(support, numpy.ndarray):
+ if not numpy.allclose(gsupport, support):
+ gfunction = interpolate.interp1d(gsupport, gstats, fill_value=1)
+ gstats = gfunction(support)
+ gsupport = support
+ if not (numpy.allclose(gsupport, fsupport)):
+ ffunction = interpolate.interp1d(fsupport, fstats, fill_value=1)
+ fstats = ffunction(gsupport)
+ fsupport = gsupport
+
+ with numpy.errstate(invalid="ignore", divide="ignore"):
+ hazard_ratio = (1 - gstats) / (1 - fstats)
+ both_zero = (gstats == 1) & (fstats == 1)
+ hazard_ratio[both_zero] = numpy.nan
+ if truncate:
+ result = _truncate(gsupport, hazard_ratio)
+ if len(result[1]) != len(hazard_ratio):
+ warnings.warn(
+ f"requested {support} bins to evaluate the J function, but"
+ f" it reaches infinity at d={result[0][-1]:.4f}, meaning only"
+ f" {len(result[0])} bins will be used to characterize the J function.",
+ stacklevel=2,
+ )
+ return result
+ else:
+ return gsupport, hazard_ratio
+
+
+
+
+[docs]
+def k(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ edge_correction=None,
+):
+ """
+ Ripley's K function
+
+ This function counts the number of pairs of points that are closer than a given distance.
+ As d increases, K approaches the number of point pairs.
+
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+
+ Returns
+ -------
+ a tuple containing the support values used to evalute the function
+ and the values of the function at each distance value in the support.
+ """
+ coordinates, support, distances, metric, hull, edge_correction = _prepare(
+ coordinates, support, distances, metric, None, edge_correction
+ )
+ n = coordinates.shape[0]
+ upper_tri_n = n * (n - 1) * 0.5
+ if distances is not None:
+ if distances.ndim == 1:
+ if distances.shape[0] != upper_tri_n:
+ raise ValueError(
+ f"Shape of inputted distances is not square, nor is the upper triangular"
+ f" matrix matching the number of input points. The shape of the input matrix"
+ f" is {distances.shape}, but required shape is ({upper_tri_n},) or ({n},{n})"
+ )
+ upper_tri_distances = distances
+ elif distances.shape[0] == distances.shape[1] == n:
+ upper_tri_distances = distances[numpy.triu_indices_from(distances, k=1)]
+ else:
+ raise ValueError(
+ f"Shape of inputted distances is not square, nor is the upper triangular"
+ f" matrix matching the number of input points. The shape of the input matrix"
+ f" is {distances.shape}, but required shape is ({upper_tri_n},) or ({n},{n})"
+ )
+ else:
+ upper_tri_distances = spatial.distance.pdist(coordinates, metric=metric)
+ n_pairs_less_than_d = (upper_tri_distances < support.reshape(-1, 1)).sum(axis=1)
+ intensity = n / _area(hull)
+ k_estimate = ((n_pairs_less_than_d * 2) / n) / intensity
+ return support, k_estimate
+
+
+
+
+[docs]
+def l(
+ coordinates,
+ support=None,
+ permutations=9999,
+ distances=None,
+ metric="euclidean",
+ edge_correction=None,
+ linearized=False,
+):
+ """
+ Ripley's L function
+
+ This is a scaled and shifted version of the K function that accounts for the K function's
+ increasing expected value as distances increase. This means that the L function, for a
+ completely random pattern, should be close to zero at all distance values in the support.
+
+ Parameters
+ ----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ linearized : bool
+ whether or not to subtract l from its expected value (support) at each
+ distance bin. This centers the l function on zero for all distances.
+ Proposed by Besag (1977)
+
+ Returns
+ -------
+ a tuple containing the support values used to evalute the function
+ and the values of the function at each distance value in the support.
+ """
+
+ support, k_estimate = k(
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ edge_correction=edge_correction,
+ )
+
+ l = numpy.sqrt(k_estimate / numpy.pi)
+
+ if linearized:
+ return support, l - support
+ return support, l
+
+
+
+# ------------------------------------------------------------#
+# Statistical Tests based on Ripley Functions #
+# ------------------------------------------------------------#
+
+FtestResult = namedtuple(
+ "FtestResult", ("support", "statistic", "pvalue", "simulations")
+)
+GtestResult = namedtuple(
+ "GtestResult", ("support", "statistic", "pvalue", "simulations")
+)
+JtestResult = namedtuple(
+ "JtestResult", ("support", "statistic", "pvalue", "simulations")
+)
+KtestResult = namedtuple(
+ "KtestResult", ("support", "statistic", "pvalue", "simulations")
+)
+LtestResult = namedtuple(
+ "LtestResult", ("support", "statistic", "pvalue", "simulations")
+)
+
+_ripley_dispatch = {
+ "F": (f, FtestResult),
+ "G": (g, GtestResult),
+ "J": (j, JtestResult),
+ "K": (k, KtestResult),
+ "L": (l, LtestResult),
+}
+
+
+def _ripley_test(
+ calltype,
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ keep_simulations=False,
+ n_simulations=9999,
+ **kwargs,
+):
+ stat_function, result_container = _ripley_dispatch.get(calltype)
+ core_kwargs = dict(
+ support=support,
+ metric=metric,
+ edge_correction=edge_correction,
+ )
+ tree = _build_best_tree(coordinates, metric=metric)
+ hull = _prepare_hull(coordinates, hull)
+ if calltype in ("F", "J"): # these require simulations
+ core_kwargs["hull"] = hull
+ # amortize to avoid doing this every time
+ empty_space_points = poisson(coordinates, size=(1000, 1))
+ if distances is None:
+ empty_space_distances, _ = _k_neighbors(tree, empty_space_points, k=1)
+ if calltype == "F":
+ distances = empty_space_distances.squeeze()
+ else: # calltype == 'J':
+ n_distances, _ = _k_neighbors(tree, coordinates, k=1)
+ distances = (n_distances.squeeze(), empty_space_distances.squeeze())
+ else:
+ pass
+ core_kwargs.update(**kwargs)
+
+ observed_support, observed_statistic = stat_function(
+ tree, distances=distances, **core_kwargs
+ )
+ core_kwargs["support"] = observed_support
+ n_observations = coordinates.shape[0]
+
+ if keep_simulations:
+ simulations = numpy.empty((len(observed_support), n_simulations)).T
+ pvalues = numpy.ones_like(observed_support)
+ for i_replication in range(n_simulations):
+ random_i = poisson(hull, size=n_observations)
+ if calltype in ("F", "J"):
+ random_tree = _build_best_tree(random_i, metric)
+ empty_distances, _ = random_tree.query(empty_space_points, k=1)
+ if calltype == "F":
+ core_kwargs["distances"] = empty_distances.squeeze()
+ else: # calltype == 'J':
+ n_distances, _ = _k_neighbors(random_tree, random_i, k=1)
+ core_kwargs["distances"] = (
+ n_distances.squeeze(),
+ empty_distances.squeeze(),
+ )
+ rep_support, simulations_i = stat_function(random_i, **core_kwargs)
+ pvalues += simulations_i >= observed_statistic
+ if keep_simulations:
+ simulations[i_replication] = simulations_i
+ pvalues /= n_simulations + 1
+ pvalues = numpy.minimum(pvalues, 1 - pvalues)
+ return result_container(
+ observed_support,
+ observed_statistic,
+ pvalues,
+ simulations if keep_simulations else None,
+ )
+
+
+
+[docs]
+def f_test(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ keep_simulations=False,
+ n_simulations=9999,
+):
+ """
+ Ripley's F function
+
+ The so-called "empty space" function, this is the cumulative density function of
+ the distances from a random set of points to the known points in the pattern.
+
+ When the estimated statistic is larger than simulated values at a given distance, then
+ the pattern is considered "dispersed" or "regular"
+
+ Parameters
+ -----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ keep_simulations: bool
+ whether or not to keep the simulation envelopes. If so,
+ will be returned as the result's simulations attribute
+ n_simulations: int
+ how many simulations to conduct, assuming that the reference pattern
+ has complete spatial randomness.
+
+ Returns
+ -------
+ a named tuple with properties
+ - support, the exact distance values used to evalute the statistic
+ - statistic, the values of the statistic at each distance
+ - pvalue, the percent of simulations that were as extreme as the observed value
+ - simulations, the distribution of simulated statistics (shaped (n_simulations, n_support_points))
+ or None if keep_simulations=False (which is the default)
+ """
+
+ return _ripley_test(
+ "F",
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ keep_simulations=keep_simulations,
+ n_simulations=n_simulations,
+ )
+
+
+
+
+[docs]
+def g_test(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ keep_simulations=False,
+ n_simulations=9999,
+):
+ """
+ Ripley's G function
+
+ The G function is computed from the cumulative density function of the nearest neighbor
+ distances between points in the pattern.
+
+ When the G function is below the simulated values, it suggests dispersion.
+
+ Parameters
+ ----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ keep_simulations: bool
+ whether or not to keep the simulation envelopes. If so,
+ will be returned as the result's simulations attribute
+ n_simulations: int
+ how many simulations to conduct, assuming that the reference pattern
+ has complete spatial randomness.
+
+ Returns
+ -------
+ a named tuple with properties
+ - support, the exact distance values used to evalute the statistic
+ - statistic, the values of the statistic at each distance
+ - pvalue, the percent of simulations that were as extreme as the observed value
+ - simulations, the distribution of simulated statistics (shaped (n_simulations, n_support_points))
+ or None if keep_simulations=False (which is the default)
+ """
+ return _ripley_test(
+ "G",
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ keep_simulations=keep_simulations,
+ n_simulations=n_simulations,
+ )
+
+
+
+
+[docs]
+def j_test(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ truncate=True,
+ keep_simulations=False,
+ n_simulations=9999,
+):
+ """
+ Ripley's J function
+
+ The so-called "spatial hazard" function, this is a function relating the F and G functions.
+
+ When the J function is consistently below 1, then it indicates clustering.
+ When consistently above 1, it suggests dispersion.
+
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ keep_simulations: bool
+ whether or not to keep the simulation envelopes. If so,
+ will be returned as the result's simulations attribute
+ n_simulations: int
+ how many simulations to conduct, assuming that the reference pattern
+ has complete spatial randomness.
+
+ Returns
+ -------
+ a named tuple with properties
+ - support, the exact distance values used to evalute the statistic
+ - statistic, the values of the statistic at each distance
+ - pvalue, the percent of simulations that were as extreme as the observed value
+ - simulations, the distribution of simulated statistics (shaped (n_simulations, n_support_points))
+ or None if keep_simulations=False (which is the default)
+ """
+ result = _ripley_test(
+ "J",
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ keep_simulations=keep_simulations,
+ n_simulations=n_simulations,
+ truncate=False,
+ )
+ if truncate:
+ result_trunc = _truncate(*result)
+ result_trunc = JtestResult(*result_trunc)
+ if len(result_trunc.statistic) != len(result.statistic):
+ warnings.warn(
+ f"requested {support} bins to evaluate the J function, but"
+ f" it reaches infinity at d={result[0][-1]:.4f}, meaning only"
+ f" {len(result[0])} bins will be used to characterize the J function.",
+ stacklevel=2,
+ )
+ return result_trunc
+
+ else:
+ return result
+
+
+
+
+[docs]
+def k_test(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ keep_simulations=False,
+ n_simulations=9999,
+):
+ """
+ Ripley's K function
+
+ This function counts the number of pairs of points that are closer than a given distance.
+ As d increases, K approaches the number of point pairs.
+
+ When the K function is below simulated values, it suggests that the pattern is dispersed.
+
+ Parameters
+ ----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ keep_simulations: bool
+ whether or not to keep the simulation envelopes. If so,
+ will be returned as the result's simulations attribute
+ n_simulations: int
+ how many simulations to conduct, assuming that the reference pattern
+ has complete spatial randomness.
+
+ Returns
+ -------
+ a named tuple with properties
+ - support, the exact distance values used to evalute the statistic
+ - statistic, the values of the statistic at each distance
+ - pvalue, the percent of simulations that were as extreme as the observed value
+ - simulations, the distribution of simulated statistics (shaped (n_simulations, n_support_points))
+ or None if keep_simulations=False (which is the default)
+ """
+ return _ripley_test(
+ "K",
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ keep_simulations=keep_simulations,
+ n_simulations=n_simulations,
+ )
+
+
+
+
+[docs]
+def l_test(
+ coordinates,
+ support=None,
+ distances=None,
+ metric="euclidean",
+ hull=None,
+ edge_correction=None,
+ linearized=False,
+ keep_simulations=False,
+ n_simulations=9999,
+):
+ """
+ Ripley's L function
+
+ This is a scaled and shifted version of the K function that accounts for the K function's
+ increasing expected value as distances increase. This means that the L function, for a
+ completely random pattern, should be close to zero at all distance values in the support.
+
+ When the L function is negative, this suggests dispersion.
+
+ Parameters
+ ----------
+ coordinates : numpy.ndarray, (n,2)
+ input coordinates to function
+ support : tuple of length 1, 2, or 3, int, or numpy.ndarray
+ tuple, encoding (stop,), (start, stop), or (start, stop, num)
+ int, encoding number of equally-spaced intervals
+ numpy.ndarray, used directly within numpy.histogram
+ distances: numpy.ndarray, (n, p) or (p,)
+ distances from every point in a random point set of size p
+ to some point in `coordinates`
+ metric: str or callable
+ distance metric to use when building search tree
+ hull: bounding box, scipy.spatial.ConvexHull, shapely.geometry.Polygon
+ the hull used to construct a random sample pattern, if distances is None
+ edge_correction: bool or str
+ whether or not to conduct edge correction. Not yet implemented.
+ keep_simulations: bool
+ whether or not to keep the simulation envelopes. If so,
+ will be returned as the result's simulations attribute
+ n_simulations: int
+ how many simulations to conduct, assuming that the reference pattern
+ has complete spatial randomness.
+
+ Returns
+ -------
+ a named tuple with properties
+ - support, the exact distance values used to evalute the statistic
+ - statistic, the values of the statistic at each distance
+ - pvalue, the percent of simulations that were as extreme as the observed value
+ - simulations, the distribution of simulated statistics (shaped (n_simulations, n_support_points))
+ or None if keep_simulations=False (which is the default)
+ """
+ return _ripley_test(
+ "L",
+ coordinates,
+ support=support,
+ distances=distances,
+ metric=metric,
+ hull=hull,
+ edge_correction=edge_correction,
+ linearized=linearized,
+ keep_simulations=keep_simulations,
+ n_simulations=n_simulations,
+ )
+
+
+
+def _truncate(support, realizations, *rest):
+ is_invalid = numpy.isinf(realizations) | numpy.isnan(realizations)
+ first_inv = is_invalid.argmax()
+ if not is_invalid.any():
+ return support, realizations, *rest
+ elif first_inv < len(realizations):
+ return (
+ support[:first_inv],
+ realizations[:first_inv],
+ *[r[:first_inv] if r is not None else None for r in rest],
+ )
+
+import numpy as np
+
+
+
+[docs]
+def plot_density(
+ data,
+ bandwidth,
+ kernel=None,
+ resolution=100,
+ levels=10,
+ fill=False,
+ margin=0.1,
+ ax=None,
+ figsize=None,
+ **kwargs,
+):
+ """Plot kernel density of a given point pattern
+
+ The KDE can be done either using :class:`statsmodels.nonparametric.KDEMultivariate`,
+ which is used when ``kernel=None``, or using :class:`KDEpy.FFTKDE` when kernel is
+ set. :class:`~KDEpy.FFTKDE` tends to be generally faster in most cases but may need
+ different than ``"gaussian"`` kernel to resolve in some cases. For small data of up
+ to 10 000 points, the difference is not noticeable. For larger data, specify
+ ``bandwidth`` to enforce the use of :class:`~KDEpy.FFTKDE`. Note that while being
+ faster, :class:`~KDEpy.FFTKDE` may in some case result in erroneous KDE.
+
+ KDE is plotted using matplotlib's :meth:`~matplotlib.pyplot.contour` or
+ :meth:`~matplotlib.pyplot.contourf` function to plot the density.
+
+ If MultiPoints are given, each point is treated as separate observation.
+
+ Parameters
+ ----------
+ data : array or geopandas object
+ Array with a shape (2, n) containing coordinates of points
+ or a geopandas object with (Multi)Point geometry. Assumes
+ projected coordinates, geographical coordinates (latitude, longitude)
+ are not supported.
+ bandwidth : float
+ bandwidth in the units of CRS in which data is
+ kernel : str | None, optional
+ The kernel function. If None, defaults to the Gaussian kernel and statsmodels
+ implementation. If set, uses KDEpy implementation. See
+ :meth:`KDEpy.FFTKDE._available_kernels.keys()` for choices.
+ resolution : int | tuple(int, int), optional
+ resolution of the grid used to evaluate the probability density
+ function. If tuple, each dimension of the grid is specified separately.
+ By default 100
+ levels : int or array-like, optional
+ Determines the number and positions of the contour lines / regions.
+ See the documentation of :meth:`~matplotlib.pyplot.contour` for details.
+ By default 10
+ fill : bool, optional
+ Fill the area between contour lines, by default False
+ margin : float, optional
+ The factor of the margin by which the extent of the data will be expanded when
+ creating the grid. 0.1 means 10% on each side, by default 0.1. Only used
+ with the ``statsmodels`` implementation.
+ ax : matplotlib.axes.Axes (default None)
+ axes on which to draw the plot
+ figsize : tuple of integers (default None)
+ Size of the resulting ``matplotlib.figure.Figure``. If the argument
+ ``ax`` is given explicitly, ``figsize`` is ignored.
+ **kwargs
+ Keyword arguments passed to :meth:`~matplotlib.pyplot.contour` or
+ :meth:`~matplotlib.pyplot.contourf` used for further
+ styling of the plot, for example ``cmap``, ``linewidths``, ``linestyles``,
+ or `alpha`. See the documentation of :meth:`~matplotlib.pyplot.contour` for
+ details.
+
+ Returns
+ -------
+ matplotlib.axes.Axes
+ matplotlib axes instance with the contour plot
+ """
+ if kernel is None:
+ try:
+ import statsmodels.api as sm
+ except ImportError as err:
+ raise ImportError(
+ "statsmodels is required for `plot_density` when kernel"
+ "is not specified."
+ ) from err
+
+ engine = "sm"
+ else:
+ try:
+ from KDEpy import FFTKDE
+ except ImportError as err:
+ raise ImportError(
+ "KDEpy is required for `plot_density` when kernel is not None."
+ ) from err
+
+ engine = "kdepy"
+
+ try:
+ import matplotlib.pyplot as plt
+ except ImportError as err:
+ raise ImportError("matplotlib is required for `plot_density`") from err
+
+ if ax is None:
+ _, ax = plt.subplots(figsize=figsize)
+
+ ax.set_aspect("equal") # bandwidth is fixed, hence aspect shall be equal
+
+ if isinstance(data, np.ndarray):
+ pass
+ else: # geopandas
+ if not data.geom_type.str.contains("Point").all():
+ raise ValueError(
+ "data contain non-point geometries. "
+ "Only (Multi)Points are supported."
+ )
+ data = data.get_coordinates().values
+
+ if engine == "sm":
+ dens_u = sm.nonparametric.KDEMultivariate(
+ data=[data[:, 0], data[:, 1]],
+ var_type="cc",
+ bw=[bandwidth, bandwidth],
+ )
+
+ xmax = data[:, 0].max()
+ xmin = data[:, 0].min()
+ ymax = data[:, 1].max()
+ ymin = data[:, 1].min()
+
+ # get margin to go beyond the extent to avoid cutting of countour lines
+ x_margin = (xmax - xmin) * margin
+ y_margin = (ymax - ymin) * margin
+
+ if isinstance(resolution, tuple):
+ x_res, y_res = resolution
+ elif isinstance(resolution, (float, int)):
+ x_res = resolution
+ y_res = resolution
+ elif resolution is None:
+ x_res = 100
+ y_res = 100
+ else:
+ raise ValueError("Unsupported option for `resolution`.")
+
+ # create mesh for predicting KDE on with more space around the points
+ x_mesh, y_mesh = np.meshgrid(
+ np.linspace(xmin - x_margin, xmax + x_margin, x_res),
+ np.linspace(ymin - y_margin, ymax + y_margin, y_res),
+ )
+
+ # get the prediction
+ pred = dens_u.pdf(np.vstack([x_mesh.flatten(), y_mesh.flatten()]).T)
+ z = pred.reshape(x_mesh.shape)
+
+ else:
+ kde = FFTKDE(bw=bandwidth, kernel=kernel)
+ grid, points = kde.fit(data).evaluate(resolution)
+ x_mesh, y_mesh = np.unique(grid[:, 0]), np.unique(grid[:, 1])
+ z = points.reshape(resolution, resolution).T
+
+ if fill:
+ ax.contourf(x_mesh, y_mesh, z, levels=levels, **kwargs)
+ else:
+ ax.contour(x_mesh, y_mesh, z, levels=levels, **kwargs)
+
+ return ax
+
+
+"""
+Planar Point Pattern Class
+
+"""
+import numpy as np
+import sys
+from libpysal.cg import KDTree
+from .centrography import hull
+from .window import as_window, poly_from_bbox
+from .util import cached_property
+import pandas as pd
+from matplotlib import pyplot as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Polygon
+
+__author__ = "Serge Rey sjsrey@gmail.com"
+__all__ = ['PointPattern']
+
+if sys.version_info[0] > 2:
+ xrange = range
+
+
+
+[docs]
+class PointPattern(object):
+ """
+ Planar Point Pattern Class 2-D.
+
+ Parameters
+ ----------
+ points: array
+ (n,p), n points with p >= 2 attributes on each
+ point. Two attributes must comprise the spatial
+ coordinate pair. Default is that the first two
+ attributes are the x and y spatial coordinates.
+ window: :class:`.Window`
+ Bounding geometric object for the point pattern.
+ If not specified, window will be set to the minimum
+ bounding rectangle of the point pattern.
+ names: list
+ The names of the attributes.
+ coord_names: list
+ The names of the attributes defining the two spatial
+ coordinates.
+
+ Examples
+ --------
+ >>> from pointpats import PointPattern
+ >>> points = [[66.22, 32.54], [22.52, 22.39], [31.01, 81.21],
+ ... [9.47, 31.02], [30.78, 60.10], [75.21, 58.93],
+ ... [79.26, 7.68], [8.23, 39.93], [98.73, 77.17],
+ ... [89.78, 42.53], [65.19, 92.08], [54.46, 8.48]]
+ >>> pp = PointPattern(points)
+ >>> pp.n
+ 12
+ >>> pp.mean_nnd
+ 21.612139802089246
+ >>> pp.lambda_mbb
+ 0.0015710507711240867
+ >>> pp.lambda_hull
+ 0.0022667153468973137
+ >>> pp.hull_area
+ 5294.00395
+ >>> pp.mbb_area
+ 7638.200000000001
+
+ """
+
+[docs]
+ def __init__(self, points, window=None, names=None, coord_names=None):
+
+ # first two series in df are x, y unless coor_names and names are
+ # specified
+
+ self.df = pd.DataFrame(points)
+ n, p = self.df.shape
+ self._n_marks = p - 2
+ if coord_names is None:
+ if names is not None:
+ coord_names = names[:2]
+ else:
+ coord_names = ['x', 'y']
+ if names is None:
+ col_names = coord_names
+ if p > 2:
+ for m in range(2, p):
+ col_names.append("mark_{}".format(m-2))
+ coord_names = coord_names[:2]
+ else:
+ col_names = names
+
+ self.coord_names = coord_names
+ self._x, self._y = coord_names
+ self.df.columns = col_names
+ self.points = self.df.loc[:, [self._x, self._y]]
+ self._n, self._p = self.points.shape
+ if window is None:
+ self.set_window(as_window(poly_from_bbox(self.mbb)))
+ else:
+ self.set_window(window)
+
+ self._facade()
+
+
+ def __len__(self):
+ """Return the number of points. Use the expression 'len(pp)'.
+
+ Returns
+ -------
+ length : int
+ The number of points in the point pattern.
+
+ Examples
+ --------
+ >>> from pointpats import PointPattern
+ >>> points = [[1, 3], [4, 5], [0,0]]
+ >>> pp = PointPattern(points)
+ >>> len(pp)
+ 3
+
+ """
+ return len(self.df)
+
+ def __contains__(self, n):
+ """Return True if n is a point (a tuple of coordinates), False otherwise.
+ Use the expression 'n in pp'.
+
+ Examples
+ --------
+ >>> from pointpats import PointPattern
+ >>> points = [[1, 3], [4, 5], [0,0]]
+ >>> pp = PointPattern(points)
+ >>> [1, 3] in pp
+ True
+ """
+ name = self.df.columns.values.tolist()
+ return ((self.df[name[0]] == n[0]) & (self.df[name[1]] == n[1])).any()
+
+
+[docs]
+ def set_window(self, window):
+ try:
+ self._window = window
+ except:
+ print("not a valid Window object")
+
+
+
+[docs]
+ def get_window(self):
+ """
+ Bounding geometry for the point pattern
+
+ :class:`.window.Window`
+ """
+ if not hasattr(self, '_window') or self._window is None:
+ # use bbox as window
+ self.set_window(as_window(poly_from_bbox(self.mbb)))
+ return self._window
+
+
+ window = property(get_window, set_window)
+
+
+[docs]
+ def summary(self):
+ '''
+ Description of the point pattern.
+ '''
+
+ print('Point Pattern')
+ print("{} points".format(self.n))
+ print("Bounding rectangle [({},{}), ({},{})]".format(*self.mbb))
+ print("Area of window: {}".format(self.window.area))
+ print("Intensity estimate for window: {}".format(self.lambda_window))
+ print(self.head())
+
+
+
+[docs]
+ def add_marks(self, marks, mark_names=None):
+ if mark_names is None:
+ nm = range(len(marks))
+ mark_names = ["mark_{}".format(self._n_marks+1+j) for j in nm]
+ for name, mark in zip(mark_names, marks):
+ self.df[name] = mark
+ self._n_marks += 1
+
+
+
+[docs]
+ def plot(self, window=False, title="Point Pattern", hull=False,
+ get_ax=False):
+ """
+ Plot function for a point pattern.
+
+ Parameters
+ ----------
+ window : boolean
+ If window is True, plot window of the point
+ pattern. If not, don't plot window.
+ title : string
+ Name of the figure.
+ hull : boolean
+ If hull is True, plot convex hull of the point
+ pattern. If not, don't plot convex hull.
+ get_ax : boolean
+ If get_ax is True, return the current plot ax.
+
+ Returns
+ -------
+ ax : matplotlib.axes._subplots.AxesSubplot
+ Current plot ax. Only return it when get_ax is True.
+
+ """
+ fig, ax = plt.subplots()
+ plt.plot(self.df[self._x], self.df[self._y], '.')
+ # plt.scatter(self.df[self._x], self.df[self._y])
+ plt.title(title)
+ if window:
+ patches = []
+ for part in self.window.parts:
+ p = Polygon(np.asarray(part))
+ patches.append(p)
+ ax.add_collection(
+ PatchCollection(patches, facecolor="none", edgecolor="k", alpha=0.3)
+ )
+ if hull:
+ patches = []
+ p = Polygon(self.hull)
+ patches.append(p)
+ ax.add_collection(
+ PatchCollection(patches, facecolor="none", edgecolor="k", alpha=0.3)
+ )
+ ax.set_aspect("equal")
+
+ # plt.plot(x, y, '.')
+ if get_ax:
+ return ax
+
+
+ def _mbb(self):
+ """
+ Minimum bounding box
+ """
+ mins = self.points.min(axis=0)
+ maxs = self.points.max(axis=0)
+ return np.hstack((mins, maxs))
+
+ mbb = cached_property(_mbb)
+
+ def _mbb_area(self):
+ """
+ Area of minimum bounding box
+ """
+
+ return np.product(self.mbb[[2, 3]]-self.mbb[[0, 1]])
+
+ mbb_area = cached_property(_mbb_area)
+
+ def _n(self):
+ """
+ Number of points
+ """
+ return self.points.shape[0]
+
+ n = cached_property(_n)
+
+ def _rot(self):
+ """
+ Ripley's rule of thumb for distance range in plotting k and related functions
+
+ One-quarter the smallest side of the mbb.
+ """
+ w, s, e, n = self.mbb
+ return 0.25 * min(e-w, n-s)
+
+ rot = cached_property(_rot)
+
+ def _lambda_mbb(self):
+ """
+ Intensity based on minimum bounding box
+ """
+ return self.n * 1. / self.mbb_area
+
+ lambda_mbb = cached_property(_lambda_mbb)
+
+ def _hull(self):
+ """
+ Points defining convex hull in counterclockwise order
+ """
+ return hull(self.points)
+
+ hull = cached_property(_hull)
+
+ def _lambda_window(self):
+ """
+ Intensity estimate based on area of window
+
+ The intensity of a point process at point :math:`s_j` can be defined
+ as:
+
+ .. math::
+
+ \\lambda(s_j) = \\lim \\limits_{|\\mathbf{A}s_j|
+ \\to 0} \\left \\{ \\frac{E(Y(\mathbf{A}s_j)}{|\mathbf{A}s_j|}
+ \\right \\}
+
+ where :math:`\\mathbf{A}s_j` is a small region surrounding location
+ :math:`s_j` with area :math:`|\\mathbf{A}s_j|`, and
+ :math:`E(Y(\\mathbf{A}s_j))` is the expected number of event points in
+ :math:`\\mathbf{A}s_j`.
+
+ The intensity is the mean number of event points per unit of area at
+ point :math:`s_j`.
+
+ """
+ return self.n / self.window.area
+
+ lambda_window = cached_property(_lambda_window)
+
+ def _hull_area(self):
+ """
+ Area of convex hull
+ """
+ h = self.hull
+ if not np.alltrue(h[0] == h[-1]):
+ # not in closed cartographic form
+ h = np.vstack((h, h[0]))
+ s = h[:-1, 0] * h[1:, 1] - h[1:, 0] * h[:-1, 1]
+ return s.sum() / 2.
+
+ hull_area = cached_property(_hull_area)
+
+ def _lambda_hull(self):
+ """
+ Intensity based on convex hull
+ """
+ return self.n * 1. / self.hull_area
+
+ lambda_hull = cached_property(_lambda_hull)
+
+ def _build_tree(self):
+ return KDTree(self.points)
+ tree = cached_property(_build_tree)
+
+
+[docs]
+ def knn(self, k=1):
+ """
+ Find k nearest neighbors for each point in the pattern
+
+ Parameters
+ ----------
+ k: int
+ number of nearest neighbors to find
+
+ Returns
+ -------
+ nn: array (n x k)
+ row i column j contains the id for i's jth nearest neighbor
+
+ nnd: array(n x k)
+ row i column j contains the distance between i and its jth
+ nearest neighbor
+ """
+ if k < 1:
+ raise ValueError('k must be at least 1')
+ nn = self.tree.query(self.tree.data, k=k+1)
+ return nn[1][:, 1:], nn[0][:, 1:]
+
+
+ def _nn_sum(self):
+ """
+ Nearest neighbor distances
+ """
+ ids, nnd = self.knn(1)
+ return nnd
+
+ nnd = cached_property(_nn_sum) # nearest neighbor distances
+
+ def _min_nnd(self):
+ """
+ Min nearest neighbor distance
+ """
+ return self.nnd.min()
+
+ min_nnd = cached_property(_min_nnd)
+
+ def _max_nnd(self):
+ """
+ Max nearest neighbor distance
+ """
+ return self.nnd.max()
+
+ max_nnd = cached_property(_max_nnd)
+
+ def _mean_nnd(self):
+ """
+ Mean nearest neighbor distance
+ """
+ return self.nnd.mean()
+
+ mean_nnd = cached_property(_mean_nnd)
+
+
+[docs]
+ def find_pairs(self, r):
+ """
+ Find all pairs of points in the pattern that are within r units of each
+ other
+
+ Parameters
+ ----------
+ r: float
+ diameter of pair circle
+
+ Returns
+ -------
+ s: set
+ pairs of points within r units of each other
+
+ """
+ return self.tree.query_pairs(r)
+
+
+
+[docs]
+ def knn_other(self, other, k=1):
+ """
+ Find k nearest neighbors in the pattern for each point in other
+
+ Parameters
+ ----------
+ other: PointPattern
+ :py:class:`pointpats.PointPattern`
+ k: int
+ number of nearest neighbors to find
+
+ Returns
+ -------
+ nn: array (n x k)
+ row i column j contains the id for i's jth nearest neighbor
+
+ nnd: array(n x k)
+ row i column j contains the distance between i and its jth
+ nearest neighbor
+ """
+ if k < 1:
+ raise ValueError('k must be at least 1')
+ try:
+ nn = self.tree.query(np.asarray(other.points), k=k)
+ except:
+ nn = self.tree.query(np.asarray(other), k=k)
+ return nn[1], nn[0]
+
+
+
+[docs]
+ def explode(self, mark):
+ """
+ Explode a marked point pattern into a sequence of individual point
+ patterns. If the mark has k unique values, then the sequence will be of
+ length k.
+
+ Parameters
+ ----------
+ mark: string
+ The label of the mark to use for the subsetting
+
+ Returns
+ -------
+ pps: list
+ sequence of :class:`PointPattern` instances
+ """
+
+ uv = np.unique(self.df[mark])
+ pps = [self.df[self.df[mark] == v] for v in uv]
+ names = self.df.columns.values.tolist()
+ cnames = self.coord_names
+ return[PointPattern(pp, names=names, coord_names=cnames) for pp in pps]
+
+
+
+[docs]
+ def unique(self):
+ """ Remove duplicate points in the point pattern.
+
+ Two points in a point pattern are deemed to be identical if their
+ coordinates are the same, and their marks are the same (if any)
+
+ Returns
+ -------
+ pp: list
+ A deduplicated :class:`PointPattern` instance
+
+ Examples
+ --------
+ >>> from pointpats import PointPattern
+ >>> points = [[1.2, 2.1], [1.2, 2.1], [0, 1], [1, 2]]
+ >>> pp = PointPattern(points)
+ >>> pp.unique().df
+ x y
+ 0 1.2 2.1
+ 2 0.0 1.0
+ 3 1.0 2.0
+ """
+ names = self.df.columns.values.tolist()
+ coord_names = self.coord_names
+ window = self.set_window
+ unique_df = self.df.drop_duplicates()
+ return PointPattern(unique_df, names=names, coord_names=coord_names,
+ window=window)
+
+
+
+[docs]
+ def superimpose(self, point_pattern):
+ """Returns a superimposed point pattern.
+
+ Parameters
+ ----------
+ point_pattern:
+ :class:`PointPattern` instance
+
+ Returns
+ -------
+ superimposed :
+ :class:`PointPattern` instance
+
+ Examples
+ --------
+ >>> from pointpats import PointPattern
+ >>> points1 = [[1, 3], [4, 5], [0, 0]]
+ >>> points2 = [[5, 6], [1, 4], [0, 0]]
+ >>> pp1 = PointPattern(points1)
+ >>> pp2 = PointPattern(points2)
+ >>> pp1.superimpose(pp2).points
+ x y
+ 0 1 3
+ 1 4 5
+ 2 0 0
+ 0 5 6
+ 1 1 4
+ """
+ names_pp1 = self.df.columns.values.tolist()
+ cnames_pp1 = self.coord_names
+ names_pp2 = point_pattern.df.columns.values.tolist()
+ cnames_pp2 = point_pattern.coord_names
+ if names_pp1 != names_pp2 or cnames_pp1 != cnames_pp2:
+ raise TypeError('Both point patterns should have similar\
+ attributes and spatial coordinates ')
+ pp = pd.concat((self.df, point_pattern.df))
+ pp = pp.drop_duplicates()
+ return PointPattern(pp, names=names_pp1, coord_names=cnames_pp1)
+
+
+
+[docs]
+ def flip_coordinates(self):
+ """ Flips the coordinates of a point pattern.
+
+ Doesn't change the structure of data frame. This function swaps
+ `_x` and `_y` variables, which are used to represent coordinates.
+ """
+ self._x, self._y = self._y, self._x
+
+
+ # Pandas facade
+ def _facade(self):
+ self.head = self.df.head
+ self.tail = self.df.tail
+
+
+"""
+Simulation of planar point processes
+
+TODO
+
+- inhibition process(es)
+- optimize draws for complex windows
+- documentation
+"""
+
+__author__ = "Serge Rey sjsrey@gmail.com"
+__all__ = ["PointProcess", "PoissonPointProcess", "PoissonClusterPointProcess"]
+
+import numpy as np
+import libpysal as ps
+from numpy.random import poisson
+from .pointpattern import PointPattern as PP
+import warnings
+
+warnings.filterwarnings(
+ "ignore", "Objects based on the `Geometry` class will", FutureWarning
+)
+
+
+def runif_in_circle(n, radius=1.0, center=(0.0, 0.0), burn=2, verbose=False):
+ """
+ Generate n points within a circle of given radius.
+
+ Parameters
+ ----------
+ n : int
+ Number of points.
+ radius : float
+ Radius of the circle.
+ center : tuple
+ Coordinates of the center.
+
+ Returns
+ -------
+ : array
+ (n+1, 2), coordinates of generated points as well as
+ the center.
+
+ """
+
+ good = np.zeros((n, 2), float)
+ c = 0
+ r = radius
+ r2 = r * r
+ it = 0
+ while c < n:
+ x = np.random.uniform(-r, r, (burn * n, 1))
+ y = np.random.uniform(-r, r, (burn * n, 1))
+ ids = np.where(x * x + y * y <= r2)
+ candidates = np.hstack((x, y))[ids[0]]
+ nc = candidates.shape[0]
+ need = n - c
+ if nc > need: # more than we need
+ good[c:] = candidates[:need]
+ else: # use them all and keep going
+ good[c : c + nc] = candidates
+ c += nc
+ it += 1
+ if verbose:
+ print("Iterations: {}".format(it))
+ return good + np.asarray(center)
+
+
+
+[docs]
+class PointProcess(object):
+ """
+ Point Process base class.
+
+ Parameters
+ ----------
+ window : :py:class:`~.window.Window`
+ Bounding geometric object to contain point process
+ realizations.
+ n : int
+ Size of each realization.
+ samples : list
+ Number of realizations.
+ asPP : bool
+ Control the data type of value in the "realizations"
+ dictionary. If True, the data type is point
+ pattern as defined in pointpattern.py; if False,
+ the data type is an two-dimensional array.
+
+ Attributes
+ ----------
+ realizations : dictionary
+ The key is the index of each realization, and the
+ value is simulated event points for each
+ realization. The data type of the value is
+ controlled by the parameter "asPP".
+ parameters : dictionary
+ Dictionary of a dictionary.
+ The key is the index of each realization, and the
+ value is a dictionary with the key 'n' and the
+ value size of each realization.
+
+ """
+
+
+[docs]
+ def __init__(self, window, n, samples, asPP=False, **args):
+ super(PointProcess, self).__init__()
+ self.window = window
+ self.n = n
+ self.samples = samples
+ self.args = args
+ self.realizations = {}
+ self.setup()
+ for sample in range(samples):
+ self.realizations[sample] = self.draw(self.parameters[sample])
+ if asPP:
+ for sample in self.realizations:
+ points = self.realizations[sample]
+ self.realizations[sample] = PP(points, window=self.window)
+ warnings.warn(
+ "These point pattern simulators are deprecated! Please replace them"
+ " with equivalent functions in pointpats.random",
+ DeprecationWarning,
+ stacklevel=2,
+ )
+
+
+
+[docs]
+ def draw(self, parameter):
+ """
+ Generate a series of point coordinates within the given window.
+
+ Parameters
+ ----------
+ parameter : dictionary
+ Key: 'n'.
+ Value: size of the realization.
+
+ Returns
+ -------
+ : array
+ A series of point coordinates.
+
+ """
+ c = 0
+ sample = []
+ n = parameter["n"]
+ while c < n:
+ pnts = self.realize(n)
+ pnts = [ps.cg.shapes.Point((x, y)) for x, y in pnts]
+ pins = self.window.filter_contained(pnts)
+ sample.extend(pins)
+ c = len(sample)
+ return np.array([np.asarray(p) for p in sample[:n]])
+
+
+
+
+
+
+
+
+
+
+
+[docs]
+class PoissonPointProcess(PointProcess):
+ """
+ Poisson point process including :math:`N`-conditioned CSR process and
+ :math:`\lambda`-conditioned CSR process.
+
+ Parameters
+ ----------
+ window : :py:class:`~.window.Window`
+ Bounding geometric object to contain point process
+ realizations.
+ n : int
+ Size of each realization.
+ samples : list
+ Number of realizations.
+ conditioning : bool
+ If True, use the :math:`\lambda`-conditioned CSR process,
+ number of events would vary across realizations;
+ if False, use the :math:`N`-conditioned CSR process.
+ asPP : bool
+ Control the data type of value in the "realizations"
+ dictionary. If True, the data type is point
+ pattern as defined in pointpattern.py; if False,
+ the data type is an two-dimensional array.
+
+ Attributes
+ ----------
+ realizations : dictionary
+ The key is the index of each realization, and the
+ value is simulated event points for each
+ realization. The data type of the value is
+ controlled by the parameter "asPP".
+ parameters : dictionary
+ Dictionary of a dictionary.
+ The key is the index of each realization, and the
+ value is a dictionary with the key 'n' and the
+ value:
+ 1. always equal to the parameter n in the case of
+ N-conditioned process.
+ For example, {0:{'n':100},1:{'n':100},2:{'n':100}}
+ 2. randomly generated from a Possion process in
+ the case of lambda-conditioned process.
+ For example, {0:{'n':97},1:{'n':100},2:{'n':98}}
+
+ Examples
+ --------
+ >>> import libpysal as ps
+ >>> import numpy as np
+ >>> from pointpats import Window
+ >>> from libpysal.cg import shapely_ext
+
+ Open the virginia polygon shapefile
+
+ >>> va = ps.io.open(ps.examples.get_path("virginia.shp"))
+
+ Create the exterior polygons for VA from the union of the county shapes
+
+ >>> polys = [shp for shp in va]
+ >>> state = shapely_ext.cascaded_union(polys)
+
+ Create window from virginia state boundary
+
+ >>> window = Window(state.parts)
+
+ 1. Simulate a :math:`N`-conditioned csr process in the same window (10
+ points, 2 realizations)
+
+ >>> np.random.seed(5)
+ >>> samples1 = PoissonPointProcess(window, 10, 2, conditioning=False, asPP=False)
+ >>> samples1.realizations[0] # the first realized event points
+ array([[-81.80326547, 36.77687577],
+ [-78.5166233 , 37.34055832],
+ [-77.21660795, 37.7491503 ],
+ [-79.30361037, 37.40467853],
+ [-78.61625258, 36.61234487],
+ [-81.43369537, 37.13784646],
+ [-80.91302108, 36.60834063],
+ [-76.90806444, 37.95525903],
+ [-76.33475868, 36.62635347],
+ [-79.71621808, 37.27396618]])
+
+ 2. Simulate a :math:`\lambda`-conditioned csr process in the same window (10
+ points, 2 realizations)
+
+ >>> np.random.seed(5)
+ >>> samples2 = PoissonPointProcess(window, 10, 2, conditioning=True, asPP=True)
+ >>> samples2.realizations[0].n # the size of first realized point pattern
+ 10
+ >>> samples2.realizations[1].n # the size of second realized point pattern
+ 13
+
+ """
+
+
+[docs]
+ def __init__(self, window, n, samples, conditioning=False, asPP=False):
+ self.conditioning = conditioning
+ super(PoissonPointProcess, self).__init__(window, n, samples, asPP)
+
+
+
+[docs]
+ def setup(self):
+ """
+ Generate the number of events for each realization. If
+ "conditioning" is False, all the event numbers are the same;
+ if it is True, the event number is a random variable
+ following a Poisson distribution.
+
+ """
+
+ self.parameters = {}
+ if self.conditioning:
+ lambdas = poisson(self.n, self.samples)
+ for i, l in enumerate(lambdas):
+ self.parameters[i] = {"n": l}
+ else:
+ for i in range(self.samples):
+ self.parameters[i] = {"n": self.n}
+
+
+
+[docs]
+ def realize(self, n):
+ """
+ Generate n points which are randomly and independently
+ distributed in the minimum bounding box of "window".
+
+ Parameters
+ ----------
+ n : int
+ Number of point events.
+
+ Returns
+ -------
+ : array
+ (n,2), n point coordinates.
+
+ """
+
+ l, b, r, t = self.window.bbox
+ xs = np.random.uniform(l, r, (n, 1))
+ ys = np.random.uniform(b, t, (n, 1))
+ return zip(xs, ys)
+
+
+
+
+
+[docs]
+class PoissonClusterPointProcess(PointProcess):
+ """
+ Poisson cluster point process (Neyman Scott).
+ Two stages:
+ 1. parent CSR process: :math:`N`-conditioned or
+ :math:`\lambda`-conditioned. If parent events follow a
+ :math:`\lambda`-conditioned CSR process,
+ the number of parent events varies across realizations.
+ 2. child process: fixed number of points in circle centered
+ on each parent.
+
+ Parameters
+ ----------
+ window : :py:class:`~.window.Window`
+ Bounding geometric object to contain point process
+ realizations.
+ n : int
+ Size of each realization.
+ parents : int
+ Number of parents.
+ radius : float
+ Radius of the circle centered on each parent.
+ samples : list
+ Number of realizations.
+ asPP : bool
+ Control the data type of value in the "realizations"
+ dictionary. If True, the data type is point
+ pattern as defined in pointpattern.py; if False,
+ the data type is an two-dimensional array.
+ conditioning : bool
+ If True, use the :math:`lambda`-conditioned CSR process
+ for parent events, leading to varied number of
+ parent events across realizations;
+ if False, use the :math:`N`-conditioned CSR process.
+
+ Attributes
+ ----------
+ children : int
+ Number of childrens centered on each parent. Can
+ be considered as local intensity.
+ num_parents : dictionary
+ The key is the index of each realization. The
+ value is the number of parent events for each
+ realization.
+ realizations : dictionary
+ The key is the index of each realization, and the
+ value is simulated event points for each
+ realization. The data type of the value is
+ controlled by the parameter "asPP".
+ parameters : dictionary
+ Dictionary of a dictionary.
+ The key is the index of each realization, and the
+ value is a dictionary with the key 'n' and the
+ value always equal to the parameter n in the
+ case of
+ N-conditioned process.
+ For example, {0:{'n':100},1:{'n':100},2:{'n':100}}
+ 2. randomly generated from a Possion process in
+ the case of lambda-conditioned process.
+ For example, {0:{'n':97},1:{'n':100},2:{'n':98}}
+
+ Examples
+ --------
+ >>> import libpysal as ps
+ >>> import numpy as np
+ >>> from pointpats import Window
+ >>> from libpysal.cg import shapely_ext
+
+ Open the virginia polygon shapefile
+
+ >>> va = ps.io.open(ps.examples.get_path("virginia.shp"))
+
+ Create the exterior polygons for VA from the union of the county shapes
+
+ >>> polys = [shp for shp in va]
+ >>> state = shapely_ext.cascaded_union(polys)
+
+ Create window from virginia state boundary
+
+ >>> window = Window(state.parts)
+
+ 1. Simulate a Poisson cluster process of size 200 with 10 parents
+ and 20 children within 0.5 units of each parent
+ (parent events: :math:`N`-conditioned CSR)
+
+ >>> np.random.seed(10)
+ >>> samples1 = PoissonClusterPointProcess(window, 200, 10, 0.5, 1, asPP=True, conditioning=False)
+ >>> samples1.parameters # number of events for the realization
+ {0: {'n': 200}}
+ >>> samples1.num_parents #number of parent events for each realization
+ {0: 10}
+ >>> samples1.children # number of children events centered on each parent event
+ 20
+
+ 2. Simulate a Poisson cluster process of size 200 with 10 parents
+ and 20 children within 0.5 units of each parent
+ (parent events: :math:`\lambda`-conditioned CSR)
+
+ >>> np.random.seed(10)
+ >>> samples2 = PoissonClusterPointProcess(window, 200, 10, 0.5, 1, asPP=True, conditioning=True)
+ >>> samples2.parameters # number of events for the realization might not be equal to 200
+ {0: {'n': 260}}
+ >>> samples2.num_parents #number of parent events for each realization
+ {0: 13}
+ >>> samples2.children # number of children events centered on each parent event
+ 20
+
+ """
+
+
+[docs]
+ def __init__(
+ self,
+ window,
+ n,
+ parents,
+ radius,
+ samples,
+ keep=False,
+ asPP=False,
+ conditioning=False,
+ ):
+ self.conditioning = conditioning
+ self.parents = parents
+ self.children = int(np.ceil(n * 1.0 / parents))
+ self.radius = radius
+ self.keep = keep
+ super(PoissonClusterPointProcess, self).__init__(window, n, samples, asPP)
+
+
+
+[docs]
+ def setup(self):
+ """
+ Generate the number of events for each realization. If
+ "conditioning" is False, all the event numbers are the same;
+ if it is True, the number of parents is a random variable
+ following a Poisson distribution, resulting in varied number
+ of events.
+
+ """
+
+ self.parameters = {}
+ self.num_parents = {}
+ if self.conditioning:
+ lambdas = poisson(self.parents, self.samples)
+ for i, l in enumerate(lambdas):
+ num = l * self.children
+ self.parameters[i] = {"n": num}
+ self.num_parents[i] = l
+ else:
+ for i in range(self.samples):
+ self.parameters[i] = {"n": self.n}
+ self.num_parents[i] = self.parents
+
+
+
+[docs]
+ def realize(self, n):
+ """
+ Generate n points which are distributed in a clustered
+ fashion in the minimum bounding box of "window".
+
+ Parameters
+ ----------
+ n : int
+ Number of point events.
+
+ Returns
+ -------
+ res : array
+ (n,2), n point coordinates.
+
+ """
+ l, b, r, t = self.window.bbox
+ d = self.radius
+ # get parent points
+ pxs = np.random.uniform(l, r, (int(n / self.children), 1))
+ pys = np.random.uniform(b, t, (int(n / self.children), 1))
+ cents = np.hstack((pxs, pys))
+ # generate children points
+ pnts = [runif_in_circle(self.children, d, center) for center in cents]
+ res = np.vstack(np.asarray(pnts))
+ if self.keep:
+ res = np.vstack((np.asarray(cents), res))
+ np.random.shuffle(res) # so we don't truncate in a biased fashion
+ return res
+
+
+
+"""
+Quadrat statistics for planar point patterns
+
+
+TODO
+
+- use patch in matplotlib to plot rectangles and hexagons
+- plot chi2 statistics in each cell
+- delete those cells that do not intersect with the window (study area)
+
+"""
+
+__author__ = "Serge Rey, Wei Kang, Hu Shao"
+__all__ = ["RectangleM", "HexagonM", "QStatistic"]
+
+from .pointpattern import PointPattern
+import numpy as np
+from matplotlib import pyplot as plt
+import math
+import scipy
+
+
+
+[docs]
+class RectangleM:
+ """
+ Rectangle grid structure for quadrat-based method.
+
+ Parameters
+ ----------
+ pp : :class:`.PointPattern`
+ Point Pattern instance.
+ count_column : integer
+ Number of rectangles in the horizontal
+ direction. Use in pair with count_row to
+ fully specify a rectangle. Incompatible with
+ rectangle_width and rectangle_height.
+ count_row : integer
+ Number of rectangles in the vertical
+ direction. Use in pair with count_column to
+ fully specify a rectangle. Incompatible with
+ rectangle_width and rectangle_height.
+ rectangle_width : float
+ Rectangle width. Use in pair with
+ rectangle_height to fully specify a rectangle.
+ Incompatible with count_column & count_row.
+ rectangle_height : float
+ Rectangle height. Use in pair with
+ rectangle_width to fully specify a rectangle.
+ Incompatible with count_column & count_row.
+
+ Attributes
+ ----------
+ pp : :class:`.PointPattern`
+ Point Pattern instance.
+ mbb : array
+ Minimum bounding box for the point pattern.
+ points : array
+ x,y coordinates of the point points.
+ count_column : integer
+ Number of columns.
+ count_row : integer
+ Number of rows.
+ num : integer
+ Number of rectangular quadrats.
+ rectangle_width : float
+ Width of a rectangular quadrat.
+ rectangle_height : float
+ Height of a rectangular quadrat.
+
+ """
+
+
+[docs]
+ def __init__(
+ self, pp, count_column=3, count_row=3, rectangle_width=0, rectangle_height=0
+ ):
+ self.mbb = pp.mbb
+ self.pp = pp
+ self.points = np.asarray(pp.points)
+ x_range = self.mbb[2] - self.mbb[0]
+ y_range = self.mbb[3] - self.mbb[1]
+ if rectangle_width and rectangle_height:
+ self.rectangle_width = rectangle_width
+ self.rectangle_height = rectangle_height
+
+ # calculate column count and row count
+ self.count_column = int(math.ceil(x_range / rectangle_width))
+ self.count_row = int(math.ceil(y_range / rectangle_height))
+ else:
+ self.count_column = count_column
+ self.count_row = count_row
+
+ # calculate the actual width and height of cell
+ self.rectangle_width = x_range / float(count_column)
+ self.rectangle_height = y_range / float(count_row)
+ self.num = self.count_column * self.count_row
+
+
+
+[docs]
+ def point_location_sta(self):
+ """
+ Count the point events in each cell.
+
+ Returns
+ -------
+ dict_id_count : dict
+ keys: rectangle id, values: number of point
+ events in each cell.
+ """
+
+ dict_id_count = {}
+ for i in range(self.count_row):
+ for j in range(self.count_column):
+ dict_id_count[j + i * self.count_column] = 0
+
+ for point in self.points:
+ index_x = (point[0] - self.mbb[0]) // self.rectangle_width
+ index_y = (point[1] - self.mbb[1]) // self.rectangle_height
+ if index_x == self.count_column:
+ index_x -= 1
+ if index_y == self.count_row:
+ index_y -= 1
+ id = index_y * self.count_column + index_x
+ dict_id_count[id] += 1
+ return dict_id_count
+
+
+
+[docs]
+ def plot(self, title="Quadrat Count"):
+ """
+ Plot rectangle tessellation as well as the number of points falling in each rectangle.
+
+ Parameters
+ ----------
+ title: str, optional
+ Title of the plot. Default is "Quadrat Count".
+
+ """
+
+ line_width_cell = 1
+ line_color_cell = "red"
+
+ x_min = self.mbb[0]
+ y_min = self.mbb[1]
+
+ # draw the point pattern along with its window
+ ax = self.pp.plot(window=True, title=title, get_ax=True)
+
+ # draw cells and counts
+ x_start_end = [x_min, x_min + self.count_column * self.rectangle_width]
+ for row in range(self.count_row + 1):
+ y = y_min + row * self.rectangle_height
+ ax.plot(x_start_end, [y, y], lw=line_width_cell, color=line_color_cell)
+ y_start_end = [y_min, y_min + self.count_row * self.rectangle_height]
+ for column in range(self.count_column + 1):
+ x = x_min + column * self.rectangle_width
+ ax.plot([x, x], y_start_end, lw=line_width_cell, color=line_color_cell)
+
+ dict_id_count = self.point_location_sta()
+ for x in range(self.count_column):
+ for y in range(self.count_row):
+ cell_id = x + y * self.count_column
+ count = dict_id_count[cell_id]
+ position_x = x_min + self.rectangle_width * (x + 0.5)
+ position_y = y_min + self.rectangle_height * (y + 0.5)
+ ax.text(position_x, position_y, str(count), ha="center", va="center")
+ return ax
+
+
+
+
+
+[docs]
+class HexagonM:
+ """
+ Hexagon grid structure for quadrat-based method.
+
+ Parameters
+ ----------
+ pp : :class:`.PointPattern`
+ Point Pattern instance.
+ lh : float
+ Hexagon length (hexagon).
+
+ Attributes
+ ----------
+ pp : :class:`.PointPattern`
+ Point Pattern instance.
+ h_length : float
+ Hexagon length (hexagon).
+ mbb : array
+ Minimum bounding box for the point pattern.
+ points : array
+ x,y coordinates of the point points.
+ h_length : float
+ Hexagon length (hexagon).
+ count_row_even : integer
+ Number of even rows.
+ count_row_odd : integer
+ Number of odd rows.
+ count_column : integer
+ Number of columns.
+ num : integer
+ Number of hexagonal quadrats.
+
+ """
+
+
+[docs]
+ def __init__(self, pp, lh):
+
+ self.points = np.asarray(pp.points)
+ self.pp = pp
+ self.h_length = lh
+ self.mbb = pp.mbb
+ range_x = self.mbb[2] - self.mbb[0]
+ range_y = self.mbb[3] - self.mbb[1]
+
+ # calculate column count
+ self.count_column = 1
+ if self.h_length / 2.0 < range_x:
+ temp = math.ceil((range_x - self.h_length / 2) / (1.5 * self.h_length))
+ self.count_column += int(temp)
+
+ # calculate row count for the even columns
+ self.semi_height = self.h_length * math.cos(math.pi / 6)
+ self.count_row_even = 1
+ if self.semi_height < range_y:
+ temp = math.ceil((range_y - self.semi_height) / (self.semi_height * 2))
+ self.count_row_even += int(temp)
+
+ # for the odd columns
+ self.count_row_odd = int(math.ceil(range_y / (self.semi_height * 2)))
+
+ # quadrat number
+ self.num = self.count_row_odd * (
+ (self.count_column // 2) + self.count_column % 2
+ ) + self.count_row_even * (self.count_column // 2)
+
+
+
+[docs]
+ def point_location_sta(self):
+ """
+ Count the point events in each hexagon cell.
+
+ Returns
+ -------
+ dict_id_count : dict
+ keys: rectangle id, values: number of point
+ events in each hexagon cell.
+ """
+ semi_cell_length = self.h_length / 2.0
+ dict_id_count = {}
+
+ # even row may be equal with odd row or 1 more than odd row
+ for i in range(self.count_row_even):
+ for j in range(self.count_column):
+ if (
+ self.count_row_even != self.count_row_odd
+ and i == self.count_row_even - 1
+ ):
+ if j % 2 == 1:
+ continue
+ dict_id_count[j + i * self.count_column] = 0
+
+ x_min = self.mbb[0]
+ y_min = self.mbb[1]
+ x_max = self.mbb[2]
+ y_max = self.mbb[3]
+ points = np.array(self.points)
+ for point in points:
+ # find the possible x index
+ intercept_degree_x = (point[0] - x_min) // semi_cell_length
+
+ # find the possible y index
+ possible_y_index_even = int(
+ (point[1] + self.semi_height - y_min) / (self.semi_height * 2)
+ )
+ possible_y_index_odd = int((point[1] - y_min) / (self.semi_height * 2))
+ if intercept_degree_x % 3 != 1:
+ center_index_x = (intercept_degree_x + 1) // 3
+ center_index_y = possible_y_index_odd
+ if center_index_x % 2 == 0:
+ center_index_y = possible_y_index_even
+ dict_id_count[center_index_x + center_index_y * self.count_column] += 1
+ else: # two columns of cells can be possible
+ center_index_x = intercept_degree_x // 3
+ center_x = center_index_x * semi_cell_length * 3 + x_min
+ center_index_y = possible_y_index_odd
+ center_y = (center_index_y * 2 + 1) * self.semi_height + y_min
+ if center_index_x % 2 == 0:
+ center_index_y = possible_y_index_even
+ center_y = center_index_y * self.semi_height * 2 + y_min
+
+ if point[1] > center_y: # compare the upper bound
+ x0 = center_x + self.h_length
+ y0 = center_y
+ x1 = center_x + semi_cell_length
+ y1 = center_y + self.semi_height
+ indicator = -(
+ point[1]
+ - (
+ (y0 - y1) / (x0 - x1) * point[0]
+ + (x0 * y1 - x1 * y0) / (x0 - x1)
+ )
+ )
+ else: # compare the lower bound
+ x0 = center_x + semi_cell_length
+ y0 = center_y - self.semi_height
+ x1 = center_x + self.h_length
+ y1 = center_y
+ indicator = point[1] - (
+ (y0 - y1) / (x0 - x1) * point[0]
+ + (x0 * y1 - x1 * y0) / (x0 - x1)
+ )
+ if indicator <= 0:
+ # we select right hexagon instead of the left
+ center_index_x += 1
+ center_index_y = possible_y_index_odd
+ if center_index_x % 2 == 0:
+ center_index_y = possible_y_index_even
+ dict_id_count[center_index_x + center_index_y * self.count_column] += 1
+ return dict_id_count
+
+
+
+[docs]
+ def plot(self, title="Quadrat Count"):
+ """
+ Plot hexagon quadrats as well as the number of points falling in each quadrat.
+
+ Parameters
+ ----------
+ title: str, optional
+ Title of the plot. Default is "Quadrat Count".
+
+ """
+ line_width_cell = 1
+ line_color_cell = "red"
+
+ # draw the point pattern along with its window
+ ax = self.pp.plot(window=True, title=title, get_ax=True)
+
+ x_min = self.mbb[0]
+ y_min = self.mbb[1]
+
+ # draw cells and counts
+ dict_id_count = self.point_location_sta()
+ for id in dict_id_count.keys():
+ index_x = id % self.count_column
+ index_y = id // self.count_column
+ center_x = index_x * self.h_length / 2.0 * 3.0 + x_min
+ center_y = index_y * self.semi_height * 2.0 + y_min
+ if index_x % 2 == 1: # for the odd columns
+ center_y = (index_y * 2.0 + 1) * self.semi_height + y_min
+ list_points_cell = []
+ list_points_cell.append([center_x + self.h_length, center_y])
+ list_points_cell.append(
+ [center_x + self.h_length / 2, center_y + self.semi_height]
+ )
+ list_points_cell.append(
+ [center_x - self.h_length / 2, center_y + self.semi_height]
+ )
+ list_points_cell.append([center_x - self.h_length, center_y])
+ list_points_cell.append(
+ [center_x - self.h_length / 2, center_y - self.semi_height]
+ )
+ list_points_cell.append(
+ [center_x + self.h_length / 2, center_y - self.semi_height]
+ )
+ list_points_cell.append([center_x + self.h_length, center_y])
+ ax.plot(
+ np.array(list_points_cell)[:, 0],
+ np.array(list_points_cell)[:, 1],
+ lw=line_width_cell,
+ color=line_color_cell,
+ )
+
+ ax.text(center_x, center_y, str(dict_id_count[id]), ha="center", va="center")
+ return ax
+
+
+
+
+
+[docs]
+class QStatistic:
+ """
+ Quadrat analysis of point pattern.
+
+ Parameters
+ ----------
+ pp : :class:`.PointPattern` or numpy.ndarray
+ Point Pattern instance, or (n_observations, 2) array
+ that can be used to construct a Point Pattern instance.
+ shape : string
+ Grid structure. Either "rectangle" or "hexagon".
+ Default is "rectangle".
+ nx : integer
+ Number of rectangles in the horizontal
+ direction. Only when shape is specified as
+ "rectangle" will nx be considered.
+ ny : integer
+ Number of rectangles in the vertical direction.
+ Only when shape is specified as "rectangle"
+ will ny be considered.
+ rectangle_width : float
+ Rectangle width. Use in pair with
+ rectangle_height to fully specify a rectangle.
+ Incompatible with nx & ny.
+ rectangle_height : float
+ Rectangle height. Use in pair with
+ rectangle_width to fully specify a rectangle.
+ Incompatible with nx & ny.
+ lh : float
+ Hexagon length (hexagon). Only when shape is
+ specified as "hexagon" will lh be considered.
+ Incompatible with nx & ny.
+ realizations : :class:`PointProcess`
+ Point process instance with more than 1 point
+ pattern realizations which would be used for
+ simulation based inference. Default is 0
+ where no simulation based inference is
+ performed.
+
+ Attributes
+ ----------
+ pp : :class:`.PointPattern`
+ Point Pattern instance.
+ mr : :class:`.RectangleM` or :class:`.HexagonM`
+ RectangleM or HexagonM instance.
+ chi2 : float
+ Chi-squared test statistic for the observed
+ point pattern pp.
+ df : integer
+ Degree of freedom.
+ chi2_pvalue : float
+ p-value based on analytical chi-squared
+ distribution.
+ chi2_r_pvalue : float
+ p-value based on simulated sampling
+ distribution. Only available when
+ realizations is correctly specified.
+ chi2_realizations : array
+ Chi-squared test statistics calculated for
+ all the simulated csr point patterns.
+ """
+
+
+[docs]
+ def __init__(self, pp, shape="rectangle", nx=3, ny=3, rectangle_width=0,
+ rectangle_height=0, lh=10, realizations=0):
+ if isinstance(pp, np.ndarray):
+ pp = PointPattern(pp)
+ self.pp = pp
+ if shape == "rectangle":
+ self.mr = RectangleM(pp, count_column=nx, count_row=ny,
+ rectangle_width=rectangle_width,
+ rectangle_height = rectangle_height)
+ elif shape == "hexagon":
+ self.mr = HexagonM(pp, lh)
+ else:
+ raise ValueError(
+ f'shape type {shape} not understood. Must be either "rectangle" or'
+ ' "hexagon"'
+ )
+
+ # calculate chi2 test statisitc for the observed point pattern
+ dict_id_count = self.mr.point_location_sta()
+ self.chi2, self.chi2_pvalue = scipy.stats.chisquare(
+ list(dict_id_count.values())
+ )
+
+ self.df = self.mr.num - 1
+
+ # when realizations is specified, perform simulation based
+ # inference.
+ if realizations:
+ reals = realizations.realizations
+ sim_n = realizations.samples
+ chi2_realizations = [] # store test statistics for all the
+ # similations
+
+ for i in range(sim_n):
+ if shape == "rectangle":
+ mr_temp = RectangleM(reals[i], count_column=nx, count_row=ny)
+ elif shape == "hexagon":
+ mr_temp = HexagonM(reals[i], lh)
+ id_count_temp = mr_temp.point_location_sta().values()
+
+ # calculate test statistics for simulated point patterns
+ chi2_sim, p = scipy.stats.chisquare(list(id_count_temp))
+ chi2_realizations.append(chi2_sim)
+ self.chi2_realizations = np.array(chi2_realizations)
+
+ # calculate pseudo pvalue
+ above_chi2 = self.chi2_realizations >= self.chi2
+ larger_chi2 = sum(above_chi2)
+ self.chi2_r_pvalue = (larger_chi2 + 1.0) / (sim_n + 1.0)
+
+
+
+[docs]
+ def plot(self, title="Quadrat Count"):
+ """
+ Plot quadrats as well as the number of points falling in each quadrat.
+
+ Parameters
+ ----------
+ title: str, optional
+ Title of the plot. Default is "Quadrat Count".
+
+ """
+
+ return self.mr.plot(title=title)
+
+
+
+import numpy
+from .geometry import (
+ spatial,
+ area as _area,
+ centroid as _centroid,
+ contains as _contains,
+ bbox as _bbox,
+ prepare_hull as _prepare_hull,
+ HULL_TYPES,
+)
+
+# ------------------------------------------------------------ #
+# Utilities #
+# ------------------------------------------------------------ #
+
+
+def parse_size_and_intensity(hull, intensity=None, size=None):
+ """
+ Given a hull, an intensity, and a size int/tuple, correctly
+ compute the resulting missing quantities. Defaults to 100 points in one
+ replication, meaning the intensity will be computed on the fly
+ if nothing is provided.
+
+ Parameters
+ ----------
+ hull : A geometry-like object
+ This encodes the "space" in which to simulate the normal pattern. All points will
+ lie within this hull. Supported values are:
+ - a bounding box encoded in a numpy array as numpy.array([xmin, ymin, xmax, ymax])
+ - an (N,2) array of points for which the bounding box will be computed & used
+ - a shapely polygon/multipolygon
+ - a scipy convexh hull
+ intensity : float
+ the number of observations per unit area in the hull to use. If provided,
+ then the number of observations is determined using the intensity * area(hull) and
+ the size is assumed to represent n_replications (if provided).
+ size : tuple or int
+ a tuple of (n_observations, n_replications), where the first number is the number
+ of points to simulate in each replication and the second number is the number of
+ total replications. So, (10, 4) indicates 10 points, 4 times.
+ If an integer is provided and intensity is None, n_replications is assumed to be 1.
+ If size is an integer and intensity is also provided, then size indicates n_replications,
+ and the number of observations is computed on the fly using intensity and area.
+ """
+ if size is None:
+ if intensity is not None:
+ # if intensity is provided, assume
+ # n_observations
+ n_observations = int(_area(hull) * intensity)
+ else:
+ # default to 100 points
+ n_observations = 100
+ intensity = n_observations / _area(hull)
+ n_simulations = 1
+ size = (n_observations, n_simulations)
+ elif isinstance(size, tuple):
+ if len(size) == 2 and intensity is None:
+ n_observations, n_simulations = size
+ intensity = n_observations / _area(hull)
+ elif len(size) == 2 and intensity is not None:
+ raise ValueError(
+ "Either intensity or size as (n observations, n simulations)"
+ " can be provided. Providing both creates statistical conflicts."
+ " between the requested intensity and implied intensity by"
+ " the number of observations and the area of the hull. If"
+ " you want to specify the intensity, use the intensity argument"
+ " and set size equal to the number of simulations."
+ )
+ else:
+ raise ValueError(
+ f"Intensity and size not understood. Provide size as a tuple"
+ f" containing (number of observations, number of simulations)"
+ f" with no specified intensity, or an intensity and size equal"
+ f" to the number of simulations."
+ f" Recieved: `intensity={intensity}, size={size}`"
+ )
+ elif isinstance(size, int):
+ # assume int size with specified intensity means n_simulations at x intensity
+ if intensity is not None:
+ n_observations = int(intensity * _area(hull))
+ n_simulations = size
+ else: # assume we have one replication at the specified number of points
+ n_simulations = 1
+ n_observations = size
+ intensity = n_observations / _area(hull)
+ else:
+ raise ValueError(
+ f"Intensity and size not understood. Provide size as a tuple"
+ f" containing (number of observations, number of simulations)"
+ f" with no specified intensity, or an intensity and size equal"
+ f" to the number of simulations."
+ f" Recieved: `intensity={intensity}, size={size}`"
+ )
+ return (n_observations, n_simulations, intensity)
+
+
+# ------------------------------------------------------------ #
+# Distributions #
+# ------------------------------------------------------------ #
+
+
+
+[docs]
+def poisson(hull, intensity=None, size=None):
+ """
+ Simulate a poisson random point process with a specified intensity.
+
+ Parameters
+ ----------
+ hull : A geometry-like object
+ This encodes the "space" in which to simulate the normal pattern. All points will
+ lie within this hull. Supported values are:
+ - a bounding box encoded in a numpy array as numpy.array([xmin, ymin, xmax, ymax])
+ - an (N,2) array of points for which the bounding box will be computed & used
+ - a shapely polygon/multipolygon
+ - a scipy convexh hull
+ intensity : float
+ the number of observations per unit area in the hull to use. If provided, then
+ size must be an integer describing the number of replications to use.
+ size : tuple or int
+ a tuple of (n_observations, n_replications), where the first number is the number
+ of points to simulate in each replication and the second number is the number of
+ total replications. So, (10, 4) indicates 10 points, 4 times.
+ If an integer is provided and intensity is None, n_replications is assumed to be 1.
+ If size is an integer and intensity is also provided, then size indicates n_replications,
+ and the number of observations is computed from the intensity.
+
+ Returns
+ --------
+ : numpy.ndarray
+ either an (n_replications, n_observations, 2) or (n_observations,2) array containing
+ the simulated realizations.
+ """
+ if isinstance(hull, numpy.ndarray):
+ if hull.shape == (4,):
+ hull = hull
+ else:
+ hull = _prepare_hull(hull)
+ n_observations, n_simulations, intensity = parse_size_and_intensity(
+ hull, intensity=intensity, size=size
+ )
+
+ result = numpy.empty((n_simulations, n_observations, 2))
+
+ bbox = _bbox(hull)
+
+ for i_replication in range(n_simulations):
+ generating = True
+ i_observation = 0
+ while i_observation < n_observations:
+ x, y = (
+ numpy.random.uniform(bbox[0], bbox[2]),
+ numpy.random.uniform(bbox[1], bbox[3]),
+ )
+ if _contains(hull, x, y):
+ result[i_replication, i_observation] = (x, y)
+ i_observation += 1
+ return result.squeeze()
+
+
+
+
+[docs]
+def normal(hull, center=None, cov=None, size=None):
+ """
+ Simulate a multivariate random normal point cluster
+
+ Parameters
+ ----------
+ hull : A geometry-like object
+ This encodes the "space" in which to simulate the normal pattern. All points will
+ lie within this hull. Supported values are:
+ - a bounding box encoded in a numpy array as numpy.array([xmin, ymin, xmax, ymax])
+ - an (N,2) array of points for which the bounding box will be computed & used
+ - a shapely polygon/multipolygon
+ - a scipy convexh hull
+ center : iterable of shape (2, )
+ A point where the simulations will be centered.
+ cov : float or a numpy array of shape (2,2)
+ either the standard deviation of an independent and identically distributed
+ normal distribution, or a 2 by 2 covariance matrix expressing the covariance
+ of the x and y for the distribution. Default is half of the width or height
+ of the hull's bounding box, whichever is larger.
+ size : tuple or int
+ a tuple of (n_observations, n_replications), where the first number is the number
+ of points to simulate in each replication and the second number is the number of
+ total replications. So, (10, 4) indicates 10 points, 4 times.
+ If an integer is provided, n_replications is assumed to be 1.
+
+ Returns
+ --------
+ : numpy.ndarray
+ either an (n_replications, n_observations, 2) or (n_observations,2) array containing
+ the simulated realizations.
+ """
+ if isinstance(hull, numpy.ndarray):
+ if hull.shape == (4,):
+ hull = hull
+ else:
+ hull = _prepare_hull(hull)
+ if center is None:
+ center = _centroid(hull)
+ n_observations, n_simulations, intensity = parse_size_and_intensity(
+ hull, intensity=None, size=size
+ )
+ if cov is None:
+ bbox = _bbox(hull)
+ width = bbox[2] - bbox[0]
+ height = bbox[3] - bbox[1]
+ cov = numpy.maximum(width / 2, height / 2) ** 2
+
+ if isinstance(cov, (int, float)):
+ sd = cov
+ cov = numpy.eye(2) * sd
+ elif isinstance(cov, numpy.ndarray):
+ if cov.ndim == 2:
+ assert cov.shape == (2, 2), "Bivariate covariance matrices must be 2 by 2"
+ elif cov.ndim == 3:
+ assert cov.shape[1:] == (
+ 2,
+ 2,
+ ), "3-dimensional covariance matrices must have shape (n_simulations, 2,2)"
+ assert (
+ cov.shape[0] == n_simulations
+ ), "3-dimensional covariance matrices must have shape (n_simulations, 2,2)"
+ else:
+ raise ValueError(
+ "`cov` argument must be a float (signifying a standard deviation)"
+ " or a 2 by 2 array expressing the covariance matrix of the "
+ " bivariate normal distribution."
+ )
+
+ result = numpy.empty((n_simulations, n_observations, 2))
+
+ bbox = _bbox(hull)
+
+ for i_replication in range(n_simulations):
+ generating = True
+ i_observation = 0
+ replication_cov = cov[i] if cov.ndim == 3 else cov
+ replication_sd = numpy.diagonal(replication_cov) ** 0.5
+ replication_cor = (1 / replication_sd) * replication_cov * (1 / replication_sd)
+
+ while i_observation < n_observations:
+ candidate = numpy.random.multivariate_normal((0, 0), replication_cor)
+ x, y = center + candidate * replication_sd
+ if _contains(hull, x, y):
+ result[i_replication, i_observation] = (x, y)
+ i_observation += 1
+ return result.squeeze()
+
+
+
+
+[docs]
+def cluster_poisson(
+ hull, intensity=None, size=None, n_seeds=2, cluster_radius=None,
+):
+ """
+ Simulate a cluster poisson random point process with a specified intensity & number of seeds.
+ A cluster poisson process is a poisson process where the center of each "cluster" is
+ itself distributed according to a spatial poisson process.
+
+ Parameters
+ ----------
+ hull : A geometry-like object
+ This encodes the "space" in which to simulate the normal pattern. All points will
+ lie within this hull. Supported values are:
+ - a bounding box encoded in a numpy array as numpy.array([xmin, ymin, xmax, ymax])
+ - an (N,2) array of points for which the bounding box will be computed & used
+ - a shapely polygon/multipolygon
+ - a scipy convexh hull
+ intensity : float
+ the number of observations per unit area in the hull to use. If provided, then
+ size must be an integer describing the number of replications to use.
+ size : tuple or int
+ a tuple of (n_observations, n_replications), where the first number is the number
+ of points to simulate in each replication and the second number is the number of
+ total replications. So, (10, 4) indicates 10 points, 4 times.
+ If an integer is provided and intensity is None, n_replications is assumed to be 1.
+ If size is an integer and intensity is also provided, then size indicates n_replications,
+ and the number of observations is computed from the intensity.
+ n_seeds : int
+ the number of sub-clusters to use.
+ cluster_radius : float or iterable
+ the radius of each cluster. If a float, the same radius is used for all clusters.
+ If an array, then there must be the same number of radii as clusters.
+ If None, 50% of the minimum inter-point distance is used, which may fluctuate across
+ replications.
+
+ Returns
+ --------
+ : numpy.ndarray
+ either an (n_replications, n_observations, 2) or (n_observations,2) array containing
+ the simulated realizations.
+ """
+ if isinstance(hull, numpy.ndarray):
+ if hull.shape == (4,):
+ hull = hull
+ else:
+ hull = _prepare_hull(hull)
+
+ if isinstance(cluster_radius, numpy.ndarray):
+ cluster_radii = cluster_radius.flatten()
+ assert len(cluster_radii) == n_seeds, (
+ f"number of radii provided ({len(cluster_radii)})"
+ f"does not match number of clusters requested"
+ f" ({n_seeds})."
+ )
+ elif isinstance(cluster_radius, (int, float)):
+ cluster_radii = [cluster_radius] * n_seeds
+
+ n_observations, n_simulations, intensity = parse_size_and_intensity(
+ hull, intensity=intensity, size=size
+ )
+
+ result = numpy.empty((n_simulations, n_observations, 2))
+ hull_area = _area(hull)
+ for i_replication in range(n_simulations):
+ seeds = poisson(hull, size=(n_seeds, n_simulations))
+ if cluster_radius is None:
+ # default cluster radius is one half the minimum distance between seeds
+ cluster_radii = [spatial.distance.pdist(seeds).min() * 0.5] * n_seeds
+ clusters = numpy.array_split(result[i_replication], n_seeds)
+ for i_cluster, radius in enumerate(cluster_radii):
+ seed = seeds[i_cluster]
+ cluster_points = clusters[i_cluster]
+ n_in_cluster = len(cluster_points)
+ if n_in_cluster == 1:
+ clusters[i_cluster] = seed
+ continue
+ if n_in_cluster < 1:
+ raise Exception(
+ "There are too many clusters requested for the "
+ " inputted number of samples. Reduce `n_seeds` or"
+ " increase the number of sampled points."
+ )
+ candidates = _uniform_circle(
+ n_in_cluster - 1, radius=radius, center=seed, hull=hull
+ )
+ clusters[i_cluster] = numpy.row_stack((seed, candidates))
+ result[i_replication] = numpy.row_stack(clusters)
+ return result.squeeze()
+
+
+
+
+[docs]
+def cluster_normal(hull, cov=None, size=None, n_seeds=2):
+ """
+ Simulate a cluster poisson random point process with a specified intensity & number of seeds.
+ A cluster poisson process is a poisson process where the center of each "cluster" is
+ itself distributed according to a spatial poisson process.
+
+ Parameters
+ ----------
+ hull : A geometry-like object
+ This encodes the "space" in which to simulate the normal pattern. All points will
+ lie within this hull. Supported values are:
+ - a bounding box encoded in a numpy array as numpy.array([xmin, ymin, xmax, ymax])
+ - an (N,2) array of points for which the bounding box will be computed & used
+ - a shapely polygon/multipolygon
+ - a scipy convexh hull
+ cov : float, int, or numpy.ndarray of shape (2,2)
+ The covariance structure for clusters. By default, this is the squared
+ average distance between cluster seeds.
+ size : tuple or int
+ a tuple of (n_observations, n_replications), where the first number is the number
+ of points to simulate in each replication and the second number is the number of
+ total replications. So, (10, 4) indicates 10 points, 4 times.
+ If an integer is provided and intensity is None, n_replications is assumed to be 1.
+ If size is an integer and intensity is also provided, then size indicates n_replications,
+ and the number of observations is computed from the intensity.
+ n_seeds : int
+ the number of sub-clusters to use.
+
+ Returns
+ --------
+ : numpy.ndarray
+ either an (n_replications, n_observations, 2) or (n_observations,2) array containing
+ the simulated realizations.
+ """
+ if isinstance(hull, numpy.ndarray):
+ if hull.shape == (4,):
+ hull = hull
+ else:
+ hull = _prepare_hull(hull)
+ n_observations, n_simulations, intensity = parse_size_and_intensity(
+ hull, intensity=None, size=size
+ )
+ result = numpy.empty((n_simulations, n_observations, 2))
+ for i_replication in range(n_simulations):
+ seeds = poisson(hull, size=(n_seeds, n_simulations))
+ if cov is None:
+ cov = spatial.distance.pdist(seeds).mean() ** 2
+ clusters = numpy.array_split(result[i_replication], n_seeds)
+ for i_cluster, seed in enumerate(seeds):
+ cluster_points = clusters[i_cluster]
+ n_in_cluster = len(cluster_points)
+ if n_in_cluster == 1:
+ clusters[i_cluster] = seed
+ continue
+ if n_in_cluster < 1:
+ raise Exception(
+ "There are too many clusters requested for the "
+ " inputted number of samples. Reduce `n_seeds` or"
+ " increase the number of sampled points."
+ )
+ candidates = normal(hull, center=seed, cov=cov, size=n_in_cluster - 1)
+ clusters[i_cluster] = numpy.row_stack((seed, candidates))
+ result[i_replication] = numpy.row_stack(clusters)
+ return result.squeeze()
+
+
+
+def _uniform_circle(n, radius=1.0, center=(0.0, 0.0), burn=2, verbose=False, hull=None):
+ """
+ Generate n points within a circle of given radius.
+
+ Parameters
+ ----------
+ n : int
+ Number of points.
+ radius : float
+ Radius of the circle.
+ center : tuple
+ Coordinates of the center.
+ burn : int
+ number of coordinates to simulate at a time. This is the "chunk"
+ size sent to numpy.random.uniform each iteration of the rejection sampler
+
+ Returns
+ -------
+ : array
+ (n, 2), coordinates of generated points
+
+ """
+
+ good = numpy.zeros((n, 2), float)
+ c = 0
+ center_x, center_y = center
+ r = radius
+ r2 = r * r
+ it = 0
+ while c < n:
+ x = numpy.random.uniform(-r, r, (burn * n, 1))
+ y = numpy.random.uniform(-r, r, (burn * n, 1))
+ if hull is None:
+ in_hull = True
+ else:
+ in_hull = numpy.asarray(
+ [
+ _contains(hull, xi + center_x, yi + center_y)
+ for xi, yi in numpy.column_stack((x, y))
+ ]
+ ).reshape(-1, 1)
+ ids, *_ = numpy.where(((x * x + y * y) <= r2) & in_hull)
+ candidates = numpy.hstack((x, y))[ids]
+ nc = candidates.shape[0]
+ need = n - c
+ if nc > need: # more than we need
+ good[c:] = candidates[:need]
+ else: # use them all and keep going
+ good[c : c + nc] = candidates
+ c += nc
+ it += 1
+ if verbose:
+ print("Iterations: {}".format(it))
+ return good + numpy.asarray(center)
+
+"""
+Methods for identifying space-time interaction in spatio-temporal event
+data.
+"""
+
+__author__ = (
+ "Eli Knaap <eknaap@sdsu.edu>",
+ "Nicholas Malizia <nmalizia@asu.edu>",
+ "Sergio J. Rey <srey@sdsu.edu>",
+ "Philip Stephens <philip.stephens@asu.edu",
+)
+
+__all__ = [
+ "SpaceTimeEvents",
+ "knox",
+ "mantel",
+ "jacquez",
+ "modified_knox",
+ "Knox",
+ "KnoxLocal",
+]
+
+import os
+from datetime import date
+from warnings import warn
+
+import geopandas as gpd
+import libpysal as lps
+import numpy as np
+import pandas
+import pandas as pd
+import scipy.stats as stats
+from libpysal import cg
+from libpysal.graph import Graph
+from pandas.api.types import is_numeric_dtype
+from scipy.spatial import KDTree
+from scipy.stats import hypergeom, poisson
+from shapely.geometry import LineString
+
+
+
+[docs]
+class SpaceTimeEvents:
+ """
+ Method for reformatting event data stored in a shapefile for use in
+ calculating metrics of spatio-temporal interaction.
+
+ Parameters
+ ----------
+ path : string
+ the path to the appropriate shapefile, including the
+ file name and extension
+ time : string
+ column header in the DBF file indicating the column
+ containing the time stamp.
+ infer_timestamp : bool, optional
+ if the column containing the timestamp is formatted as
+ calendar dates, try to coerce them into Python datetime
+ objects (the default is False).
+
+ Attributes
+ ----------
+ n : int
+ number of events.
+ x : array
+ (n, 1), array of the x coordinates for the events.
+ y : array
+ (n, 1), array of the y coordinates for the events.
+ t : array
+ (n, 1), array of the temporal coordinates for the events.
+ space : array
+ (n, 2), array of the spatial coordinates (x,y) for the
+ events.
+ time : array
+ (n, 2), array of the temporal coordinates (t,1) for the
+ events, the second column is a vector of ones.
+
+ Examples
+ --------
+ Read in the example shapefile data, ensuring to omit the file
+ extension. In order to successfully create the event data the .dbf file
+ associated with the shapefile should have a column of values that are a
+ timestamp for the events. This timestamp may be a numerical value
+ or a date. Date inference was added in version 1.6.
+
+ >>> import libpysal as lps
+ >>> path = lps.examples.get_path("burkitt.shp")
+ >>> from pointpats import SpaceTimeEvents
+
+ Create an instance of SpaceTimeEvents from a shapefile, where the
+ temporal information is stored in a column named "T".
+
+ >>> events = SpaceTimeEvents(path,'T')
+
+ See how many events are in the instance.
+
+ >>> events.n
+ 188
+
+ Check the spatial coordinates of the first event.
+
+ >>> events.space[0]
+ array([300., 302.])
+
+ Check the time of the first event.
+
+ >>> events.t[0]
+ array([413.])
+
+ Calculate the time difference between the first two events.
+
+ >>> events.t[1] - events.t[0]
+ array([59.])
+
+ New, in 1.6, date support:
+ Now, create an instance of SpaceTimeEvents from a shapefile, where the
+ temporal information is stored in a column named "DATE".
+
+ >>> events = SpaceTimeEvents(path,'DATE')
+
+ See how many events are in the instance.
+
+ >>> events.n
+ 188
+
+ Check the spatial coordinates of the first event.
+
+ >>> events.space[0]
+ array([300., 302.])
+
+ Check the time of the first event. Note that this value is equivalent to
+ 413 days after January 1, 1900.
+
+ >>> events.t[0][0]
+ datetime.date(1901, 2, 16)
+
+ Calculate the time difference between the first two events.
+
+ >>> (events.t[1][0] - events.t[0][0]).days
+ 59
+ """
+
+
+[docs]
+ def __init__(self, path, time_col, infer_timestamp=False):
+ shp = lps.io.open(path)
+ head, tail = os.path.split(path)
+ dbf_tail = tail.split(".")[0] + ".dbf"
+ dbf = lps.io.open(lps.examples.get_path(dbf_tail))
+
+ # extract the spatial coordinates from the shapefile
+ x = [coords[0] for coords in shp]
+ y = [coords[1] for coords in shp]
+
+ self.n = n = len(shp)
+ x = np.array(x)
+ y = np.array(y)
+ self.x = np.reshape(x, (n, 1))
+ self.y = np.reshape(y, (n, 1))
+ self.space = np.hstack((self.x, self.y))
+
+ # extract the temporal information from the database
+ if infer_timestamp:
+ col = dbf.by_col(time_col)
+ if isinstance(col[0], date):
+ day1 = min(col)
+ col = [(d - day1).days for d in col]
+ t = np.array(col)
+ else:
+ print(
+ "Unable to parse your time column as Python datetime \
+ objects, proceeding as integers."
+ )
+ t = np.array(col)
+ else:
+ t = np.array(dbf.by_col(time_col))
+ line = np.ones((n, 1))
+ self.t = np.reshape(t, (n, 1))
+ self.time = np.hstack((self.t, line))
+
+ # close open objects
+ dbf.close()
+ shp.close()
+
+
+
+
+def knox(s_coords, t_coords, delta, tau, permutations=99, debug=False):
+ """
+ Knox test for spatio-temporal interaction. :cite:`Knox:1964`
+
+ Parameters
+ ----------
+ s_coords : array
+ (n, 2), spatial coordinates.
+ t_coords : array
+ (n, 1), temporal coordinates.
+ delta : float
+ threshold for proximity in space.
+ tau : float
+ threshold for proximity in time.
+ permutations : int, optional
+ the number of permutations used to establish pseudo-
+ significance (the default is 99).
+ debug : bool, optional
+ if true, debugging information is printed (the default is
+ False).
+
+ Returns
+ -------
+ knox_result : dictionary
+ contains the statistic (stat) for the test and the
+ associated p-value (pvalue).
+ stat : float
+ value of the knox test for the dataset.
+ pvalue : float
+ pseudo p-value associated with the statistic.
+ counts : int
+ count of space time neighbors.
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import libpysal as lps
+ >>> from pointpats import SpaceTimeEvents, knox
+
+ Read in the example data and create an instance of SpaceTimeEvents.
+
+ >>> path = lps.examples.get_path("burkitt.shp")
+ >>> events = SpaceTimeEvents(path,'T')
+
+ Set the random seed generator. This is used by the permutation based
+ inference to replicate the pseudo-significance of our example results -
+ the end-user will normally omit this step.
+
+ >>> np.random.seed(100)
+
+ Run the Knox test with distance and time thresholds of 20 and 5,
+ respectively. This counts the events that are closer than 20 units in
+ space, and 5 units in time.
+
+ >>> result = knox(events.space, events.t, delta=20, tau=5, permutations=99)
+
+ Next, we examine the results. First, we call the statistic from the
+ results dictionary. This reports that there are 13 events close
+ in both space and time, according to our threshold definitions.
+
+ >>> result['stat'] == 13
+ True
+
+ Next, we look at the pseudo-significance of this value, calculated by
+ permuting the timestamps and rerunning the statistics. In this case,
+ the results indicate there is likely no space-time interaction between
+ the events.
+
+ >>> print("%2.2f"%result['pvalue'])
+ 0.17
+ """
+ warn("This function is deprecated. Use Knox", DeprecationWarning, stacklevel=2)
+
+ # Do a kdtree on space first as the number of ties (identical points) is
+ # likely to be lower for space than time.
+
+ kd_s = cg.KDTree(s_coords)
+ neigh_s = kd_s.query_pairs(delta)
+ tau2 = tau * tau
+ ids = np.array(list(neigh_s))
+
+ # For the neighboring pairs in space, determine which are also time
+ # neighbors
+
+ d_t = (t_coords[ids[:, 0]] - t_coords[ids[:, 1]]) ** 2
+ n_st = sum(d_t <= tau2)
+
+ knox_result = {"stat": n_st[0]}
+
+ if permutations:
+ joint = np.zeros((permutations, 1), int)
+ for p in range(permutations):
+ np.random.shuffle(t_coords)
+ d_t = (t_coords[ids[:, 0]] - t_coords[ids[:, 1]]) ** 2
+ joint[p] = np.sum(d_t <= tau2)
+
+ larger = sum(joint >= n_st[0])
+ if (permutations - larger) < larger:
+ larger = permutations - larger
+ p_sim = (larger + 1.0) / (permutations + 1.0)
+ knox_result["pvalue"] = p_sim
+ return knox_result
+
+
+
+[docs]
+def mantel(
+ s_coords, t_coords, permutations=99, scon=1.0, spow=-1.0, tcon=1.0, tpow=-1.0
+):
+ """
+ Standardized Mantel test for spatio-temporal interaction. :cite:`Mantel:1967`
+
+ Parameters
+ ----------
+ s_coords : array
+ (n, 2), spatial coordinates.
+ t_coords : array
+ (n, 1), temporal coordinates.
+ permutations : int, optional
+ the number of permutations used to establish pseudo-
+ significance (the default is 99).
+ scon : float, optional
+ constant added to spatial distances (the default is 1.0).
+ spow : float, optional
+ value for power transformation for spatial distances
+ (the default is -1.0).
+ tcon : float, optional
+ constant added to temporal distances (the default is 1.0).
+ tpow : float, optional
+ value for power transformation for temporal distances
+ (the default is -1.0).
+
+ Returns
+ -------
+ mantel_result : dictionary
+ contains the statistic (stat) for the test and the
+ associated p-value (pvalue).
+ stat : float
+ value of the knox test for the dataset.
+ pvalue : float
+ pseudo p-value associated with the statistic.
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import libpysal as lps
+ >>> from pointpats import SpaceTimeEvents, mantel
+
+ Read in the example data and create an instance of SpaceTimeEvents.
+
+ >>> path = lps.examples.get_path("burkitt.shp")
+ >>> events = SpaceTimeEvents(path,'T')
+
+ Set the random seed generator. This is used by the permutation based
+ inference to replicate the pseudo-significance of our example results -
+ the end-user will normally omit this step.
+
+ >>> np.random.seed(100)
+
+ The standardized Mantel test is a measure of matrix correlation between
+ the spatial and temporal distance matrices of the event dataset. The
+ following example runs the standardized Mantel test without a constant
+ or transformation; however, as recommended by :cite:`Mantel:1967`, these
+ should be added by the user. This can be done by adjusting the constant
+ and power parameters.
+
+ >>> result = mantel(events.space, events.t, 99, scon=1.0, spow=-1.0, tcon=1.0, tpow=-1.0)
+
+ Next, we examine the result of the test.
+
+ >>> print("%6.6f"%result['stat'])
+ 0.048368
+
+ Finally, we look at the pseudo-significance of this value, calculated by
+ permuting the timestamps and rerunning the statistic for each of the 99
+ permutations. According to these parameters, the results indicate
+ space-time interaction between the events.
+
+ >>> print("%2.2f"%result['pvalue'])
+ 0.01
+
+ """
+
+ t = t_coords
+ s = s_coords
+ n = len(t)
+
+ # calculate the spatial and temporal distance matrices for the events
+ distmat = cg.distance_matrix(s)
+ timemat = cg.distance_matrix(t)
+
+ # calculate the transformed standardized statistic
+ timevec = (timemat[np.tril_indices(timemat.shape[0], k=-1)] + tcon) ** tpow
+ distvec = (distmat[np.tril_indices(distmat.shape[0], k=-1)] + scon) ** spow
+ stat = stats.pearsonr(timevec, distvec)[0].sum()
+
+ # return the results (if no inference)
+ if not permutations:
+ return stat
+
+ # loop for generating a random distribution to assess significance
+ dist = []
+ for _i in range(permutations):
+ trand = _shuffle_matrix(timemat, np.arange(n))
+ timevec = (trand[np.tril_indices(trand.shape[0], k=-1)] + tcon) ** tpow
+ m = stats.pearsonr(timevec, distvec)[0].sum()
+ dist.append(m)
+
+ # establish the pseudo significance of the observed statistic
+ distribution = np.array(dist)
+ greater = np.ma.masked_greater_equal(distribution, stat)
+ count = np.ma.count_masked(greater)
+ pvalue = (count + 1.0) / (permutations + 1.0)
+
+ # report the results
+ mantel_result = {"stat": stat, "pvalue": pvalue}
+ return mantel_result
+
+
+
+
+[docs]
+def jacquez(s_coords, t_coords, k, permutations=99):
+ """
+ Jacquez k nearest neighbors test for spatio-temporal interaction.
+ :cite:`Jacquez:1996`
+
+ Parameters
+ ----------
+ s_coords : array
+ (n, 2), spatial coordinates.
+ t_coords : array
+ (n, 1), temporal coordinates.
+ k : int
+ the number of nearest neighbors to be searched.
+ permutations : int, optional
+ the number of permutations used to establish pseudo-
+ significance (the default is 99).
+
+ Returns
+ -------
+ jacquez_result : dictionary
+ contains the statistic (stat) for the test and the
+ associated p-value (pvalue).
+ stat : float
+ value of the Jacquez k nearest neighbors test for the
+ dataset.
+ pvalue : float
+ p-value associated with the statistic (normally
+ distributed with k-1 df).
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import libpysal as lps
+ >>> from pointpats import SpaceTimeEvents, jacquez
+
+ Read in the example data and create an instance of SpaceTimeEvents.
+
+ >>> path = lps.examples.get_path("burkitt.shp")
+ >>> events = SpaceTimeEvents(path,'T')
+
+ The Jacquez test counts the number of events that are k nearest
+ neighbors in both time and space. The following runs the Jacquez test
+ on the example data and reports the resulting statistic. In this case,
+ there are 12 instances where events are nearest neighbors in both space
+ and time.
+ # turning off as kdtree changes from scipy < 0.12 return 13
+
+ >>> np.random.seed(100)
+ >>> result = jacquez(events.space, events.t ,k=3,permutations=99)
+ >>> print(result['stat'])
+ 12
+
+ The significance of this can be assessed by calling the p-
+ value from the results dictionary, as shown below. Again, no
+ space-time interaction is observed.
+
+ >>> result['pvalue'] < 0.01
+ False
+
+ """
+ time = t_coords
+ space = s_coords
+ n = len(time)
+
+ # calculate the nearest neighbors in space and time separately
+ knnt = lps.weights.KNN.from_array(time, k)
+ knns = lps.weights.KNN.from_array(space, k)
+
+ nnt = knnt.neighbors
+ nns = knns.neighbors
+ knn_sum = 0
+
+ # determine which events are nearest neighbors in both space and time
+ for i in range(n):
+ t_neighbors = nnt[i]
+ s_neighbors = nns[i]
+ check = set(t_neighbors)
+ inter = check.intersection(s_neighbors)
+ count = len(inter)
+ knn_sum += count
+
+ stat = knn_sum
+
+ # return the results (if no inference)
+ if not permutations:
+ return stat
+
+ # loop for generating a random distribution to assess significance
+ dist = []
+ for _p in range(permutations):
+ j = 0
+ trand = np.random.permutation(time)
+ knnt = lps.weights.KNN.from_array(trand, k)
+ nnt = knnt.neighbors
+ for i in range(n):
+ t_neighbors = nnt[i]
+ s_neighbors = nns[i]
+ check = set(t_neighbors)
+ inter = check.intersection(s_neighbors)
+ count = len(inter)
+ j += count
+
+ dist.append(j)
+
+ # establish the pseudo significance of the observed statistic
+ distribution = np.array(dist)
+ greater = np.ma.masked_greater_equal(distribution, stat)
+ count = np.ma.count_masked(greater)
+ pvalue = (count + 1.0) / (permutations + 1.0)
+
+ # report the results
+ jacquez_result = {"stat": stat, "pvalue": pvalue}
+ return jacquez_result
+
+
+
+
+[docs]
+def modified_knox(s_coords, t_coords, delta, tau, permutations=99):
+ """
+ Baker's modified Knox test for spatio-temporal interaction.
+ :cite:`Baker:2004`
+
+ Parameters
+ ----------
+ s_coords : array
+ (n, 2), spatial coordinates.
+ t_coords : array
+ (n, 1), temporal coordinates.
+ delta : float
+ threshold for proximity in space.
+ tau : float
+ threshold for proximity in time.
+ permutations : int, optional
+ the number of permutations used to establish pseudo-
+ significance (the default is 99).
+
+ Returns
+ -------
+ modknox_result : dictionary
+ contains the statistic (stat) for the test and the
+ associated p-value (pvalue).
+ stat : float
+ value of the modified knox test for the dataset.
+ pvalue : float
+ pseudo p-value associated with the statistic.
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import libpysal as lps
+ >>> from pointpats import SpaceTimeEvents, modified_knox
+
+ Read in the example data and create an instance of SpaceTimeEvents.
+
+ >>> path = lps.examples.get_path("burkitt.shp")
+ >>> events = SpaceTimeEvents(path, 'T')
+
+ Set the random seed generator. This is used by the permutation based
+ inference to replicate the pseudo-significance of our example results -
+ the end-user will normally omit this step.
+
+ >>> np.random.seed(100)
+
+ Run the modified Knox test with distance and time thresholds of 20 and 5,
+ respectively. This counts the events that are closer than 20 units in
+ space, and 5 units in time.
+
+ >>> result = modified_knox(events.space, events.t, delta=20, tau=5, permutations=99)
+
+ Next, we examine the results. First, we call the statistic from the
+ results dictionary. This reports the difference between the observed
+ and expected Knox statistic.
+
+ >>> print("%2.8f" % result['stat'])
+ 2.81016043
+
+ Next, we look at the pseudo-significance of this value, calculated by
+ permuting the timestamps and rerunning the statistics. In this case,
+ the results indicate there is likely no space-time interaction.
+
+ >>> print("%2.2f" % result['pvalue'])
+ 0.11
+
+ """
+ s = s_coords
+ t = t_coords
+ n = len(t)
+
+ # calculate the spatial and temporal distance matrices for the events
+ sdistmat = cg.distance_matrix(s)
+ tdistmat = cg.distance_matrix(t)
+
+ # identify events within thresholds
+ spacmat = np.ones((n, n))
+ spacbin = sdistmat <= delta
+ spacmat = spacmat * spacbin
+ timemat = np.ones((n, n))
+ timebin = tdistmat <= tau
+ timemat = timemat * timebin
+
+ # calculate the observed (original) statistic
+ knoxmat = timemat * spacmat
+ obsstat = knoxmat.sum() - n
+
+ # calculate the expectated value
+ ssumvec = np.reshape((spacbin.sum(axis=0) - 1), (n, 1))
+ tsumvec = np.reshape((timebin.sum(axis=0) - 1), (n, 1))
+ expstat = (ssumvec * tsumvec).sum()
+
+ # calculate the modified stat
+ stat = (obsstat - (expstat / (n - 1.0))) / 2.0
+
+ # return results (if no inference)
+ if not permutations:
+ return stat
+ distribution = []
+
+ # loop for generating a random distribution to assess significance
+ for _p in range(permutations):
+ rtdistmat = _shuffle_matrix(tdistmat, list(range(n)))
+ timemat = np.ones((n, n))
+ timebin = rtdistmat <= tau
+ timemat = timemat * timebin
+
+ # calculate the observed knox again
+ knoxmat = timemat * spacmat
+ obsstat = knoxmat.sum() - n
+
+ # calculate the expectated value again
+ ssumvec = np.reshape((spacbin.sum(axis=0) - 1), (n, 1))
+ tsumvec = np.reshape((timebin.sum(axis=0) - 1), (n, 1))
+ expstat = (ssumvec * tsumvec).sum()
+
+ # calculate the modified stat
+ tempstat = (obsstat - (expstat / (n - 1.0))) / 2.0
+ distribution.append(tempstat)
+
+ # establish the pseudo significance of the observed statistic
+ distribution = np.array(distribution)
+ greater = np.ma.masked_greater_equal(distribution, stat)
+ count = np.ma.count_masked(greater)
+ pvalue = (count + 1.0) / (permutations + 1.0)
+
+ # return results
+ modknox_result = {"stat": stat, "pvalue": pvalue}
+ return modknox_result
+
+
+
+def _shuffle_matrix(X, ids):
+ """
+ Random permutation of rows and columns of a matrix
+
+ Parameters
+ ----------
+ X : array
+ (k, k), array to be permuted.
+ ids : array
+ range (k, ).
+
+ Returns
+ -------
+ : array
+ (k, k) with rows and columns randomly shuffled.
+ """
+ np.random.shuffle(ids)
+ return X[ids, :][:, ids]
+
+
+def _knox(s_coords, t_coords, delta, tau, permutations=99, keep=False):
+ """
+ Parameters
+ ==========
+
+ s_coords: array-like
+ spatial coordinates
+ t_coords: array-like
+ temporal coordinates
+ delta: float
+ distance threshold
+ tau: float
+ temporal threshold
+ permutations: int
+ number of permutations
+ keep: bool
+ return values from permutations (default False)
+
+
+ Returns
+ =======
+
+ summary table observed
+ summary table h0
+
+ ns
+ nt
+ nst
+ n
+ p-value
+ """
+
+ n = s_coords.shape[0]
+
+ stree = KDTree(s_coords)
+ ttree = KDTree(t_coords)
+ sneighbors = stree.query_ball_tree(stree, r=delta)
+ sneighbors = [
+ set(neighbors).difference([i]) for i, neighbors in enumerate(sneighbors)
+ ]
+ tneighbors = ttree.query_ball_tree(ttree, r=tau)
+ tneighbors = [
+ set(neighbors).difference([i]) for i, neighbors in enumerate(tneighbors)
+ ]
+
+ # number of spatial neighbor pairs
+ ns = np.array([len(neighbors) for neighbors in sneighbors]) # by i
+
+ NS = ns.sum() / 2 # total
+
+ # number of temporal neigbor pairs
+ nt = np.array([len(neighbors) for neighbors in tneighbors])
+ NT = nt.sum() / 2
+
+ # s-t neighbors (list of lists)
+ stneighbors = [
+ sneighbors_i.intersection(tneighbors_i)
+ for sneighbors_i, tneighbors_i in zip(sneighbors, tneighbors)
+ ]
+
+ # number of spatio-temporal neigbor pairs
+ nst = np.array([len(neighbors) for neighbors in stneighbors])
+ NST = nst.sum() / 2
+
+ all_pairs = []
+ pairs = {}
+ for i, neigh in enumerate(stneighbors):
+ if len(neigh) > 0:
+ all_pairs.extend([sorted((i, j)) for j in neigh])
+ st_pairs = {tuple(l) for l in all_pairs}
+
+ # ENST: expected number of spatio-temporal neighbors under HO
+ pairs = n * (n - 1) / 2
+ ENST = NS * NT / pairs
+
+ # observed table
+ observed = np.zeros((2, 2))
+
+ NS_ = NS - NST # spatial only
+ NT_ = NT - NST # temporal only
+
+ observed[0, 0] = NST
+ observed[0, 1] = NS_
+ observed[1, 0] = NT_
+ observed[1, 1] = pairs - NST - NS_ - NT_
+
+ # expected table
+
+ expected = np.zeros((2, 2))
+ expected[0, 0] = ENST
+ expected[0, 1] = NS - expected[0, 0]
+ expected[1, 0] = NT - expected[0, 0]
+ expected[1, 1] = pairs - expected.sum()
+
+ p_value_poisson = 1 - poisson.cdf(NST, expected[0, 0])
+
+ results = {}
+ results["ns"] = ns.sum() / 2
+ results["nt"] = nt.sum() / 2
+ results["nst"] = nst.sum() / 2
+ results["pairs"] = pairs
+ results["expected"] = expected
+ results["observed"] = observed
+ results["p_value_poisson"] = p_value_poisson
+ results["st_pairs"] = st_pairs
+ results["sneighbors"] = sneighbors
+ results["tneighbors"] = tneighbors
+ results["stneighbors"] = stneighbors
+
+ if permutations > 0:
+ exceedence = 0
+ n = len(sneighbors)
+ ids = np.arange(n)
+ if keep:
+ ST = np.zeros(permutations)
+
+ for perm in range(permutations):
+ st = 0
+ rids = np.random.permutation(ids)
+ for i in range(n):
+ ri = rids[i]
+ tni = tneighbors[ri]
+ rjs = [rids[j] for j in sneighbors[i]]
+ sti = [j for j in rjs if j in tni]
+ st += len(sti)
+ st /= 2
+ if st >= results["nst"]:
+ exceedence += 1
+ if keep:
+ ST[perm] = st
+ results["p_value_sim"] = (exceedence + 1) / (permutations + 1)
+ results["exceedence"] = exceedence
+ if keep:
+ results["st_perm"] = ST
+
+ return results
+
+
+
+[docs]
+class Knox:
+ """Global Knox statistic for space-time interactions
+
+ Parameters
+ ----------
+ s_coords: array-like
+ spatial coordinates of point events
+ t_coords: array-like
+ temporal coordinates of point events (floats or ints, not dateTime)
+ delta: float
+ spatial threshold defining distance below which pairs are spatial
+ neighbors
+ tau: float
+ temporal threshold defining distance below which pairs are temporal
+ neighbors
+ permutations: int
+ number of random permutations for inference
+ keep: bool
+ whether to store realized values of the statistic under permutations
+
+
+ Attributes
+ ----------
+ s_coords: array-like
+ spatial coordinates of point events
+ t_coords: array-like
+ temporal coordinates of point events (floats or ints, not dateTime)
+ delta: float
+ spatial threshold defining distance below which pairs are spatial
+ neighbors
+ tau: float
+ temporal threshold defining distance below which pairs are temporal
+ neighbors
+ permutations: int
+ number of random permutations for inference
+ keep: bool
+ whether to store realized values of the statistic under permutations
+ nst: int
+ number of space-time pairs
+ p_poisson: float
+ Analytical p-value under Poisson assumption
+ p_sim: float
+ Pseudo p-value based on random permutations
+ expected: array
+ Two-by-two array with expected counts under the null of no space-time
+ interactions. [[NST, NS_], [NT_, N__]] where NST is the expected number
+ of space-time pairs, NS_ is the expected number of spatial (but not also
+ temporal) pairs, NT_ is the number of expected temporal (but not also
+ spatial pairs), N__ is the number of pairs that are neighor spatial or
+ temporal neighbors.
+ observed: array
+ Same structure as expected with the observed pair classifications
+ sim: array
+ Global statistics from permutations (if keep=True)
+
+ Notes
+ -----
+ Technical details can be found in :cite:`Rogerson:2001`
+
+ Examples
+ --------
+ >>> import libpysal
+ >>> path = libpysal.examples.get_path('burkitt.shp')
+ >>> import geopandas
+ >>> df = geopandas.read_file(path)
+ >>> from pointpats.spacetime import Knox
+ >>> global_knox = Knox(df[['X', 'Y']], df[["T"]], delta=20, tau=5)
+ >>> global_knox.statistic_
+ 13
+ >>> global_knox.p_poisson
+ 0.14624558197140414
+ >>> global_knox.observed
+ array([[1.300e+01, 3.416e+03],
+ [3.900e+01, 1.411e+04]])
+ >>> global_knox.expected
+ array([[1.01438161e+01, 3.41885618e+03],
+ [4.18561839e+01, 1.41071438e+04]])
+ >>> hasattr(global_knox, 'sim')
+ False
+ >>> import numpy
+ >>> numpy.random.seed(12345)
+ >>> global_knox = Knox(df[['X', 'Y']], df[["T"]], delta=20, tau=5, keep=True)
+ >>> hasattr(global_knox, 'sim')
+ True
+ >>> global_knox.p_sim
+ 0.21
+ """
+
+
+[docs]
+ def __init__(self, s_coords, t_coords, delta, tau, permutations=99, keep=False):
+ self.s_coords = s_coords
+ self.t_coords = t_coords
+ self.delta = delta
+ self.tau = tau
+ self.permutations = permutations
+ self.keep = keep
+ results = _knox(s_coords, t_coords, delta, tau, permutations, keep)
+ self.nst = int(results["nst"])
+ if permutations > 0:
+ self.p_sim = results["p_value_sim"]
+ if keep:
+ self.sim = results["st_perm"]
+
+ self.p_poisson = results["p_value_poisson"]
+ self.observed = results["observed"]
+ self.expected = results["expected"]
+ self.statistic_ = self.nst
+
+
+
+[docs]
+ @classmethod
+ def from_dataframe(
+ cls,
+ dataframe,
+ time_col: int,
+ delta: int,
+ tau: int,
+ permutations: int = 99,
+ keep: bool = False,
+ ):
+ """Compute a Knox statistic from a dataframe of Point observations
+
+ Parameters
+ ----------
+ dataframe : geopandas.GeoDataFrame
+ geodataframe holding observations. Should be in a projected coordinate
+ system with geometries stored as Point
+ time_col : str
+ column in the dataframe storing the time values (integer coordinate)
+ for each observation. For example if the observations are stored with
+ a timestamp, the time_col should be converted to a series of integers
+ representing, e.g. hours, days, seconds, etc.
+ delta : int
+ delta parameter defining the spatial neighbor threshold measured in the
+ same units as the dataframe CRS
+ tau : int
+ tau parameter defining the temporal neihgbor threshold (in the units
+ measured by `time_col`)
+ permutations : int, optional
+ permutations to use for computation inference, by default 99
+ keep : bool
+ whether to store realized values of the statistic under permutations
+
+ Returns
+ -------
+ pointpats.spacetime.Knox
+ a fitted Knox class
+
+ """
+ s_coords, t_coords = _spacetime_points_to_arrays(dataframe, time_col)
+
+ return cls(s_coords, dataframe[[time_col]], delta, tau, permutations, keep)
+
+
+
+
+def _knox_local(s_coords, t_coords, delta, tau, permutations=99, keep=False):
+ """
+
+ Parameters
+ ----------
+ s_coords: array (n,2)
+ spatial coordinates
+ t_coords: array (n,1)
+ temporal coordinates
+ delta: numeric
+ spatial threshold distance for neighbor relation
+ tau: numeric
+ temporal threshold distance for neighbor relation
+ permutations: int
+ number of permutations for conditional randomization inference
+ keep: bool
+ whether to store local statistics from the permtuations
+
+ """
+ # think about passing in the global object as an option to avoid recomputing the trees
+ res = _knox(s_coords, t_coords, delta, tau, permutations=permutations)
+ sneighbors = {i: tuple(ns) for i, ns in enumerate(res["sneighbors"])}
+ tneighbors = {i: tuple(nt) for i, nt in enumerate(res["tneighbors"])}
+
+ n = len(s_coords)
+ ids = np.arange(n)
+ res["nsti"] = np.zeros(n) # number of observed st_pairs for observation i
+ res["nsi"] = [len(r) for r in res["sneighbors"]]
+ res["nti"] = [len(r) for r in res["tneighbors"]]
+ for pair in res["st_pairs"]:
+ i, j = pair
+ res["nsti"][i] += 1
+ res["nsti"][j] += 1
+
+ nsti = res["nsti"]
+ nsi = res["nsi"]
+ res["nti"]
+
+ # rather than do n*permutations, we reuse the permutations
+ # ensuring that each permutation is conditional on a focal unit i
+ # for each of the permutations we loop over i and swap labels between the
+ # label at index i in the current permutation and the label at the index
+ # assigned i in the permutation.
+
+ if permutations > 0:
+ exceedence = np.zeros(n)
+ if keep:
+ STI = np.zeros((n, permutations))
+ for perm in range(permutations):
+ rids = np.random.permutation(ids)
+ for i in range(n):
+ rids_i = rids.copy()
+ # set observed value of focal unit i
+ # swap with value assigned to rids[i]
+ # example
+ # 0 1 2 (ids)
+ # 2 0 1 (rids)
+ # i=0
+ # 0 2 1 (rids_i)
+ # i=1
+ # 2 1 0 (rids_i)
+ # i=2
+ # 1 0 2 (rids_i)
+
+ rids_i[rids == i] = rids[i]
+ rids_i[i] = i
+
+ # calculate local stat
+ rjs = [rids_i[j] for j in sneighbors[i]]
+ tni = tneighbors[i]
+ sti = [j for j in rjs if j in tni]
+ count = len(sti)
+ if count >= res["nsti"][i]:
+ exceedence[i] += 1
+ if keep:
+ STI[i, perm] = count
+
+ if keep:
+ res["sti_perm"] = STI
+ res["exceedence_pvalue"] = (exceedence + 1) / (permutations + 1)
+ res["exceedences"] = exceedence
+
+ # analytical inference
+ ntjis = [len(r) for r in res["tneighbors"]]
+ n1 = n - 1
+ hg_pvalues = [
+ 1 - hypergeom.cdf(nsti[i] - 1, n1, ntjis, nsi[i]).mean() for i in range(n)
+ ]
+ res["hg_pvalues"] = np.array(hg_pvalues)
+
+ # identification of hot spots
+
+ adjlist = []
+ for i, j in res["st_pairs"]:
+ adjlist.append([i, j])
+ adjlist.append([j, i])
+ adjlist = pandas.DataFrame(data=adjlist, columns=["focal", "neighbor"])
+ adjlist = adjlist.sort_values(by=["focal", "neighbor"])
+ adjlist.reset_index(drop=True, inplace=True)
+
+ adjlist["orientation"] = ""
+ for index, row in adjlist.iterrows():
+ focal = row["focal"]
+ neighbor = row["neighbor"]
+ ft = t_coords[focal]
+ nt = t_coords[neighbor]
+ if ft < nt:
+ adjlist.iloc[index, 2] = "lead"
+ elif ft > nt:
+ adjlist.iloc[index, 2] = "lag"
+ else:
+ adjlist.iloc[index, 2] = "coincident"
+
+ res["stadjlist"] = adjlist
+ return res
+
+
+
+[docs]
+class KnoxLocal:
+ """Local Knox statistics for space-time interactions
+
+ Parameters
+ ----------
+ s_coords: array (nx2)
+ spatial coordinates of point events
+ t_coords: array (nx1)
+ temporal coordinates of point events (floats or ints, not dateTime)
+ delta: float
+ spatial threshold defining distance below which pairs are spatial
+ neighbors
+ tau: float
+ temporal threshold defining distance below which pairs are temporal
+ neighbors
+ permutations: int
+ number of random permutations for inference
+ keep: bool
+ whether to store realized values of the statistic under permutations
+ conditional: bool
+ whether to include conditional permutation inference
+ crit: float
+ signifcance level for local statistics
+ crs: str (optional)
+ coordinate reference system string for s_coords
+
+ Attributes
+ ----------
+ s_coords: array (nx2)
+ spatial coordinates of point events
+ t_coords: array (nx1)
+ temporal coordinates of point events (floats or ints, not dateTime)
+ delta: float
+ spatial threshold defining distance below which pairs are spatial
+ neighbors
+ tau: float
+ temporal threshold defining distance below which pairs are temporal
+ neighbors
+ permutations: int
+ number of random permutations for inference
+ keep: bool
+ whether to store realized values of the statistic under permutations
+ nst: int
+ number of space-time pairs (global)
+ p_poisson: float
+ Analytical p-value under Poisson assumption (global)
+ p_sim: float
+ Pseudo p-value based on random permutations (global)
+ expected: array
+ Two-by-two array with expected counts under the null of no space-time
+ interactions. [[NST, NS_], [NT_, N__]] where NST is the expected number
+ of space-time pairs, NS_ is the expected number of spatial (but not also
+ temporal) pairs, NT_ is the number of expected temporal (but not also
+ spatial pairs), N__ is the number of pairs that are neighor spatial or
+ temporal neighbors. (global)
+ observed: array
+ Same structure as expected with the observed pair classifications (global)
+ sim: array
+ Global statistics from permutations (if keep=True and keep=True) (global)
+ p_sims: array
+ Local psuedo p-values from conditional permutations (if permutations>0)
+ sims: array
+ Local statistics from conditional permutations (if keep=True and
+ permutations>0)
+ nsti: array
+ Local statistics
+ p_hypergeom: array
+ Analytical p-values based on hypergeometric distribution
+
+ Notes
+ -----
+ Technical details can be found in :cite:`Rogerson:2001`. The conditional
+ permutation inference is unique to pysal.pointpats.
+
+ Examples
+ -------
+ >>> import libpysal
+ >>> path = libpysal.examples.get_path('burkitt.shp')
+ >>> import geopandas
+ >>> df = geopandas.read_file(path)
+ >>> from pointpats.spacetime import Knox
+ >>> import numpy
+ >>> numpy.random.seed(12345)
+ >>> local_knox = KnoxLocal(df[['X', 'Y']], df[["T"]], delta=20, tau=5, keep=True)
+ >>> local_knox.statistic_.shape
+ (188,)
+ >>> lres = local_knox
+ >>> gt0ids = numpy.where(lres.nsti>0)
+ >>> gt0ids # doctest: +NORMALIZE_WHITESPACE
+ (array([ 25, 26, 30, 31, 35, 36, 41, 42, 46, 47, 51, 52, 102,
+ 103, 116, 118, 122, 123, 137, 138, 139, 140, 158, 159, 162, 163]),)
+ >>> lres.nsti[gt0ids]
+ array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
+ 1., 1., 1., 1., 1., 1., 1., 1., 1.])
+ >>> lres.p_hypergeom[gt0ids]
+ array([0.1348993 , 0.14220663, 0.07335085, 0.08400282, 0.1494317 ,
+ 0.21524073, 0.0175806 , 0.04599869, 0.17523687, 0.18209188,
+ 0.19111321, 0.16830444, 0.13734428, 0.14703242, 0.06796364,
+ 0.03192559, 0.13734428, 0.17523687, 0.12998154, 0.1933476 ,
+ 0.13244507, 0.13244507, 0.12502644, 0.14703242, 0.12502644,
+ 0.12998154])
+ >>> lres.p_sims[gt0ids]
+ array([0.3 , 0.33, 0.11, 0.17, 0.3 , 0.42, 0.06, 0.06, 0.33, 0.34, 0.36,
+ 0.38, 0.3 , 0.29, 0.41, 0.19, 0.31, 0.39, 0.18, 0.39, 0.48, 0.41,
+ 0.22, 0.41, 0.39, 0.32])
+ """
+
+
+[docs]
+ def __init__(
+ self,
+ s_coords,
+ t_coords,
+ delta,
+ tau,
+ permutations=99,
+ keep=False,
+ crit=0.05,
+ crs=None,
+ ids=None,
+ ):
+ if not isinstance(t_coords, np.ndarray):
+ raise ValueError("t_coords should be numpy.ndarray type")
+ if not isinstance(s_coords, np.ndarray):
+ raise ValueError("s_coords should be numpy.ndarray type")
+ n_s, k = s_coords.shape
+ rangeids = list(range(n_s))
+ if k < 2:
+ raise ValueError("s_coords shape required to be nx2")
+ n_t, k = t_coords.shape
+ if n_s != n_t:
+ raise ValueError("t_coords and s_coords need to be same length")
+ if ids is not None:
+ if len(ids) != n_s:
+ raise ValueError("`ids` must have the same length as the inputs")
+ else:
+ ids = rangeids
+
+ self._ids = ids
+ self.s_coords = s_coords
+ self.t_coords = t_coords
+ self.delta = delta
+ self.tau = tau
+ self.permutations = permutations
+ self.keep = keep
+ self.crit = crit
+ results = _knox_local(s_coords, t_coords, delta, tau, permutations, keep)
+ self.adjlist = results["stadjlist"]
+ self.nst = int(results["nst"])
+ if permutations > 0:
+ self.p_sim = results["p_value_sim"]
+ if keep:
+ self.sim = results["sti_perm"]
+
+ self.p_poisson = results["p_value_poisson"]
+ self.observed = results["observed"]
+ self.expected = results["expected"]
+ self.p_hypergeom = results["hg_pvalues"]
+ if permutations > 0:
+ self.p_sims = results["exceedence_pvalue"]
+ if keep:
+ self.sims = results["sti_perm"]
+ self.nsti = results["nsti"]
+ # self.hotspots = results["hotspots"]
+ self._crs = crs
+ self.statistic_ = self.nsti
+ self._id_map = dict(zip(rangeids, self._ids))
+ self.adjlist["focal"] = self.adjlist["focal"].replace(self._id_map)
+ self.adjlist["neighbor"] = self.adjlist["neighbor"].replace(self._id_map)
+ # reconstruct df
+ geom = gpd.points_from_xy(self.s_coords[:, 0], self.s_coords[:, 1])
+ _gdf = gpd.GeoDataFrame(geometry=geom, crs=self._crs, index=self._ids)
+ _gdf["time"] = self.t_coords
+ if permutations > 0:
+ _gdf["p_sim"] = self.p_sims
+ _gdf["p_hypergeom"] = self.p_hypergeom
+ self._gdf_static = _gdf
+
+
+ @property
+ def _gdf(self):
+ return self._gdf_static.copy()
+
+
+[docs]
+ @classmethod
+ def from_dataframe(
+ cls,
+ dataframe,
+ time_col: str,
+ delta: int,
+ tau: int,
+ permutations: int = 99,
+ keep: bool = False,
+ ):
+ """Compute a set of local Knox statistics from a dataframe of Point observations
+
+ Parameters
+ ----------
+ dataframe : geopandas.GeoDataFrame
+ dataframe holding observations. Should be in a projected coordinate system
+ with geometries stored as Points
+ time_col : str
+ column in the dataframe storing the time values (integer coordinate)
+ for each observation. For example if the observations are stored with
+ a timestamp, the time_col should be converted to a series of integers
+ representing, e.g. hours, days, seconds, etc.
+ delta : int
+ delta parameter defining the spatial neighbor threshold measured in the
+ same units as the dataframe CRS
+ tau : int
+ tau parameter defining the temporal neihgbor threshold (in the units
+ measured by `time_col`)
+ permutations : int, optional
+ permutations to use for computational inference, by default 99
+ keep : bool
+ whether to store realized values of the statistic under permutations
+
+ Returns
+ -------
+ pointpats.spacetime.LocalKnox
+ a fitted KnoxLocal class
+
+ """
+ s_coords, t_coords = _spacetime_points_to_arrays(dataframe, time_col)
+
+ return cls(
+ s_coords,
+ t_coords,
+ delta,
+ tau,
+ permutations,
+ keep,
+ crs=dataframe.crs,
+ ids=dataframe.index.values,
+ )
+
+
+
+[docs]
+ def hotspots(self, crit=0.05, inference="permutation", keep_neighbors=True):
+ """Table of significant space-time clusters that define local hotspots.
+
+ Parameters
+ ----------
+ crit : float, optional
+ critical value for statistical inference, by default 0.05
+ inference : str, optional
+ whether p-values should use permutation or analutical inference, by default
+ "permutation"
+ keep_neighbors: bool
+ whether to included nonsignificant members of hotspots. While these
+ observations are not themselves significant, these still define the spatial
+ extent of the cluster, and the the focal observation cannot become
+ significant without their presence. If True, return all members of a
+ significant hotspot, else return only the significant locations
+
+ Returns
+ -------
+ pandas.DataFrame
+ dataframe of significant hotspots
+
+ Raises
+ ------
+ ValueError
+ if `inference` is not in {'permutation', 'analytic'}
+ """
+ if inference == "permutation":
+ if not hasattr(self, "p_sim"):
+ warn(
+ "Pseudo-p values not availalable. Permutation-based p-values require "
+ "fitting the KnoxLocal class using `permutations` set to a large "
+ "number. Using analytic p-values instead",
+ stacklevel=1,
+ )
+ col = "p_hypergeom"
+ else:
+ col = "p_sim"
+ elif inference == "analytic":
+ col = "p_hypergeom"
+ else:
+ raise ValueError("inference must be either `permutation` or `analytic`")
+ # determine hot spots
+ pdf_sig = self._gdf[self._gdf[col] <= crit][[col, "time"]].rename(
+ columns={col: "pvalue", "time": "focal_time"}
+ )
+
+ # if keep_neighbors, we want to include a 'cluster' column denoting which
+ # cluster nonsig observations belong to. Need to use a graph for that
+ temp_neighbors = self.adjlist[
+ (self.adjlist.focal.isin(pdf_sig.index.values))
+ | self.adjlist.neighbor.isin(pdf_sig.index.values)
+ ]
+ pdf_sig = pd.concat([pdf_sig,
+ self._gdf[self._gdf.index.isin(temp_neighbors.neighbor.values)
+
+ ][[col, "time"]].rename(
+ columns={col: "pvalue", "time": "focal_time"}
+ )
+ ])
+
+ pdf_sig = pdf_sig.merge(
+ temp_neighbors, how='outer', left_index=True, right_on="focal"
+ ).reset_index(drop=True)
+ # significant focals can be neighbors of others (dupes)
+ pdf_sig = pdf_sig.groupby("focal").first().reset_index()
+ graph = Graph.from_adjacency(pdf_sig.assign(weight=1))
+ pdf_sig["cluster"] = graph.component_labels.values
+ if not keep_neighbors :
+ pdf_sig = pdf_sig[pdf_sig.pvalue<=crit]
+
+ return self._gdf[["geometry"]].merge(
+ pdf_sig.copy(), left_index=True, right_on="focal"
+ )
+
+
+
+[docs]
+ def plot(
+ self,
+ colors: dict = {"focal": "red", "neighbor": "yellow", "nonsig": "grey"},
+ crit: float = 0.05,
+ inference: str = "permutation",
+ point_kwargs: dict = None,
+ plot_edges: bool = True,
+ edge_color: str = "black",
+ edge_kwargs: dict = None,
+ ax=None,
+ ):
+ """plot hotspots
+
+ Parameters
+ ----------
+ colors : dict, optional
+ mapping of colors to hotspot values, by default
+ {"focal": "red", "neighbor": "yellow", "nonsig": "grey"}
+ crit : float, optional
+ critical value for assessing statistical sgifnicance, by default 0.05
+ inference : str, optional
+ whether to use permutation or analytic inference, by default "permutation"
+ point_kwargs : dict, optional
+ additional keyword arguments passsed to point plot, by default None
+ plot_edges : bool, optional
+ whether to plot edges connecting members of the same hotspot subgraph,
+ by default True
+ edge_color : str, optional
+ color of edges when plot_edges is True, by default 'black'
+ edge_kwargs : dict, optional
+ additional keyword arguments passsed to edge plot, by default None
+ ax : matplotlib.axes.Axes, optional
+ axes object on which to create the plot, by default None
+
+ Returns
+ -------
+ matplotlib.axes.Axes
+ plot of local space-time hotspots
+ """
+
+ if point_kwargs is None:
+ point_kwargs = dict()
+ if edge_kwargs is None:
+ edge_kwargs = dict()
+ g = self._gdf.copy()
+
+ g["color"] = colors["nonsig"]
+ g["pvalue"] = self.p_hypergeom
+ if inference == "permutation":
+ if not hasattr(self, "p_sims"):
+ warn(
+ "Pseudo-p values not availalable. Permutation-based p-values require "
+ "fitting the KnoxLocal class using `permutations` set to a large "
+ "number. Using analytic p-values instead"
+ )
+ g["pvalue"] = self.p_hypergeom
+ else:
+ g["pvalue"] = self.p_sims
+ elif inference == "analytic":
+ g["pvalue"] = self.p_hypergeom
+ else:
+ raise ValueError("inference must be either `permutation` or `analytic`")
+
+ mask = g[g.pvalue <= crit].index.values
+ neighbors = self.adjlist[self.adjlist.focal.isin(mask)].neighbor.unique()
+ g.loc[neighbors, "color"] = colors["neighbor"]
+ g.loc[g.pvalue <= crit, "color"] = colors["focal"]
+ m = g[g.color == colors["nonsig"]].plot(
+ color=colors["nonsig"], ax=ax, **point_kwargs
+ )
+ g[g.color == colors["neighbor"]].plot(
+ ax=m, color=colors["neighbor"], **point_kwargs
+ )
+ g[g.color == colors["focal"]].plot(ax=m, color=colors["focal"], **point_kwargs)
+
+ if plot_edges:
+ # edges between hotspot and st-neighbors
+ ghs = self.hotspots(crit=crit, inference=inference)
+ ghs = ghs.dropna()
+ origins = g.loc[ghs.focal].geometry
+ destinations = g.loc[ghs.neighbor].geometry
+ ods = zip(origins, destinations)
+ lines = gpd.GeoSeries([LineString(od) for od in ods], crs=g.crs)
+ lines.plot(ax=m, color=edge_color, **edge_kwargs)
+
+ return m
+
+
+
+[docs]
+ def explore(
+ self,
+ crit: float = 0.05,
+ inference: str = "permutation",
+ radius: int = 5,
+ style_kwds: dict = None,
+ tiles: str = "CartoDB Positron",
+ plot_edges: bool = True,
+ edge_weight: int = 2,
+ edge_color: str = "black",
+ colors: dict = {"focal": "red", "neighbor": "yellow", "nonsig": "grey"},
+ ):
+ """Interactive plotting for space-time hotspots.
+
+ Parameters
+ ----------
+ crit : float, optional
+ critical value for statistical inference, by default 0.05
+ inference : str, optional
+ which p-value to use for determining hotspots. Either "permutation" or
+ "analytic", by default "permutation"
+ radius : int, optional
+ radius of the circlemarker plotted by folium, passed to
+ geopandas.GeoDataFrame.explore style_kwds as a convenience. Ignored if
+ `style_kwds` is passed directly, by default 5
+ style_kwds : dict, optional
+ additional style kewords passed to GeoDataFrame.explore, by default None
+ tiles : str, optional
+ tileset passed to GeoDataFrame.explore `tiles` argument, by default
+ "CartoDB Positron"
+ plot_edges : bool, optional
+ Whether to include lines drawn between members of a singnificant hotspot, by
+ default True
+ edge_weight : int, optional
+ line thickness when `plot_edges=True`, by default 2
+ edge_color : str, optional
+ color of line when `plot_edges=True`, by default "black"
+ colors : dict, optional
+ mapping of observation type to color,
+ by default {"focal": "red", "neighbor": "yellow", "nonsig": "grey"}
+
+ Returns
+ -------
+ folium.Map
+ an interactive map showing locally-significant spacetime hotspots
+ """
+ if style_kwds is None:
+ style_kwds = {"radius": radius}
+ g = self._gdf.copy()
+
+ g["color"] = colors["nonsig"]
+
+ if inference == "permutation":
+ if not hasattr(self, "p_sims"):
+ warn(
+ "Pseudo-p values not availalable. Permutation-based p-values require "
+ "fitting the KnoxLocal class using `permutations` set to a large "
+ "number. Using analytic p-values instead"
+ )
+ g["pvalue"] = self.p_hypergeom
+ else:
+ g["pvalue"] = self.p_sims
+ elif inference == "analytic":
+ g["pvalue"] = self.p_hypergeom
+ else:
+ raise ValueError("inference must be either `permutation` or `analytic`")
+
+ mask = g[g.pvalue <= crit].index.values
+ neighbors = self.adjlist[self.adjlist.focal.isin(mask)].neighbor.unique()
+
+ # this is clunky, but enforces plotting order so significance is prioritized
+ g.loc[neighbors, "color"] = colors["neighbor"]
+ g.loc[g.pvalue <= crit, "color"] = colors["focal"]
+ nbs = self.adjlist.groupby("focal").agg(list)["neighbor"]
+ g = g.reset_index().merge(nbs, left_on="index", right_index=True, how="left")
+
+ m = g[g.color == colors["nonsig"]].explore(
+ color="grey", style_kwds=style_kwds, tiles=tiles
+ )
+ blues = g[g.color == colors["neighbor"]]
+ if blues.shape[0] == 0:
+ warn("empty neighbor set.")
+ else:
+ m = blues.explore(m=m, color=colors["neighbor"], style_kwds=style_kwds)
+ m = g[g.color == colors["focal"]].explore(
+ m=m, color=colors["focal"], style_kwds=style_kwds
+ )
+
+ if plot_edges:
+ # edges between hotspot and st-neighbors
+ g = g.set_index("index")
+ ghs = self.hotspots(crit=crit, inference=inference)
+ ghs = ghs.dropna()
+ origins = g.loc[ghs.focal].geometry
+ destinations = g.loc[ghs.neighbor].geometry
+ ods = zip(origins, destinations)
+ lines = gpd.GeoSeries([LineString(od) for od in ods], crs=g.crs)
+ lines.explore(m=m, color=edge_color, style_kwds={"weight": edge_weight})
+
+ return m
+
+
+ def _gdfhs(self, crit=0.05, inference="permutation"):
+ # merge df with self.hotspots
+ return gpd.GeoDataFrame(
+ self._gdf.merge(
+ self.hotspots(crit=crit, inference=inference),
+ left_index=True,
+ right_on="focal",
+ )
+ )
+
+
+
+def _spacetime_points_to_arrays(dataframe, time_col):
+ """convert long-form geodataframe into arrays for kdtree
+
+ Parameters
+ ----------
+ dataframe : geopandas.GeoDataFrame
+ geodataframe with point geometries
+ time_col : str
+ name of the column on dataframe that stores time values
+
+ Returns
+ -------
+ tuple
+ two numpy arrays holding spatial coodinates s_coords (n,2)
+ and temporal coordinates t_coords (n,1)
+
+ """
+ if dataframe.crs is None:
+ warn(
+ "There is no CRS set on the dataframe. The KDTree will assume coordinates "
+ "are stored in Euclidean distances"
+ )
+ else:
+ if dataframe.crs.is_geographic:
+ raise ValueError(
+ "The input dataframe must be in a projected coordinate system."
+ )
+
+ assert dataframe.geom_type.unique().tolist() == [
+ "Point"
+ ], "The Knox statistic is only defined for Point geometries"
+
+ # kdtree wont operate on datetime
+ if is_numeric_dtype(dataframe[time_col].dtype) is False:
+ raise ValueError(
+ "The time values must be stored as "
+ f"a numeric dtype but the column {time_col} is stored as "
+ f"{dataframe[time_col].dtype}"
+ )
+
+ s_coords = np.vstack((dataframe.geometry.x.values, dataframe.geometry.y.values)).T
+ t_coords = np.vstack(dataframe[time_col].values)
+
+ return s_coords, t_coords
+
+"""
+Window class for point patterns
+"""
+
+__author__ = "Serge Rey sjsrey@gmail.com"
+
+import warnings
+
+import libpysal as ps
+import numpy as np
+
+__all__ = ["as_window", "poly_from_bbox", "to_ccf", "Window"]
+
+
+warnings.filterwarnings(
+ "ignore", "Objects based on the `Geometry` class will", FutureWarning
+)
+
+
+
+[docs]
+def poly_from_bbox(bbox):
+ l, b, r, t = bbox
+ c = [(l, b), (l, t), (r, t), (r, b), (l, b)]
+ return ps.cg.shapes.Polygon(c)
+
+
+
+
+
+
+
+
+[docs]
+def as_window(pysal_polygon):
+ """
+ Convert a libpysal polygon to a Window.
+
+ Parameters
+ ----------
+ pysal_polygon: libpysal.cg.shapes.Polygon
+ libpysal Polygon instance.
+
+ Returns
+ -------
+ Window
+ A Window instance.
+ """
+
+ if pysal_polygon.holes == [[]]:
+ return Window(pysal_polygon.parts)
+ else:
+ return Window(pysal_polygon.parts, pysal_polygon.holes)
+
+
+
+
+[docs]
+class Window(ps.cg.Polygon):
+ """
+ Geometric container for point patterns.
+
+ A window is used to define the area over which the pattern is observed.
+ This area is used in estimating the intensity of the point pattern.
+ See :attr:`PointPattern.lambda_window`.
+
+ Parameters
+ ----------
+ parts: sequence
+ A sequence of rings which bound the positive space point
+ pattern.
+ holes: sequence
+ A sequence of rings which bound holes in the polygons that bound the
+ point pattern.
+
+ """
+
+
+[docs]
+ def __init__(self, parts, holes=[]):
+ if holes:
+ super(Window, self).__init__(parts, holes)
+ else:
+ super(Window, self).__init__(parts)
+
+
+
+[docs]
+ def filter_contained(self, points):
+ return [np.asarray(pnt) for pnt in points if self.contains_point(pnt)]
+
+
+
+ 
+ 
+