Error when trying to mask raster to shapefile after interpolation using R - "sp supports Z dimension only for POINT and MULTIPOINT"

Question

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
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

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