<jupyter_start><jupyter_text>Python for spatial analysis  Part 3 - Point pattern analysis - Distance Based Statistical Method (G function, F function, Ripley K)- Point simulationI recommend an open-source python library for the statistical analysis of point patterns - [pointpats](https://pointpats.readthedocs.io/en/master/point-pattern-analysis-pointpats). Although it's a subpackage of PySAL, you need to install it specially via 'pip install pointpats'<jupyter_code>from pointpats import PointPattern, distance_statistics, PoissonPointProcess
import geopandas as gpd<jupyter_output><empty_output><jupyter_text>1. Statistical analysis<jupyter_code>point_data = gpd.read_file('../Data/BJ6Ring_POI_point.shp',encoding='utf-8')
point_data = point_data.sample(1000,random_state=1)
pattern = PointPattern(point_data, coord_names=['wgs84_x','wgs84_y'])
# The first parameter of PointPattern can be (n,p) array, n points with p>=2 attributes. 
# The default is that the first two columns are x and y coordinates.
# Or, you can the name of the two attributes representing spatial coordinates.
pattern.summary()<jupyter_output>Point Pattern
1000 points
Bounding rectangle [(116.0976455,39.704590036), (116.69697417,40.162781676)]
Area of window: 0.27460738620632696
Intensity estimate for window: 3641.5626462743703
           wgs84_x    wgs84_y                   mark_0  \
375496  116.271069  39.921151                  Leisure   
110068  116.359717  39.971339              Residential   
249749  116.274537  39.887298       Daily Life Service   
361390  116.301160  39.977463       Daily Life Service   
260324  116.403919  39.888717  Government & Enterprise   

                            mark_1  
375496  POINT (116.27107 39.92115)  
110068  POINT (116.35972 39.97134)  
249749  POINT (116.27454 39.88730)  
361390  POINT (116.30116 39.97746)  
260324  POINT (116.40392 39.88872)<jupyter_text>Mean Nearest Neighbor Distance Statistics The mean distance ($\overline d_{min}$) between points and their nearest neighbour is assumed to be normally distributed for completely spatial random (CSR) distribution. We can test whether the target distribution is the outcome of CSR or not.<jupyter_code>pattern.knn() # the default is the one nearest neighbor<jupyter_output><empty_output><jupyter_text>Outcome: - tuple[0]: the id of the first nearest neighbour of each data point;- tuple[1]: corresponding distances<jupyter_code>pattern.knn(2)   # the first two nearest neighbors and corresponding distances
pattern.nnd   # tuple[1]
pattern.min_nnd
pattern.max_nnd
pattern.mean_nnd
pattern.plot()<jupyter_output><empty_output><jupyter_text>Nearest Neighbor Analysis Functions **$G$-function (event-event)** measures the proportion of nearest neighbor distances that are less than $d$. The target point $point_i$ only include those whose distance to the marginal is larger than $d$.$$G(d)=\frac{\sum_i^N I_i(d)}{N},\ I_i(d)=\left\{              \begin{array}{cl}             1, & if\ \overline d_{min}(point_i)<d; \\             0, & else             \end{array}\right.$$ For CSR pattern, the expectation of $G-$function is:$$E[G(d)]=1-e^{-\lambda\pi d^2}$$ Then we can calculate the max difference of the $G-$function between the target and the CSR distributions:$$D_d = max|E[G(d)]-G(d)|$$<jupyter_code>gf = distance_statistics.G(pattern, intervals=20)
gf.plot()
gf.plot(qq=True)   #quantile-quantile plot
# x values in the plot
gf.d
# y values in the plot
gf.G<jupyter_output><empty_output><jupyter_text>It is obvious that the $G$ grows quicker than the CSR when $d$ is small, and then slowly increases at large $d$. It then comes a reasonable assumption that the points are clustered compared to CSR. Then we need to simulate CSR many times (such as >1000) to come to a final conclusion based on significant tests. We can use **Simulation Envelopes** to achieve such simulations and significance test.<jupyter_code>realizations = PoissonPointProcess(pattern.window, pattern.n, 100, asPP=True) # simulate CSR 100 times
genv = distance_statistics.Genv(pattern, intervals=20, realizations=realizations) # call Genv to generate simulation envelope
genv
genv.observed
genv.plot()<jupyter_output><empty_output><jupyter_text>The $G$ curve is obviously outside and upon the CSR significant curves, so we can conclude that the points are outcomes of cluster point process. **$F$-function (point-event)** is similar to the $G-$function, but the distance is not to the points in target events, but to the points coming from CSR process.<jupyter_code>ff = distance_statistics.F(pattern, intervals=20)
ff.plot()
ff.plot(qq=True)
fenv = distance_statistics.Fenv(pattern, intervals=20, realizations=realizations) # call Genv to generate simulation envelope
fenv
fenv.plot()<jupyter_output><empty_output><jupyter_text>**$J-$function** is based on these two functions. One of its advantage is that it is based on the non-CSR null hypothesis. When $J1$, it is a regular process.$$J(d) = \frac{1-G(d)}{1-F(d)}$$<jupyter_code>jf = distance_statistics.J(pattern, intervals=50)
jf.plot()  
# The warning may due to F(d) becoming 1 in CSR process, but I do not find out how to fix it until now. This also influence the belowed Jenv calculations. 
# If you know, please contact me.
# jenv = distance_statistics.Jenv(pattern, intervals=50, realizations=realizations)
# jenv.plot()<jupyter_output><empty_output><jupyter_text>The $J(d)$ is obviously under $J(d)=1$, so we can conclude that the underlining process has a large propability to be a cluster process. There also has **$S-$function** as $S(d)=\sum_{i=1}^N[F(d_i)-G(d_i)]$ or $S(d)=\sum_{i=1}^N|F(d_i)-G(d_i)|$, which is quite similar to $J-$function to compare the $F(d)$ and $G(d)$. The latter is better under most circumstances because the positive and negative values can not be added to be zero. It seems that the $S-$function is not included in the package.  Second-order point pattern analysis  The $F-$function, $G-$function, and their combinations only consider the nearest distances which rely on a single spatial scale. We can call them `the first-order measurement'. Acturally, we can use all the points including higher order neighbors to consider distances between all target points.  $K-$function is related to the expect number of points within a specific circle.$$K(d)=A\sum_{i=1,i\neq j}^N\sum_{j=1}^Nh_i(d)I_d(i,j)/N^2$$$A$ is the area of target zones; $h_i(d)=1$ if the circle with point $i$ as centre and $d$ as its radius located completely within the study area, otherwise $h_i(d)$ equals to the reciprocal intersected areas; $I_d(i,j)$ is 1 when $d_{ij}\leq d$, else 0; $N$ is the total number of points.If the target points are generated from a CSR process, the expect point number will follow Poisson distribution. The expect number of points within a circle with $d$ as its radius is $n\pit^2/A$.We can also compare target process with CSR via $L-$function:$$L(d)=\sqrt{K(d)/\pi}-d$$<jupyter_code>kf = distance_statistics.K(pattern)
kf.plot()
kenv = distance_statistics.Kenv(pattern, intervals=20, realizations=realizations)
kenv.plot()
lf = distance_statistics.L(pattern)
lf.plot()
lenv = distance_statistics.Lenv(pattern, intervals=20, realizations=realizations)
lenv.plot()<jupyter_output><empty_output><jupyter_text>2. Simulating point process - Random pattern ($N$-conditioned, $\lambda$-conditioned)- Cluster pattern- Regular patternAll refer to https://nbviewer.jupyter.org/github/pysal/pointpats/blob/master/notebooks/process.ipynb.<jupyter_code>from pointpats import PoissonPointProcess, PoissonClusterPointProcess, Window, poly_from_bbox, PointPattern
import libpysal as ps
from libpysal.cg import shapely_ext
import numpy as np
# open the virginia polygon shapefile
va = ps.io.open('../Data/Beijing6Ring_Polygon.shp')
polys = [shp for shp in va]
# Create the exterior polygons for VA from the union of the county shapes
state = shapely_ext.cascaded_union(polys)
# create window from virginia state boundary
window = Window(state.parts)<jupyter_output>SHAPEFILE WARNING: Polygon 3 topology has been fixed. (ccw -> cw)<jupyter_text>Random pattern $N$-conditioned<jupyter_code>np.random.seed(5)
samples = PoissonPointProcess(window, 200, 1, conditioning=False, asPP=False)  
# by specifying "asPP" false, we can generate a point series
# by specifying "asPP" true, or we will directly get a point pattern, samples.realizations[0] will be like <pointpats.pointpattern.PointPattern at 0x7fb5f41af208> instead of an array series. 
samples.realizations[0] # simulated event points
# build a point pattern from the simulated point series
pp_csr = PointPattern(samples.realizations[0])
pp_csr.plot(window=True, hull=True, title='Random Point Pattern')  # hull, window - whether draw the maximum bounding polygon and rectangle
pp_csr.plot(window=False, hull=False, title='Random Point Pattern')<jupyter_output><empty_output><jupyter_text>$\lambda$-conditioned (The function will calculate the intensity $\lambda$ according to your input N(=200, in our example), we will obtain points with similar intensities, but not strictly with similar point numbers.)<jupyter_code># by specifying "conditioning" True, we can simulate a lamda-conditioned CSR
np.random.seed(5)
samples = PoissonPointProcess(window, 200, 1, conditioning=True, asPP=True)
samples
pp_csr = samples.realizations[0] # simulated point pattern
pp_csr
pp_csr.plot(window=True, hull=True, title='Random Point Pattern')
pp_csr.n
samples = PoissonPointProcess(window, 100, 1, conditioning=True, asPP=True)
samples
pp_csr = samples.realizations[0] # simulated point pattern
pp_csr
pp_csr.plot(window=True, hull=True, title='Random Point Pattern')
pp_csr.n<jupyter_output><empty_output><jupyter_text>Cluster pattern There are two sources of clustering:-Contagion: presence of events at one location affects probability of events at another location (**correlated point process**)-Heterogeneity: intensity $\lambda$ varies with location (**heterogeneous Poisson point process**)The PoissonClusterPointProcess seems mainly focus on the first one. First generate parent points using CSR generation, based on both $N-$conditioned and $\lambda-$conditioned by set `conditioning' parameter True or False. Then we generate offspring events independently around the parents within a circle of radius $r$ centered on the parent.<jupyter_code>csamples = PoissonClusterPointProcess(window, 200, 5, 0.05, 1, asPP=True, conditioning=False)
# parameter: total number-200, parents-5, offsprings-200/5, radius-0.05, realization nums-1
csamples
csamples.num_parents 
csamples.children
pp_pcp = csamples.realizations[0]
pp_pcp
pp_pcp.plot(window=True, hull=True, title='Clustered Point Pattern') #plot the first realization<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09716147045137273, 'fraction_numerical': 0.03946021405304793, 'mean_word_length': 5.4078711985688726, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 55, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '173763', 'license': 'MIT', 'n_tokens_mistral': 3568, 'n_tokens_neox': 3221, 'n_words': 1145, 'path': 'Part 3 - Point pattern analysis.ipynb', 'repo_name': 'ClaireXing/Python-for-geospatial-data-analysis_updating', 'size': '10745'}