critwire · Learn · Learn's data

Point patterns s2

Prerequisites: Bias in records (s1), Grids (g6)

Are points more bunched together than chance would put them? Four ways to ask, and what each assumes: quadrats, the nearest-neighbour index, DBSCAN and kernel density, on the Tasmanian devil's records. The SQL runs in your own PostgreSQL with PostGIS and the kernel density in the Python environment, on the lab tables (the setup page).

Concept

Are points more bunched together than chance would put them, and what does each way of asking assume?

Worked example

Quadrats. The cells are the equal-area hexagons over Tasmania's main island that g6 drew, at edges of 5 and 10 km. Global Moran's I, s4 builds them as the table s4_cells, one row per cell with its edge (edge_m), its number of devil records (n) and its centroid (cx, cy); its first block, s4-cells, makes them. Every download of this module runs that block first, right after connecting, so s2.sql and s2.py work on their own; to run the blocks one at a time, run s4-cells once first. The variance-to-mean ratio of records per hexagon, at both sizes:

Checked 2026-11-02 on PostgreSQL 16.15, PostGIS 3.4.2, GEOS 3.12.1, GDAL 3.8.4; Python 3.12.3 with requirements.txt sha256 d5bbe28024cc.

s2-quadrats · PostgreSQL with PostGIS 

SELECT edge_m, count(*) AS cells, sum(n) AS records,
       round(avg(n), 3) AS mean,
       round(var_samp(n), 2) AS variance,
       round(var_samp(n) / avg(n), 1) AS vmr
FROM s4_cells
GROUP BY edge_m
ORDER BY edge_m;
-- result: s2-quadrats · 2 rows · sha256 0cefd3b143f3
edge_m | cells | records | mean | variance | vmr
5000 | 1094 | 1972 | 1.803 | 115.68 | 64.2
10000 | 296 | 1990 | 6.723 | 934.23 | 139.0

In s2.sql, with the blocks it runs after.

At 5 km the 1,094 cells hold 1.803 records on average, with a variance of 115.68: the variance is 64.2 times the mean, where randomness gives about 1. At 10 km the ratio is 139.0. var_samp is the sample variance (dividing by n − 1); the ratio grows with the cell size because larger cells gather more of each cluster into one count, as in g6.

Lab

  1. The nearest-neighbour index for the records inside the main island, counting every record and then each distinct location once. <-> with LIMIT 1 finds each point's nearest other point (Spatial indexes, g5); DISTINCT ON keeps one record per pair of coordinates:

    s2-nni · PostgreSQL with PostGIS 

    WITH island AS (
  SELECT (part).geom
  FROM (SELECT ST_Dump(ST_Transform(geom, 3577)) AS part
        FROM admin1_shapes WHERE country_code = 'AU' AND name = 'Tasmania') AS parts
  ORDER BY ST_Area((part).geom) DESC
  LIMIT 1
), devils AS (
  SELECT o.id, ST_Transform(o.geom, 3577) AS geom
  FROM occurrences AS o JOIN species AS s ON s.id = o.species_id, island
  WHERE s.scientific_name = 'Sarcophilus harrisii'
    AND ST_Contains(island.geom, ST_Transform(o.geom, 3577))
), sets AS (
  SELECT 'every record' AS points, id, geom FROM devils
  UNION ALL
  (SELECT DISTINCT ON (ST_X(geom), ST_Y(geom)) 'distinct locations', id, geom
   FROM devils ORDER BY ST_X(geom), ST_Y(geom), id)
), nearest AS (
  SELECT a.points,
         (SELECT ST_Distance(a.geom, b.geom) FROM sets AS b
          WHERE b.points = a.points AND b.id <> a.id
          ORDER BY a.geom <-> b.geom LIMIT 1) AS d_m
  FROM sets AS a
), stats AS (
  SELECT nr.points, count(*) AS n, count(*) FILTER (WHERE nr.d_m = 0) AS zero_m,
         avg(nr.d_m) AS observed_m, ST_Area(island.geom) AS area_m2
  FROM nearest AS nr, island
  GROUP BY nr.points, island.geom
)
SELECT points, n, zero_m,
       round(observed_m::numeric, 1) AS observed_m,
       round((0.5 / sqrt(n / area_m2))::numeric, 1) AS expected_m,
       round((observed_m / (0.5 / sqrt(n / area_m2)))::numeric, 3) AS r,
       round(((observed_m - 0.5 / sqrt(n / area_m2)) / (0.26136 / sqrt(n ^ 2 / area_m2)))::numeric, 1) AS z
