Spatial Values- Spatial Estimation

Eduardo Serna, Wolfram Research

Madrid, 14 Oct 2021

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

ABSTRACT

Spatial data sets consisting of a set of measured values at specific locations are becoming increasingly important. Examples include temperature, elevation, concentration of minerals, etc.

We will preview upcoming Wolfram Language functionality to perform estimation of missing values in a region of interest, modelling the underlying quantity. Traditionally, this procedure has been known as kriging.

First Example

Let’s look at a spatial data set  of a continuous quantity with large gaps, ozone readings over contiguous United States:

locs=;
ozone={};

plot1=PointValuePlot[locs->ozone,ColorFunction->"Rainbow"]

SpatialEstimate is a very powerful estimator that can combine both large scale trends and local variation.

Let’s create a SpatialEstimatorFunction with a large scale quadratic trend:

sf=SpatialEstimate[locs->ozone,SpatialTrendFunction->2(*Order of polynomial fitted*)]

Evaluate spatial predictor function at uniform random locations in the observation region:

reg = RipleyRassonRegion[locs];
pts=RandomGeoPosition[reg,3000];

plot2 =PointValuePlot[pts->sf[pts],ColorFunction->"Rainbow"];

Compare the two plots:

Row@{Legended[plot1,Placed["Data", Above]], Legended[plot2,Placed[SpatialEstimate, Above]]}

Find estimated value for a particular location:

sf[(Entity["Boston", Entity["City",{"Boston","Massachusetts","UnitedStates"}], "Entity[\"City\", {\"Boston\", \"Massachusetts\", \"UnitedStates\"}]", "city"])["Position"]]

The model-warmup

Let’s warm up with the familiar model of Ordinary Least Squares:

(Subscript[y, i])=f((Subscript[x, i]))+(Subscript[δ, i])

Where δ is an i.i.d zero mean Gaussian random variable  and f is the function that we are trying to parametrically estimate.

Once f is estimated (OverHat[f]), we can infer the value of (OverHat[y]) at a location (Subscript[x, 0]) as:

(Subscript[(OverHat[y]), 0])=(OverHat[f])((Subscript[x, 0]))

Let’s look at an example:

locs=Range[0,1,0.05];
vals=;

dataPlot=ListPlot[Transpose[{locs,vals}],PlotStyle->Black]

fm=LinearModelFit[Transpose[{locs,vals}],{1,x,x^2,x^3},x]

Show[dataPlot,Plot[fm["Function"][x],{x,0,1}]]

The model-SpatialEstimate

The model of Spatial Estimate is the same as (Generalized-more on this later) Least Squares

(Subscript[y, i])=f((Subscript[x, i]))+(Subscript[δ, i])

However, instead of using (Subscript[(OverHat[y]), 0])=(OverHat[f])((Subscript[x, 0])) as an estimator it will use a linear combination of the (Subscript[y, i])

(Subscript[(OverHat[y]), 0])=∑(Subscript[y, i])(Subscript[λ, i])((Subscript[x, 0]))

The (Subscript[λ, i])((Subscript[x, 0])) depend on the point of interest and are chosen to minimise  E[((((Subscript[y, 0])-(Subscript[(OverHat[y]), 0])))^2)]. It is a unbiased minimum variance linear estimator.

Let’s look at the previous example.

sef=SpatialEstimate[locs->vals]

Show[dataPlot,
Plot[{fm["Function"][x], sef[x]},{x,0,1},PlotLegends->{FittedModel, SpatialEstimatorFunction}]
]

The model-SpatialNoiseLevel

Let’s look again at the components of our model

(Subscript[y, i])=f((Subscript[x, i]))+(Subscript[δ, i])

The zero mean random component (Subscript[δ, i]),  is split up into a correlated component (Subscript[η, i]) and a “white noise” i.i.d component (Subscript[ε, i]).

(Subscript[y, i])=f((Subscript[x, i]))+(Subscript[η, i])+(Subscript[ε, i])

The variance of (Subscript[ε, i]) is specified by SpatialNoiseLevel. It is a measure of how noisy the data is or, in other words, how much our estimator should smooth things out (if data is noisier we are less sure of our estimator being close to the data).

