I have a SpatialPointsDataFrame with elevation data and 51'257 elements. I also have a RasterLayer object of dimensions c(114,210) covering roughly the same area (all the spatial points are within the raster). I now would like to interpolate the spatial points to a grid of the same size and resolution as my RasterLayer object. In particular, I would like to use the automap::autoKrige() function to do so.
Here is how I proceed:
# generate dummy data:
x.vals <- runif(n = 51257, min = 2603330, max = 2608098)
y.vals <- runif(n = 51257, min = 1135929, max = 1138178)
elevation.vals <- rnorm(n = 51257, mean = 0, sd = 1)
make a SpatialPointsDataFrame from dummy data:
elevation.spdf <- SpatialPointsDataFrame(coords = cbind(x.vals,y.vals),
data = data.frame(elevation.vals),
proj4string = CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
generate a dummy raster:
my.raster <- raster(nrows = 114, ncols = 210, xmn = 383187.5, xmx = 388437.5,
ymn = 5136362, ymx = 5139212, vals = NULL,
crs = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
reproject the SpatialPointsDataFrame to the target raster projection
elevation.spdf.reproj <- spTransform(x = elevation.spdf, CRSobj = CRS(my.raster@crs@projargs))
create a new_data object of type sp to be used in the autoKrige() function:
spatial.points <- SpatialPoints(coords = coordinates(my.raster),
proj4string = CRS(my.raster@crs@projargs))
spatial.grid <- SpatialPixels(points = spatial.points)
now do the interpolation:
automap::autoKrige(formula = elevation.vals ~ x.vals + y.vals,
input_data = elevation.spdf.reproj,
new_data = spatial.grid)
However, when I do that, I get the following error from the autoKrige function:
Error in locs[valid.pattern, ] : (subscript) logical subscript too long
In addition: Warning messages:
1: 'newdata' had 23940 rows but variables found have 51257 rows
2: 'newdata' had 23940 rows but variables found have 51257 rows
It seems that my SpatialPointsDataFrame has too many elements. Alternatively, I tried the approach described here, producing a convex hull around my spatial points to interpolate to. However, I get a similar error message:
automap::autoKrige(formula = bedrock.vals ~ x.vals + y.vals,
input_data = elevation.spdf.reproj,
new_data = create_new_data(obj = elevation.spdf.reproj))
Error in locs[valid.pattern, ] : (subscript) logical subscript too long
In addition: Warning messages:
1: 'newdata' had 5006 rows but variables found have 51257 rows
2: 'newdata' had 5006 rows but variables found have 51257 rows
My intuition is that my SpatialPointsDataFrame is too large (larger than the grid I want to interpolate to) and therefore the Kriging fails. I tried randomly thinning my SpatialPointsDataFrame, but the same problem still occurs when the number of elements is lower than the number of cells in my.raster. Curiously, the autoKrige seems to work when I thin the SpatialPointsDataFrame to exactly 23940 elements, but the computations take forever (I aborted after 36 Minutes).
I am aware of the explanations provided here and here, but they both do not use the autoKrige() function.
Does someone have an idea of what is going on here and how to solve it?
System information: R version 3.6.2 // RStudio version 1.2.5033 // Mac OS 10.14.6
SpatialPixelsfunction instead ofSpatialPixelsDataFrame. – Andre Silva Jun 12 '20 at 17:01