Search code examples
rrasterr-sfr-raster

How to replace values that correspond to polygons by NA?


I have a raster and a shapefile:

library(cartography)
library(sf)
library(raster)

r <- raster(matrix(rnorm(10*12), nrow=10), xmn = -180, xmx= 180, ymn = -90, ymx= 90)
mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), quiet = TRUE)

I would like to intersect the raster r with the shapefile mtq and make the corresponding pixels to the all polygons as NA (replace the values of the pixels in the raster by NA) and return the raster.


Solution

  • You are likely looking for mask; it lives in both oldish {raster} and shiny new {terra}.

    Note that I had to rewrite your r object a bit, as it was not quite compatible with the Martinique vector object from {cartography}.

    Edit: if, as seems to be indicated in the comments, you are looking for replacing with NAs the values inside the polygon (and not outside) my answer is still raster::mask(), only with a little tweaking of the masking object (you need the inverse of the polygon over the extent of your raster).

    library(cartography)
    library(sf)
    library(raster)
    
    mtq <- st_read(system.file("gpkg/mtq.gpkg", package="cartography"), quiet = TRUE) %>% 
      dplyr::summarise() # dissolve internal boundaries
    
    r <- raster(matrix(rnorm(10*12), nrow=10), 
                xmn = st_bbox(mtq)["xmin"], 
                xmx= st_bbox(mtq)["xmax"], 
                ymn = st_bbox(mtq)["ymin"], 
                ymx= st_bbox(mtq)["ymax"], 
                crs = st_crs(mtq))
    
    plot(r) # original raster - full extent + range
    

    enter image description here

    # the masking object:
    mask <- st_bbox(r) %>% # take extent of your raster... 
      st_as_sfc() %>% # make it a sf object
      st_set_crs(st_crs(mtq)) %>% # in CRS of your polygon 
      st_difference(mtq) %>% # intersect with the polygon object
      st_as_sf() # interpret as sf (and not sfc) object
    
    result <- r %>% 
      mask(mask)
    
    plot(result)
    plot(st_geometry(mtq), add = T)
    

    enter image description here