Spatial Values- Spatial Prediction

Eduardo Serna, Wolfram Research

Madrid, 19 Oct 2022

(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 existing 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:

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

locs=;
ozone={};

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

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 follows:

(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 the 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,and as SpatialNoiseLevel -> ∞, SpatialEstimate approaches (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 }]
]

We can also query the estimator for the trend.

estSquare["SymbolicTrend"]

Plot[{estSquare["GLSTrend"][x], estCube["GLSTrend"][x]}, {x, 0, 1}]

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 the points being positive definite.

SpatialNoiseLevel, revisited

Initialization

BlockRandom[
m=10;
data={locs,vals}=RandomReal[1,{2,m}];
target = RandomReal[1, 100];
,
RandomSeeding->"joiuviuañjabv uha"
]

Example

Let’s look at the Gaussian Variogram, the flat slope at the origin causes it to suffer from conditioning problems if data points are too close together.

gaussVar0 = VariogramModel["Gaussian", {1, 10}, SpatialNoiseLevel -> 0];
gaussVarAuto = VariogramModel["Gaussian", {1, 10}];

gaussVarAuto["Visualization"]

gaussSef=SpatialEstimate[locs->vals, VariogramFunction->gaussVar0]
gaussSefAuto=SpatialEstimate[locs->vals, VariogramFunction->gaussVarAuto]

data2 = SetPrecision[data, 16];

gaussSefArb0=SpatialEstimate[Rule@@data2, VariogramFunction->gaussVar0]
gaussSefArbAuto=SpatialEstimate[Rule@@data2, VariogramFunction->gaussVarAuto]

Regularised (Automatic) is smooth and reasonable, un-regularised fixed precision is nonsense, arbitrary precision is smooth but overfitted.

Through[{gaussSef, gaussSefAuto, gaussSefArb0, gaussSefArbAuto}["Visualization"]]

Close up

{gaussSefArb0["Visualization"], gaussSefArb0["Visualization", PlotRange->{0, 2}]}

An Example-Stars

Initialization

stars=;
angles=;
mags=;

Example

stars=EntityList@EntityClass["Star","ApparentMagnitude"->TakeSmallest[100]]

(*First two are angles in spherical coordinates, 3rd one is distance from the coordinate origin*)
angles=Table[AstroPosition[star, "ICRS"]["Data"], {star, stars}][[All,;;2]]

(*Use ApparentMagnitude to capture the brightness from Earth, which is proportional to Exp[-ApparentMagnitude]*)

mags=Table[star["ApparentMagnitude"], {star, stars}]

ListPlot[Exp[-mags]]

(*Geodesic in the unit sphere*)
SpherDist[{ϕ1_, λ1_}, {ϕ2_, λ2_}]=
	ArcCos[Cos[λ1] Cos[λ2] Cos[ϕ1-ϕ2]+Sin[λ1] Sin[λ2]];

sef=SpatialEstimate[angles->Exp[-mags], DistanceFunction->SpherDist, SpatialTrendFunction->2]

target=Flatten[Table[{ϕ, λ}, {ϕ,  0., 2Pi, Pi/20}, {λ,-Pi/2, Pi/2., Pi/20}], 1];
interpolated=sef[target];
triplets=ArrayFlatten[{{target, Partition[interpolated, 1]}}];
viz=
Show[{
ListDensityPlot[triplets, PlotRange->{{0, 2 Pi}, {-Pi/2, Pi/2}, All}, ColorFunction->"SunsetColors", Frame->False],
PointValuePlot[angles->Exp[-mags], ColorFunction->"SunsetColors"]
}];

ParametricPlot3D[{Cos[ϕ]*Cos[λ],Sin[ϕ]*Cos[λ],Sin[λ]},{ϕ,0,2 π},{λ,-π/2,π/2},Mesh->None,PlotStyle->Texture[viz],Lighting->"Neutral",Boxed->False,Axes->False,SphericalRegion->True,ViewAngle->Pi/10]

Conclusion

This functionality has been in the system since 13.0, and we hope it has been useful to you guys. We have been making some quality of life improvements. Keep using it and let us now what features you are interested in.

Initialization

SLDebugEchoSwitch[False];