0

I am trying to mask a raster to a shapefile boundary, but I am getting an error. How can I correctly perform this mask?

The raw data can be found here, entitled "data_for_question.txt." It is formatted so that users can copy and paste (from the web app) the text directly into an R window and generate a data frame. Otherwise, if one doesn't want to generate the data, the output raster (example_raster.tif) and shapefile (field_boundary.shp) can both also be found in the same link.

Here is what I have tried:

#Import necessary libraries
library(pacman)
p_load(sf,
       spatstat,
       maptools,
       tidyverse,
       ggplot2, 
       gstat,
       sp,
       rgdal,
       raster,
       spdep)

#Read shapefile
shp <- st_read("field_boundary.shp")

#Generate data to run interpolation on and project it to the desired CRS
data_sp <- SpatialPointsDataFrame(coords, 
                                      data[, c("OM", "data2")], 
                                      proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))


#Perform an IDW interpolation:
grd <- SpatialPixels(SpatialPoints(makegrid(data_sp, n=10000)), proj4string = proj4string(data_sp)) #Generate grid for interpolation

plot(grd)

interp <- idw(formula = OM ~ 1, data_sp, grd, idp = 0.5, nmax = 12)
plot(interp) #Makes for a very pretty picture!

#Convert to raster
rast <- raster(interp)
plot(rast)
shp <- st_transform(shp, crs(rast))

#Crop and mask the raster
crop_rast <- crop(rast, shp)
crop_om <- mask(crop_rast, mask = shp)

The error occurs 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
s_o_c_account
  • 111
  • 8
  • Not seeing where `crop_om` is generated before it is used in `mask(crop_om`. – Chris Mar 30 '22 at 18:15
  • Thanks @Chris. I have updated the code to correct that error. The same error is produced. – s_o_c_account Mar 30 '22 at 18:27
  • Same error here. Without bothering with the crop, step, and relying on the error to guide, and guessing one or the other entry would be wrapped in `st_zm`, being more like to be something spatial than raster, – Chris Mar 30 '22 at 19:42
  • and with some renaming...`mask_interp_rast <- mask(interp_rast, mask = sf::st_zm(field_shp))`, where `interp_rast` is your `rast`, `field_shp` your `shp`, gives the masked raster you're looking for, I think. As to why the error, goes to methods and calls therein when signatures are raster and sf. I couldn't directly track down the error message in code (usually good practice)[sp](https://github.com/edzer/sp/tree/main/R), and welcome to Stackoverflow. – Chris Mar 30 '22 at 20:16

0 Answers0