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.shp). Both are in a Google Drive folder and are completely safe to download.
Here is what I have tried:
#Import libraries
library(pacman)
p_load(sf,
spatstat,
maptools,
tidyverse,
ggplot2,
gstat,
sp,
rgdal,
raster,
spdep)
#Set workplace directory to wherever the files were downloaded
setwd("C:\\Users\\path_to_files")
#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)
plot(var.data7)
#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)
plot(wd.data7_krg)
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