# Spatial prediction on surface: very fine grid vs coarse grid + quick interpolation

Once I have fitted a spatial model (point-referenced data), I need to make a prediction map.

1. A natural approach is to make prediction over a fine grid over the region. However, the required resolution is very high, resulting in >10,000,000 points where I need to make a prediction.
2. Alternatively, I could make prediction on a coarser grid, so that only, say 10,000 predictions are required, and then interpolate it using some 'quick and easy' method (linear interpolation, spline etc.)

The first approach is theoretically optimal, assuming the model gives some sort of BLUP. The second approach might not be as theoretically sound but is computationally feasible (since there is typically no need to deal with large matrices).

Question:

1. Am I losing a lot by adopting approach 2 instead of 1?
2. At what point is it ok to switch from approach 1 to approach 2?
3. In approach 2, how would you determine the resolution of the coarse grid? (apart from whatever computational time you can afford)

(I understand all these depend on context, but general guidelines help)

p.s. I use R, so I wouldn't mind a bit of R in your answer, if you need to.

1. that depends on how many points you have, how they are distributed over space, and how their values vary over space. If you have few points very far apart, than you loose little. If you have many points, but they follow a very smooth pattern, you also loose a little. If you have many points with rough or irregular variations (here smooth, there rough), or few points that show large differences at short separation distances, you may loose too much.
2. at the point where the difference between the two is within acceptable limits. Only you can tell what is acceptable. You might be able to find this out by experimenting with less than 10M points.
3. See 1: this really depends on the characteristics of your point data set.

The following R script (since you asked):

n = 1e3
p = 1e7
x = data.frame(x = runif(n), y = runif(n), z = runif(n))
library(sp)
coordinates(x) = ~x+y
pred = makegrid(x, p)
coordinates(pred)= ~x1+x2
library(gstat)
k = idw(z~1, x, pred)


actually does interpolate from 1000 points to 10,000,000 new points, and runs in 90 seconds using less than 2 Gb of memory.