-
Notifications
You must be signed in to change notification settings - Fork 8
/
sim.Krig.html
391 lines (339 loc) · 14.8 KB
/
sim.Krig.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml"><head><title>R: Conditional simulation of a spatial process</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<link rel="stylesheet" type="text/css" href="R.css" />
</head><body>
<table width="100%" summary="page for sim.spatialProcess {fields}"><tr><td>sim.spatialProcess {fields}</td><td style="text-align: right;">R Documentation</td></tr></table>
<h2>Conditional simulation of a spatial process</h2>
<h3>Description</h3>
<p>Generates exact (or approximate) random draws from the conditional
distribution of a spatial process given specific observations. This is a
useful way to characterize the uncertainty in the predicted process from
data. This is known as conditional simulation in geostatistics or
generating an ensemble prediction in the geosciences. sim.Krig.grid can
generate a conditional sample for a large regular grid but is restricted
to stationary correlation functions.
</p>
<h3>Usage</h3>
<pre>
sim.spatialProcess(object, xp, M = 1, verbose = FALSE, ...)
sim.Krig(object, xp, M = 1, verbose = FALSE, ...)
sim.Krig.approx(object, grid.list = NULL, M = 1, nx = 40, ny = 40,
verbose = FALSE, extrap = FALSE,...)
sim.mKrig.approx(mKrigObject, predictionPoints = NULL,
predictionPointsList = NULL, simulationGridList =
NULL, gridRefinement = 5, gridExpansion = 1 + 1e-07, M
= 1, nx = 40, ny = 40, nxSimulation = NULL,
nySimulation = NULL, delta = NULL, verbose = FALSE,...)
sim.fastTps.approx(fastTpsObject,
predictionPointsList, simulationGridList =
NULL, gridRefinement = 5, gridExpansion = 1 + 1e-07, M
= 1, delta = NULL, verbose=FALSE,...)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
<tr valign="top"><td><code>delta</code></td>
<td>
<p>If the covariance has compact support the simulation method can
take advantage of this. This is the amount of buffer added for the simulation domain in the circulant embedding method.
A minimum size would be <code>theta</code> for the Wendland but a multiple of this maybe needed to obtain a positive definite
circulant covariance function. </p>
</td></tr>
<tr valign="top"><td><code>extrap</code></td>
<td>
<p> If FALSE conditional process is not evaluated outside
the convex hull of observations. </p>
</td></tr>
<tr valign="top"><td><code>fastTpsObject</code></td>
<td>
<p>The output object returned by fastTps</p>
</td></tr>
<tr valign="top"><td><code>grid.list</code></td>
<td>
<p>Grid information for evaluating the conditional
surface as a grid.list.</p>
</td></tr>
<tr valign="top"><td><code>gridRefinement</code></td>
<td>
<p>Amount to increase the number of grid points
for the simulation grid.</p>
</td></tr>
<tr valign="top"><td><code>gridExpansion</code></td>
<td>
<p>Amount to increase the size of teh simulation
grid. This is used to increase the simulation domain so that the
circulant embedding algorithm works.</p>
</td></tr>
<tr valign="top"><td><code>mKrigObject</code></td>
<td>
<p>An mKrig Object</p>
</td></tr>
<tr valign="top"><td><code>M</code></td>
<td>
<p>Number of draws from conditional distribution.</p>
</td></tr>
<tr valign="top"><td><code>nx</code></td>
<td>
<p> Number of grid points in prediction locations for x coordinate.</p>
</td></tr>
<tr valign="top"><td><code>ny</code></td>
<td>
<p> Number of grid points in prediction locations for x coordinate.</p>
</td></tr>
<tr valign="top"><td><code>nxSimulation</code></td>
<td>
<p> Number of grid points in the circulant embedding simulation x coordinate.</p>
</td></tr>
<tr valign="top"><td><code>nySimulation</code></td>
<td>
<p> Number of grid points in the circulant embedding simulation x coordinate.</p>
</td></tr>
<tr valign="top"><td><code>object</code></td>
<td>
<p>A Krig object.</p>
</td></tr>
<tr valign="top"><td><code>predictionPoints</code></td>
<td>
<p>A matrix of locations defining the
points for evaluating the predictions.</p>
</td></tr>
<tr valign="top"><td><code>predictionPointsList</code></td>
<td>
<p> A <code>grid.list</code> defining the
rectangular grid for evaluating the predictions.</p>
</td></tr>
<tr valign="top"><td><code>simulationGridList</code></td>
<td>
<p> A <code>gridlist</code> describing grid for
simulation. If missing this is created from the range of the
locations, <code>nx</code>, <code>ny</code>, <code>gridRefinement</code>, and <code>gridExpansion</code>
or from the range and and <code>nxSimulation</code>, <code>nySimulation</code>.</p>
</td></tr>
<tr valign="top"><td><code>xp</code></td>
<td>
<p>Same as predictionPoints above.</p>
</td></tr>
</table>
<table summary="R argblock">
<tr valign="top"><td><code>...</code></td>
<td>
<p>Any other arguments to be passed to the predict function.
Usually this is the <code>Z</code> or <code>drop.Z</code> argument when there are
additional covariates in the fixed part of the model.
(See example below.) </p>
</td></tr>
<tr valign="top"><td><code>verbose</code></td>
<td>
<p>If true prints out intermediate information. </p>
</td></tr>
</table>
<h3>Details</h3>
<p>These functions generate samples from a conditional multivariate
distribution, or an approximate one, that describes the uncertainty in the estimated spatial
process under Gaussian assumptions. An important assumption throughout
these functions is that all covariance parameters are fixed at their
estimated or prescribed values from the passed object.
</p>
<p>Given a spatial process h(x)= P(x) + g(x) observed at
</p>
<p>Y.k = Z(x.k)d + P(x.k) + g(x.k) + e.k
</p>
<p>where P(x) is a low order, fixed polynomial and g(x) a Gaussian spatial
process and Z(x.k) is a vector of covariates that are also indexed by space (such as elevation). Z(x.k)d is a linear combination of the the covariates with the parameter vector d being a component of the fixed part of the
model and estimated in the usual way by generalized least squares.
</p>
<p>With Y= Y.1, ..., Y.N,
the goal is to sample the conditional distribution of the process.
</p>
<p>[h(x) | Y ] or the full prediction Z(x)d + h(x)
</p>
<p>For fixed a covariance this is just a multivariate normal sampling
problem. <code>sim.Krig.standard</code> samples this conditional process at
the points <code>xp</code> and is exact for fixed covariance parameters.
<code>sim.Krig.grid</code> also assumes fixed covariance parameters and does
approximate sampling on a grid.
</p>
<p>The outline of the algorithm is
</p>
<p>0) Find the spatial prediction at the unobserved locations based on
the actual data. Call this h.hat(x) and this is the conditional mean.
</p>
<p>1) Generate an unconditional spatial process and from this process
simluate synthetic observations. At this point the approximation is
introduced where the field at the observation locations is
approximated using interpolation from the nearest grid points.
</p>
<p>2) Use the spatial prediction model ( using the true covariance) to
estimate the spatial process at unobserved locations.
</p>
<p>3) Find the difference between the simulated process and its
prediction based on synthetic observations. Call this e(x).
</p>
<p>4) h.hat(x) + e(x) is a draw from [h(x) | Y ].
</p>
<p><code>sim.spatialProcess</code> Follows this algorithm exactly. For the case of an
addtional covariate this of course needs to be included. For a model
with covariates use <code>drop.Z=TRUE</code> for the function to ignore prediction
using the covariate and generate conditional samples for just the spatial
process and any low order polynomial. Finally, it should be noted that this
function will also work with an <code>mKrig</code> object because the essential
prediction information in the mKrig and spatialProcess objects are the same.
The naming is through convenience.
</p>
<p><code>sim.Krig</code> Also follows this algorithm exactly but for the older <code>Krig</code> object. Note the inclusion of
drop.Z=TRUE or FALSE will determine whether the conditional simulation
includes the covariates Z or not. (See example below.)
</p>
<p><code>sim.Krig.approx</code> and <code>sim.mKrig.approx</code> evaluate the
conditional surface on grid and simulates the values of h(x) off the
grid using bilinear interpolation of the four nearest grid
points. Because of this approximation it is important to choose the
grid to be fine relative to the spacing of the observations. The
advantage of this approximation is that one can consider conditional
simulation for large grids – beyond the size possible with exact
methods. Here the method for simulation is circulant embedding and so
is restricted to stationary fields. The circulant embedding method is
known to fail if the domain is small relative to the correlation
range. The argument <code>gridExpansion</code> can be used to increase the
size of the domain to make the algorithm work.
</p>
<p><code>sim.fastTps.approx</code> Is optimized for the approximate thin plate
spline estimator in two dimensions and <code>k=2</code>. For efficiency the
ensemble prediction locations must be on a grid.
</p>
<h3>Value</h3>
<p><code>sim.Krig and sim.spatialProcess</code>
a matrix with rows indexed by the locations in <code>xp</code> and columns being the
<code>M</code> independent draws.
</p>
<p><code>sim.Krig.approx</code> a list with components <code>x</code>, <code>y</code> and <code>z</code>.
x and y define the grid for the simulated field and z is a three dimensional array
with dimensions <code>c(nx, ny, M)</code> where the
first two dimensions index the field and the last dimension indexes the draws.
</p>
<p><code>sim.mKrig.approx</code> a list with <code>predictionPoints</code> being the
locations where the field has been simulated.If these have been created from a
grid list that information is stored in the <code>attributes</code> of <code>predictionPoints</code>.
<code>Ensemble</code> is a
matrix where rows index the simulated values of the field and columns
are the different draws, <code>call</code> is the calling sequence. Not that if <code>predictionPoints</code>
has been omitted in the call or is created beforehand using <code>make.surface.grid</code> it is
easy to reformat the results into an image format for ploting using <code>as.surface</code>.
e.g. if <code>simOut</code> is the output object then to plot the 3rd draw:
</p>
<pre>
imageObject<- as.surface(simOut$PredictionGrid, simOut$Ensemble[,3] )
image.plot( imageObject)
</pre>
<p><code>sim.fastTps.approx</code> is a wrapper function that calls <code>sim.mKrig.approx</code>.
</p>
<h3>Author(s)</h3>
<p>Doug Nychka</p>
<h3>See Also</h3>
<p> sim.rf, Krig, spatialProcess</p>
<h3>Examples</h3>
<pre>
## Not run:
# conditional simulation with covariates
# colorado climate example
data(COmonthlyMet)
fit1E<- spatialProcess(CO.loc,CO.tmin.MAM.climate, Z=CO.elev )
# conditional simulation at missing data
good<- !is.na(CO.tmin.MAM.climate )
infill<- sim.spatialProcess( fit1E, xp=CO.loc[!good,],
Z= CO.elev[!good], M= 10)
# get an elevation grid ... NGRID<- 50 gives a nicer image but takes longer
NGRID <- 25
# get elevations on a grid
COGrid<- list( x=seq( -109.5, -101, ,NGRID), y= seq(39, 41.5,,NGRID) )
COGridPoints<- make.surface.grid( COGrid)
# elevations are a bilinear interpolation from the 4km
# Rocky Mountain elevation fields data set.
data( RMelevation)
COElevGrid<- interp.surface( RMelevation, COGridPoints )
# NOTE call to sim.Krig treats the grid points as just a matrix
# of locations the plot has to "reshape" these into a grid
# to use with image.plot
SEout<- sim.spatialProcess( fit1E, xp=COGridPoints, Z= COElevGrid, M= 30)
# for just the smooth surface in lon/lat
# SEout<- sim.spatialProcess( fit1E, xp=COGridPoints, drop.Z=TRUE, M= 30)
# in practice M should be larger to reduce Monte Carlo error.
surSE<- apply( SEout, 2, sd )
image.plot( as.surface( COGridPoints, surSE))
points( fit1E$x, col="magenta", pch=16)
## End(Not run)
data( ozone2)
set.seed( 399)
# fit to day 16 from Midwest ozone data set.
out<- Krig( ozone2$lon.lat, ozone2$y[16,], Covariance="Matern",
theta=1.0,smoothness=1.0, na.rm=TRUE)
# NOTE theta =1.0 is not the best choice but
# allows the sim.rf circulant embedding algorithm to
# work without increasing the domain.
#six missing data locations
xp<- ozone2$lon.lat[ is.na(ozone2$y[16,]),]
# 5 draws from process at xp given the data
# this is an exact calculation
sim.Krig( out,xp, M=5)-> sim.out
# Compare: stats(sim.out)[3,] to Exact: predictSE( out, xp)
# simulations on a grid
# NOTE this is approximate due to the bilinear interpolation
# for simulating the unconditional random field.
# also more grids points ( nx and ny) should be used
sim.Krig.approx(out,M=5, nx=20,ny=20)-> sim.out
# take a look at the ensemble members.
predictSurface( out, grid= list( x=sim.out$x, y=sim.out$y))-> look
zr<- c( 40, 200)
set.panel( 3,2)
image.plot( look, zlim=zr)
title("mean surface")
for ( k in 1:5){
image( sim.out$x, sim.out$y, sim.out$z[,,k], col=tim.colors(), zlim =zr)
}
## Not run:
data( ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]
O3.fit<- mKrig( x,y, Covariance="Matern", theta=.5,smoothness=1.0, lambda= .01 )
set.seed(122)
O3.sim<- sim.mKrig.approx( O3.fit, nx=100, ny=100, gridRefinement=3, M=5 )
set.panel(3,2)
surface( O3.fit)
for ( k in 1:5){
image.plot( as.surface( O3.sim$predictionPoints, O3.sim$Ensemble[,k]) )
}
# conditional simulation at missing data
xMissing<- ozone2$lon.lat[!good,]
O3.sim2<- sim.mKrig.approx( O3.fit, xMissing, nx=80, ny=80,
gridRefinement=3, M=4 )
## End(Not run)
## Not run:
#An example for fastTps:
data(ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]
O3FitMLE<- fastTpsMLE( x,y, theta=1.5 )
O3Obj<- fastTps( x,y, theta=1.5, lambda=O3FitMLE$lambda.MLE)
# creating a quick grid list based on ranges of locations
grid.list<- fields.x.to.grid( O3Obj$x, nx=100, ny=100)
O3Sim<- sim.fastTps.approx( O3Obj,predictionPointsList=grid.list,M=5)
# controlling the grids
xR<- range( x[,1], na.rm=TRUE)
yR<- range( x[,2], na.rm=TRUE)
simulationGridList<- list( x= seq(xR[1],xR[2],,400), y= seq( yR[1],yR[2], ,400))
# very fine localized prediction grid
O3GridList<- list( x= seq( -90.5,-88.5,,200), y= seq( 38,40,,200))
O3Sim<- sim.fastTps.approx( O3Obj, M=5, predictionPointsList=O3GridList,
simulationGridList = simulationGridList)
# check
plot( O3Obj$x)
US( add=TRUE)
image.plot( as.surface( O3GridList,O3Sim$Ensemble[,1] ), add=TRUE)
points( O3Obj$x, pch=16, col="magenta")
## End(Not run)
</pre>
<hr /><div style="text-align: center;">[Package <em>fields</em> version 9.9 <a href="00Index.html">Index</a>]</div>
</body></html>