I am trying to mask an interpolated raster to a shapefile with the same coordinate reference system using R, but I keep getting an error. The error is shown below. How can I overcome this error and successfully mask this raster to the shapefile?

To reproduce this error, the input csv file is provided here (example_data.csv) and the shapefile used for masking is provided here (field_boundary data). Both are in a Google Drive folder and are completely safe to download.

Here is what I have tried:

#Import libraries

#Set workplace directory to wherever the files were downloaded

#Load data
shp <- st_read('field_boundary.shp')
data <- read.csv('example_data.csv')

#Project field boundary into WGS 84
shp <- st_transform(shp, st_crs('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))

#Generate coordinate matrix for interpolation
lon <- data$longitude
lat <- data$latitude
coords <- cbind(lon, lat)

#Generate spatial DataFrame
data_sp <- SpatialPointsDataFrame(coords, data[, c("data1", "data2", "data3", "data4", "data5", "data6", "data7")], proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))

#Remove any duplicate points
data_sp <- remove.duplicates(data_sp)

#Generate grid to interpolate over
grd <- SpatialPixels(SpatialPoints(makegrid(data_sp, n=50000)), proj4string = proj4string(data_sp))

coordnames(grd) <- c("X", "Y")

plot(grd) #If desired to insepct outcome

#Define formula for kriging
f.data7 <- as.formula(data7 ~ X + Y)
data_sp$X <- coordinates(data_sp[ ,1])
data_sp$Y <- coordinates(data_sp[ ,2])

#Generate variogram
var.data7 <- variogram(f.data7, data = data_sp, cutoff = 0.3)

#Fit and plot variogram
wd.data7 <- fit.variogram(var.data7, vgm('Mat')) #Mattern fit
plot(var.data7, pch = 20, cex = 1.5, col = 'black', ylab = 'Semivariance', xlab = 'Distance (m)', model = wd.data7)

#Ordinary kriging interpolation
wd.data7_krg <- krige(data7~1, data_sp, grd, model = wd.data7)

The kriging interpolation was successful, but now I need to mask it to the shapefile, which is a field boundary. To do this, I transformed the kriged map into a raster.

data7_rast <- raster(wd.data7_krg)
shp <- st_transform(shp, crs(data7_rast)) #Set to same reference system

#Crop and mask
crop_rast <- crop(data7_rast, shp)
data7_mask <- mask(crop_rast, mask = shp)

The following error is thrown here:

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'addAttrToGeom': sp supports Z dimension only for POINT and MULTIPOINT.
use `st_zm(...)` to coerce to XY dimensions

Edit: I have updated the link to provide all of the additional shapefile details, as requested in the comments below.

  • 1
    You need to share the other parts of the shapefile (.dbf, .shx. .prj) which you could put in a zip file. Or save to a geopackage and share that.
    – Spacedman
    Apr 7, 2022 at 17:33
  • And that is why I hate shp…
    – dieghernan
    Apr 7, 2022 at 17:38
  • Have you tried the suggestion in the error to drop the Z and M components with the st_zm function?
    – Spacedman
    Apr 7, 2022 at 17:39
  • @Spacedman Sorry about that! I have updated the hyperlink of the shapefile to include all of the necessary information. To answer your question, yes I have, but I'm not quite sure that I'm doing it right. Apr 8, 2022 at 20:06


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.