Spatial Statistics Data

Gosia Konwerska, Wolfram Research

Champaign, 7 Oct 2020

(RowDefault[ , (Spacer1), , (Spacer1), ])   #WolframTechConf

ABSTRACT

Spatial point data, also known as spatial point patterns, refers to collections of points (or events) in space. Examples include trees in a forest, gold deposits, positions of stars, earthquakes, crime locations, animal sightings, etc.

Spatial data analysis, as a statistical exploration of point patterns, aims to answer questions about spatial randomness, point intensity and interpoint interactions.

In this talk, we will explore the existing and upcoming Wolfram Language functionality designed for the analysis and modeling of spatial point data.

Spatial Point Data

SpatialPointData is a container for points and observation region. There are three main ways to obtain SpatialPointData:

create SpatialPointData from points, or from points and region, or from points with annotations with/without region

generate SpatialPointData from a point process

obtain SpatialPointData from curated collection of ResourceObjects

Create SpatialPointData from Your Own Data

Fetch data:

cities=CountryData["Canada","LargestCities"]

Input entities directly:

spd=SpatialPointData[cities,reg=Entity["Country","Canada"]]

GeoGraphics[Point[spd["Points"]]]

Available properties:

spd["Properties"]

Quick information about the data:

spd["Summary"]

Extract region information:

spd["ObservationRegion"]

Locations are stored in a compact way:

spd["Points"]

Now this data has no information beyond the locations:

spd["Annotations"]

We can store information about points within SpatialPointData:

