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

  • I'm not sure we can help without a sample data set, or code to construct one. Can you put your data somewhere we can get it, or make a sample data set that illustrates the problem and reproduces your errors? – Spacedman Jun 10 '20 at 14:32
  • Thanks for your comment. I edited the question to add code for dummy data, with which the error should be reproducible. I also added my system information. – Abdirizak Jun 12 '20 at 11:43
  • 1
    Did not go over your entire Q/code, but take a look a this tutorial: gis.stackexchange.com/questions/158021/…. Maybe it is something related to your prediction grid (and not the observed points). I saw you used SpatialPixels function instead of SpatialPixelsDataFrame. – Andre Silva Jun 12 '20 at 17:01
  • 1
    Thanks for the hint! This solved the issue. – Abdirizak Jun 15 '20 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.