est0=SpatialEstimate[locs->vals,SpatialNoiseLevel->0];

estSmall=SpatialEstimate[locs->vals,SpatialNoiseLevel->0.0005];

estLarge=SpatialEstimate[locs->vals,SpatialNoiseLevel->0.005];

Show[dataPlot,Plot[{est0[x],estSmall[x],estLarge[x]},{x,0,1},PlotLegends->{"zero noise level","smaller noise level","larger noise level"}]]

Only with SpatialNoiseLevel -> 0 is SpatialEstimate an exact interpolator.

The model-SpatialTrendFunction

The trend f((Subscript[x, i])) is estimated according to (Generalized) Least Squares, as SpatialNoiseLevel -> ∞, (Subscript[(OverHat[y]), 0])->(OverHat[f])((Subscript[x, 0]))

estCube=SpatialEstimate[locs->vals,SpatialNoiseLevel->10,SpatialTrendFunction->3];

estSquare=SpatialEstimate[locs->vals,SpatialNoiseLevel->10,SpatialTrendFunction->2];

Show[dataPlot,
Plot[{estCube[x],estSquare[x],fm["Function"][x], estCube[x]-fm["Function"][x]},{x,0,1},PlotLegends->{"Cubic Trend", "Square Trend", FittedModel, "Cubic Trend" -FittedModel }]
]

The model-VariogramFunction

The remaining part of our model (Subscript[η, i]) is specified by the Variogram γ(r), which for stationary random processes is a function of the distance.

E[(((η((Subscript[x, i]))-η((Subscript[x, j]))))^2)]=(Subscript[γ, ij])=γ(r)

Let’s see how the Variogram affects the estimate by looking at an example:

data=Quiet@ResourceData["Sample Data: Swiss Rainfall"]

Show[Region[data["ObservationRegion"]],PointValuePlot[data,ColorFunction->"Rainbow"]]

We can specify the values by specifying the Annotation key:

data[{"Annotations","Keys"}]

est=SpatialEstimate[data->"rainfall"]

SpatialEstimate will automatically estimate the Variogram:

vm=est[VariogramFunction]

vm["Visualization"]

The data points above are the Binned (empirical) Variogram.

ListPlot@BinnedVariogramList[data->"rainfall"]

By default a suitable model is chosen, but we can specify a VariogramModel to fit

vm1=EstimatedVariogramModel[data->"rainfall","Linear"]

The fit is worse:

vm1["Visualization"]

{vm["WeightedSumOfResiduals"],vm1["WeightedSumOfResiduals"]}

est1=SpatialEstimate[data->"rainfall",VariogramFunction->vm1]

The estimator for the better fit captures more detail:

pts=RandomPoint[data["ObservationRegion"],3000];
{
Legended[Show[Region[data["ObservationRegion"]],PointValuePlot[pts->est[pts],PlotStyle->PointSize[0.015],ColorFunction->"Rainbow"]],
Placed["Cubic",Above]],
Legended[
Show[Region[data["ObservationRegion"]],PointValuePlot[pts->est1[pts],PlotStyle->PointSize[0.015],ColorFunction->"Rainbow"]],
Placed["Linear",Above]]
}

Models must have a property that amounts to the covariance of points being positive definite.

Estimation of the error

Having a model allows us to compute an estimate of the error (StandardDeviation).

est[RandomPoint[data["ObservationRegion"],2],Around]

est["Visualization", "Bands" -> True,Method-> "3D"]

Large datasets

Let’s look at a relatively large data set.

locs=GeoPosition[]

vals=

GeoDensityPlot[locs->vals, ColorFunction->"Rainbow"]

SpatialEstimate solves linear systems of the order of size the data-set, an operation which is 𝒪((n^3)) .

AbsoluteTiming[SpatialEstimate[locs->vals, Method -> {"MaxPoints"-> All}]]

For large data sets, by default we will solve a local system,  so the cost becomes 𝒪(n) on the size of the data.

AbsoluteTiming[sef=SpatialEstimate[locs->vals]]

sef["Visualization"]

Conclusion

This functionality will be in the upcoming release, and we hope it will be useful to you guys. Take it for a spin and let us now what features you are interested in.

Initialization

SLDebugEchoSwitch[False];