names=Map[#["Name"]&,cities];

spd1=SpatialPointData[cities->{"Names"->names},reg]

Now the data has annotations with key “Names”:

spd1["Summary"]

We can use annotations in visualization:

locs=MapThread[Tooltip[Point[#1],Style[#2,14]]&,{Thread[spd1["Points"]],spd1["Annotations"]["Names"]}];

GeoGraphics[{Polygon[spd1["ObservationRegion"]],locs},GeoRange->spd1["ObservationRegion"]]

(side note) Observation Region

Not every region can be a valid observation region - it has to be parameter free, full dimensional, and bounded. This can be tested with SpatialObservationRegionQ:

reg=Entity["Country","Canada"];

SpatialObservationRegionQ[reg]

If no region is provided, the observation region is inferred from the points via Rippley-Rasson estimator RippleyRassonRegion (dilated convex hull):

spd2=SpatialPointData[cities]

The computed observation region:

spd2["ObservationRegion"]

Visualize the points and the observation region:

GeoGraphics[{spd2["ObservationRegion"],Point[spd2["Points"]]}]

Compute RipleyRassonRegion for random points in a disk:

pts=RandomPoint[Disk[],50];
RRreg=RipleyRassonRegion[pts];

Visualize points with computed and original observation regions:

Show[RegionPlot[{Disk[],RRreg},PlotLegends->{"original region","Ripley-Rasson region"}],ListPlot[pts,PlotStyle->Black]]

SpatialPointData from Simulation

Generate a random point collection from a point process using RandomPointConfiguration:

proc=InhomogeneousPoissonPointProcess[Function[{x,y},Exp[10x^2+5y]],2];

sample=RandomPointConfiguration[proc,Disk[]]

Visualize the points within the region:

Show[RegionPlot[sample["ObservationRegion"]],ListPlot[sample]]

SpatialPointData from a ResourceObject

Another way to obtain spatial point data will be to get a curated spatial point data set, which will be available as a ResourceObject:

ro=ResourceObject[];

Extract the SpatialPointData from the ResourceObject:

data=ResourceData[ro]

Quick information about the data:

data["Summary"]

Plot just the points:

ListPlot[data]

Extract annotations:

data["Annotations"]

Create list of types for Species:

ants=data["Annotations"]["Species"]//DeleteDuplicates

Use SpatialPointSelect to create two subsets of original data, one for each species type:

cataglyphis=SpatialPointSelect[data,#Species=="Cataglyphis"&]

messor=SpatialPointSelect[data,#Species=="Messor"&]

Sanity check:

cataglyphis["Annotations"]

Plot both sets:

ListPlot[{cataglyphis["Points"],messor["Points"]},PlotLegends->Sort[ants]]

Point Density

Point density gives a function that describes how the number of points n varies per length, area and volume in the observation region ℛ.

The integral of the desnity function over the region is the total number of points n=(∫_(x∈ℛ))μ(x)dx.

In WL there are three functions dedicated to compute point density using three distinct approaches:

PointDensity
HistogramPointDensity
SmoothPointDensity

PointDensity

PointDensity is an adaptive density estimator (one point per bin):

(*recall ants data*)
data

Compute point density using Voronoi tessellations:

pd=PointDensity[data,"Voronoi"]

Visualize density with points:

Show[pd["DensityVisualization"],ListPlot[data,PlotStyle->Black]]

PointDensityFunction is a function of location and we can plot the density surface:

Plot3D[pd[{x,y}],{x,y}∈data["ObservationRegion"],ColorFunction->"DarkRainbow",MeshFunctions->(#3&)]

HistogramPointDensity

HistogramPointDensity is a partition based estimator where the bin specification controls the smoothing:

hd=HistogramPointDensity[data]

Show[hd["DensityVisualization"],ListPlot[data,PlotStyle->Black]]

The size and shape of the bins can be specified:

hd6=HistogramPointDensity[data,"Hexagon"]

Specify the number of bins:

hd6s=Table[HistogramPointDensity[data,{"Hexagon","Count"->d}],{d,{10,30,60}}];

MapThread[Show[#["DensityVisualization"],ListPlot[data,PlotStyle->Black],PlotLabel->Row[{"bin count = ", #2}]]&,{hd6s,{10,30,60}}]

SmoothPointDensity

SmoothPointDensity is a kernel based estimator of the point density, where bandwidth and kernel functions are used to control the smoothing:

sd=SmoothPointDensity[data]

The SmoothPointDensity specific properties:

sd["Properties"]

The distribution here is the nonparametric DataDistribution, which is either KernelMixtureDistribution or SmoothKernelDistribution:

sd["Distribution"]

This is well defined distribution in Wolfram Language; its PDF an be visualized over the observation region:

Plot3D[PDF[sd["Distribution"],{x,y}],{x,y}∈data["ObservationRegion"]]

Complete Spatial Randomness

Complete spatial randomness means that collection of points in question is independently (from each other) and uniformly distributed over the observation region.

This means that the point density is constant across the observation region and the numbers of points in any subregion are an iid random variables.

The model for complete spatial randomness is PoissonPointProcess:

proc=PoissonPointProcess[30,2];

reg=Disk[];

sample=RandomPointConfiguration[proc,reg]

Visualize the points:

ListPlot[sample,AspectRatio->1]

The point density for PoissonPointProcess is constant everywhere:

PointDensity[PoissonPointProcess[μ,2],{x,y}]

Compare with the sample point density:

PointDensity[sample,"Voronoi"]["DensityVisualization"]

The density visualization does not give a computable answer - therefore we introduced SpatialRandomnessTest:

SpatialRandomnessTest[sample](*returns p-value*)

This test behaves as all other hypothesis tests in WL:

SpatialRandomnessTest[sample,"TestConclusion"]

SpatialRandomnessTest is a super test  - the list of available tests:

SpatialRandomnessTest[sample,"AllTests"]

Use different test:

SpatialRandomnessTest[sample,{"TestDataTable","BesagL"}]

SpatialRandomnessTest for Clustered Data (aka why things are difficult here)

Simulate a point collection in a geographical region from two cluster processes:

meanparentdensity=.006;
clustermean=.3; clustermean1=1;
clusterradius=.1;
dimension=2;

mproc=MaternPointProcess[meanparentdensity,clustermean,clusterradius,dimension];

mproc1=MaternPointProcess[meanparentdensity,clustermean1,clusterradius,dimension];

Define the region for simulation and simulate:

georeg=(Entity["Poland", Entity["Country","Poland"], "Entity[\"Country\", \"Poland\"]", "country"]);

SeedRandom[1];
msample=RandomPointConfiguration[mproc,georeg]

SeedRandom[1];
msample1=RandomPointConfiguration[mproc1,georeg]

Visualize the samples:

{GeoListPlot[msample["Points"]],GeoListPlot[msample1["Points"]]}

Test both samples for spatial randomness:

SpatialRandomnessTest[msample,"TestConclusion"]

SpatialRandomnessTest[msample1,"TestConclusion"]

Conclusion: cluster data is not always evident and can be interpreted as uniform.

Homogeneity Measures

Homogeneity measures are second order measures of how uniform the points are distributed within the observation window. These are the functions available in WL (that can be grouped as follows):

RipleyK, BesagL (*BesagL is a scaled RipleyK*)

EmptySpaceF, NearestNeighborG, SpatialJ (*SpatialJ is a function of EmptySpaceF and NearestNeighborG*)

PairCorrelationG

RipleyK and BesagL

The product of RipleyK[r] and the mean point intensity gives the expected number of points within distance r of a typical point not counting the point itself.

BesagL[r] is defined as (((RipleyK [r](Subscript[b, d])))^(1/d)), where d is the spatial dimension and (Subscript[b, d]) is the volume of a unit ball in (ℝ^d).

Usually RipleyK and BesagL are used to compare the data with the ideal model of PoissonPointProcess:

ideal=PoissonPointProcess[μ,2];

{RipleyK[ideal,r],BesagL[ideal,r]}

Compute RipleyK for a cluster process and compare with the Poisson model:

sample=RandomPointConfiguration[MaternPointProcess[10,30,.1,2],Disk[]];

Any homogeneity measure without explicit radius input will return a PointStatisticFunction object which is a precomputed function to be used as univariate function of the radius:

psf=RipleyK[sample]

Plot to compare with the ideal model:

DiscretePlot[{psf[r],RipleyK[PoissonPointProcess[10,2],r]},{r,0,psf["MaxRadius"],.01},PlotLegends->{"cluster sample","Poisson PP"}]

PointStatisticsFunction stores the following information:

psf["Properties"]

The boundary correction is used to account for the points on the boundary and close to it:

psf["BoundaryCorrection"]

The method can be set using option SpatialBoundaryCorrection.

EmptySpaceF, NearestNeighborG, and SpatialJ

EmptySpaceF[r] gives the probability of finding a point within distance r of an arbitrary location which is typically not a point of the collection.

NearestNeighborG[r] gives the probability of finding another point within distance r of a typical point in the collection.

SpatialJ[r] is the ratio  (1-NearestNeighborG[r])/(1-EmptySpaceF[r]).

If the point data/process exhibit the complete spatial randomness EmptySpaceF and NearestNeighborG are the same:

EmptySpaceF[PoissonPointProcess[μ,d],r]

NearestNeighborG[PoissonPointProcess[μ,d],r]

%%-%

In this situation SpatialJ is a constant function with value 1:

SpatialJ[PoissonPointProcess[μ,d],r]

EmptySpaceF & NearestNeighborG Application Example

Consider a forest data (region scaled: 1 corresponds to (Quantity[924, "ft", feet, "Feet"])):

forestData=SpatialPointData[];

scale= (Quantity[924, "ft", feet, "Feet"]);

ListPlot[forestData,AspectRatio->1]

Compute NearestNeighborG:

nng=NearestNeighborG[forestData]

Compute EmptySpaceF:

esf=EmptySpaceF[forestData]

Compare probability of finding a tree with radius r feet from another tree and from a random location:

DiscretePlot[{nng[Quantity[r,"Feet"]/scale],esf[Quantity[r,"Feet"]/scale]},{r,1,10,.5},AxesLabel->{"r in feet","Probability"},PlotLegends->{"from a data point","from a location"}]

Conclusion: for small radius it is more likely to find another tree from a tree then from a location. So this may hint that data is clustered.

PairCorrelationG

The pair correlation function is also called radial distribution function and it is proportional to the probability density of finding two points at a given distance apart.

The pair correlation function can be understood as a derivative of NearestNeighborG.

PairCorrelationG is not a correlation in a usual sense with values varying between -1 and 1. The range is from 0 to ∞, with 1 corresponding to complete spatial randomness:

PairCorrelationG[PoissonPointProcess[μ,d],r]

PairCorrelationG can indicate hard core radius in a point collection:

sample=RandomPointConfiguration[HardcorePointProcess[100,.1,2],Disk[]]

pcg=PairCorrelationG[sample]

DiscretePlot[pcg[r],{r,0,.6,.01}]

Spatial Binned Point Data

SpatialBinnedPointData is an object to store aggregated spatial point data.

An exercise in data cleanup:

raw=Import["https://dph.illinois.gov/sitefiles/COVIDHistoricalTestResults.json", "Data"];

raw

raw[[1,2]]

raw[[2,2]][[1,2]]

counties=raw[[2,2]][[1,2]][[2;;,1,2]]

The list of counties have to be treated so it is understood by WL:

pos1=Position[raw[[2,2]][[1,2]][[2;;,1,2]],"Out Of State"]

pos2=Position[raw[[2,2]][[1,2]][[2;;,1,2]],"Unassigned"]

Now create bins:

bins=Map[If[#==="Chicago", Entity["City",{"Chicago","Illinois","UnitedStates"}],Entity["AdministrativeDivision",{#~~"County","Illinois","UnitedStates"}]]&,Delete[StringReplace[counties,{"De Witt"->"DeWitt","Jo Daviess"->"JoDaviess","Rock Island"->"RockIsland","St. Clair"->"StClair"}],Join[pos1,pos2]]];

Check that bins looks as expected:

bins[[3]]

Now we need to drop the values corresponding to “Out of State” and “Unassigned”:

confirmed=Delete[raw[[2,2]][[1,2]][[2;;,2,2]],Join[pos1,pos2]]

Check we have the same lengths:

Length/@{confirmed,bins}

Finally create a clean and compact SpatialBinnedPointData:

data=SpatialBinnedPointData[bins->confirmed]

Visualize the data:

{GeoRegionValuePlot[Thread[bins->Log@confirmed],PlotLabel->"log of confirmed cases",PlotLegends->Automatic,ImageSize->250],Plot[Exp[x],{x,0,10},PlotLabel->"Exp[x]",GridLines->All,ImageSize->250]}

Thank You!

Questions?