FROM stats
ORDER BY n DESC;
    -- result: s2-nni · 2 rows · sha256 0c88009e49e6
points | n | zero_m | observed_m | expected_m | r | z
every record | 1784 | 1109 | 985.8 | 3000.5 | 0.329 | -54.3
distinct locations | 913 | 0 | 2298.4 | 4194.3 | 0.548 | -26.1

    In s2.sql, with the blocks it runs after.

    Counting every record, 1,109 of the 1,784 records on the main island sit on exactly the same spot as another record, so their nearest-neighbour distance is 0. The mean distance is 985.8 m where randomness gives 3,000.5 m: R = 0.329, z −54.3. Counting the 913 distinct locations once, R = 0.548 (z −26.1). Both are far below 1, so the records are clustered either way; the duplicates make it look twice as strong. The 1,784 are fewer than the quadrats' 1,972 because a hexagon that touches land reaches past the outline: 188 records lie in a cell but outside the island, at most 4.2 km from it. A z-score beyond ±1.96 means a mean distance this far from the expected, in either direction, would happen by chance less than 5 % of the time.

  2. DBSCAN with eps 2 km and minpoints 20 (ArcGIS Pro's Density-based Clustering, DBSCAN method). ST_ClusterDBSCAN is a window function: it gives each point its cluster's number, or NULL for noise:

    s2-dbscan · PostgreSQL with PostGIS 

    WITH island AS (
  SELECT (part).geom
  FROM (SELECT ST_Dump(ST_Transform(geom, 3577)) AS part
        FROM admin1_shapes WHERE country_code = 'AU' AND name = 'Tasmania') AS parts
  ORDER BY ST_Area((part).geom) DESC
  LIMIT 1
), devils AS (
  SELECT o.id, ST_Transform(o.geom, 3577) AS geom
  FROM occurrences AS o JOIN species AS s ON s.id = o.species_id, island
  WHERE s.scientific_name = 'Sarcophilus harrisii'
    AND ST_Contains(island.geom, ST_Transform(o.geom, 3577))
), labelled AS (
  SELECT id, geom, ST_ClusterDBSCAN(geom, eps := 2000, minpoints := 20) OVER () AS cluster
  FROM devils
), clusters AS (
  SELECT cluster, count(*) AS points, count(DISTINCT (ST_X(geom), ST_Y(geom))) AS locations
  FROM labelled
  GROUP BY cluster
)
SELECT count(*) FILTER (WHERE cluster IS NOT NULL) AS clusters,
       sum(points) FILTER (WHERE cluster IS NOT NULL) AS clustered,
       sum(points) FILTER (WHERE cluster IS NULL) AS noise,
       max(points) FILTER (WHERE cluster IS NOT NULL) AS largest,
       count(*) FILTER (WHERE cluster IS NOT NULL AND locations = 1) AS one_location_clusters
FROM clusters;
    -- result: s2-dbscan · 1 row · sha256 34fd8cce6ca6
clusters | clustered | noise | largest | one_location_clusters
10 | 871 | 913 | 357 | 1

    In s2.sql, with the blocks it runs after.

    DBSCAN finds 10 clusters holding 871 records and calls 913 noise (the same number as the distinct locations above, by chance); the largest cluster holds 357 records. One "cluster" is a single location repeated 20 times or more: by DBSCAN's rule, 20 records at one spot are 20 points within 2 km of each other. count(DISTINCT (ST_X(geom), ST_Y(geom))) counts a cluster's distinct locations. Change eps or minpoints and the clusters change; neither has a right value, so report them with the result.

  3. Kernel density in Python (scikit-learn), a Gaussian kernel with a 5 km bandwidth, read at the 5 km hexagons' centroids. The coordinates are in metres, so the density comes out per m² and sums to 1 over the plane; multiplying by the number of records and by 10⁶ turns it into records per km². Summed over the cells it should give back roughly the records in them:

    s2-kde · Python

    Runs after: py-connect

    import numpy as np
from sklearn.neighbors import KernelDensity

with connect() as conn:
    points = np.array(conn.execute(
        """SELECT ST_X(g), ST_Y(g)
           FROM (SELECT ST_Transform(o.geom, 3577) AS g
                 FROM occurrences AS o JOIN species AS s ON s.id = o.species_id
                 WHERE s.scientific_name = 'Sarcophilus harrisii') AS devils
           ORDER BY 1, 2"""
    ).fetchall())
    cells = conn.execute(
        "SELECT id, n, cx, cy, ST_Area(geom) FROM s4_cells WHERE edge_m = 5000 ORDER BY id"
    ).fetchall()

kde = KernelDensity(bandwidth=5000, kernel="gaussian").fit(points)
centroids = np.array([(c[2], c[3]) for c in cells])
per_km2 = np.exp(kde.score_samples(centroids)) * len(points) * 1e6
cell_km2 = cells[0][4] / 1e6

print(f"records: {len(points)}  cells: {len(cells)}  cell area: {cell_km2:.2f} km2")
print(f"density x area, summed over the cells: {(per_km2 * cell_km2).sum():.0f} records")
print("densest cells: id | records in cell | records per km2")
for k in np.argsort(-per_km2, kind="stable")[:5]:
    print(f"{cells[k][0]} | {cells[k][1]} | {per_km2[k]:.3f}")
    -- result: s2-kde · python · sha256 130460cdd472
records: 2000  cells: 1094  cell area: 64.95 km2
density x area, summed over the cells: 1895 records
densest cells: id | records in cell | records per km2
917 | 256 | 2.123
947 | 42 | 0.989
948 | 49 | 0.880
888 | 111 | 0.809
916 | 45 | 0.787

    In s2.py and s2.ipynb, with the blocks it runs after.

    The surface peaks at 2.12 records per km² over cell 917, the cell with 256 records (s4's fullest). Cells 947 and 948 come next with far fewer records of their own (42 and 49): they border cell 917 and each other, so the bumps from its records spill into them. Cell 888, with 111 records, ranks only fourth: apart from cell 917, its neighbours hold 6 records between them. Summed over the 1,094 cells of 64.95 km², density times area gives back 1,895 of the 2,000 records; the rest spill into the sea or lie off the main island. score_samples returns the logarithm of the density, hence np.exp.

Result

Both tests say clustered, DBSCAN and the kernel surface show where, and duplicates inflate all four. The quadrats' variance is 64.2 times the mean at 5 km and 139.0 times at 10 km (1 for random). Inside the main island, 1,109 of the 1,784 records sit on exactly the same spot as another record, so their nearest-neighbour distance is 0: R = 0.329 (z −54.3) counting every record, R = 0.548 (z −26.1) counting the 913 distinct locations once; both far below 1. DBSCAN finds 10 clusters holding 871 records and calls 913 noise; the largest holds 357 records, and one "cluster" is a single location repeated 20 times or more. The kernel surface peaks at 2.12 records per km² over cell 917, the cell with 256 records (s4's fullest); summed over the 1,094 cells of 64.95 km² it gives back 1,895 of the 2,000 records (the rest spill into the sea or lie off the main island). Clustered means the records bunch together; s1 says why: they were recorded by a few teams, at places they returned to.

Exercise

The nearest-neighbour index needs an area, and s2-nni used the main island. For the 913 distinct locations, compute R and its z-score with four areas: the main island, the rectangle around the locations along the axes (ST_Envelope), the smallest rectangle around them at any angle (ST_OrientedEnvelope), and their convex hull (ST_ConvexHull, the smallest shape with no inward bends that holds them all). Show each area in km². How much does the area change R, and which area answers the question?

Show answer

s2-ex · PostgreSQL with PostGIS 

WITH island AS (
  SELECT (part).geom
  FROM (SELECT ST_Dump(ST_Transform(geom, 3577)) AS part
        FROM admin1_shapes WHERE country_code = 'AU' AND name = 'Tasmania') AS parts
  ORDER BY ST_Area((part).geom) DESC
  LIMIT 1
), locations AS (
  SELECT DISTINCT ON (ST_X(g), ST_Y(g)) g AS geom
  FROM (SELECT ST_Transform(o.geom, 3577) AS g
        FROM occurrences AS o JOIN species AS s ON s.id = o.species_id, island
        WHERE s.scientific_name = 'Sarcophilus harrisii'
          AND ST_Contains(island.geom, ST_Transform(o.geom, 3577))) AS devils
  ORDER BY ST_X(g), ST_Y(g)
), observed AS (
  SELECT count(*) AS n,
         avg((SELECT ST_Distance(a.geom, b.geom) FROM locations AS b
              WHERE NOT ST_Equals(b.geom, a.geom)
              ORDER BY a.geom <-> b.geom LIMIT 1)) AS observed_m,
         ST_Collect(geom) AS points
  FROM locations AS a
), areas AS (
  SELECT 'main island' AS area, ST_Area(island.geom) AS area_m2 FROM island
  UNION ALL
  SELECT 'bounding rectangle', ST_Area(ST_Envelope(points)) FROM observed
  UNION ALL
  SELECT 'smallest rectangle', ST_Area(ST_OrientedEnvelope(points)) FROM observed
  UNION ALL
  SELECT 'convex hull', ST_Area(ST_ConvexHull(points)) FROM observed
)
SELECT a.area, round((a.area_m2 / 1e6)::numeric) AS km2, o.n,
       round(o.observed_m::numeric, 1) AS observed_m,
       round((0.5 / sqrt(o.n / a.area_m2))::numeric, 1) AS expected_m,
       round((o.observed_m / (0.5 / sqrt(o.n / a.area_m2)))::numeric, 3) AS r,
       round(((o.observed_m - 0.5 / sqrt(o.n / a.area_m2)) / (0.26136 / sqrt(o.n ^ 2 / a.area_m2)))::numeric, 1) AS z
FROM areas AS a, observed AS o
ORDER BY a.area_m2 DESC;
-- result: s2-ex · 4 rows · sha256 ff1aa0fff429
area | km2 | n | observed_m | expected_m | r | z
bounding rectangle | 100881 | 913 | 2298.4 | 5255.8 | 0.437 | -32.5
smallest rectangle | 95274 | 913 | 2298.4 | 5107.7 | 0.450 | -31.8
convex hull | 70537 | 913 | 2298.4 | 4394.9 | 0.523 | -27.6
main island | 64246 | 913 | 2298.4 | 4194.3 | 0.548 | -26.1

In s2.sql, with the blocks it runs after.

Check

The answer: s2-nni's distinct locations, one observed mean distance, and four areas. ST_Collect gathers the locations into one geometry so each shape can be built around them; since the locations are distinct, NOT ST_Equals leaves out only the point itself.

The observed mean is 2,298.4 m in every row, as in s2-nni; only the expected distance changes, with the square root of the area. The main island (64,246 km²) gives R 0.548 and z −26.1, the convex hull (70,537 km²) 0.523 and −27.6, the smallest rectangle (95,274 km²) 0.450 and −31.8, and the rectangle along the axes (100,881 km²) 0.437 and −32.5. Same points, four answers: every area says clustered, but the rectangles add 31,000 to 37,000 km² of sea, where no devil can be recorded, and that empty space alone makes the points look more bunched. The island is the area a record could fall in, so it answers the question best; even then it includes mountains and forest nobody visited, so R below 1 still mixes devils with recording effort. ArcGIS Pro's Average Nearest Neighbor uses a rectangle around all the points when it isn't given an area, so give it the island's.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone