1

I want to extract those grid cells that are overlapped with the polygons. I used "crop" and "mask", but it not working properly and did not detect all the grid cells correctly.

#US shapefile
Us <- readOGR(dsn ="C:/shapefile",layer="boundaries") # read shapefile

#ozone ncdf file
data <- brick(path="C:/",pattern = "nc", varname = "o3")
data_us <- crop(data, Us)
data_us2 <- mask(data_us, Us)

#urban shapefile
urban <- readOGR("C:/shapefile_urban") # read shapefile
urban_T <- st_as_sf(spTransform(urban, CRS("+proj=longlat +datum=WGS84 +no_defs 
+ellps=WGS84 +towgs84=0,0,0"))) # Transform projection

#mask with urban
data_us3 <- crop(data_us2, urban_T)
data_us4 <- mask(data_us3, urban_T)

df_i <- as.data.frame(data_us4, xy = TRUE)
names(df_i) <- c("Longitude", "Latitude", "value")
df_i <- na.omit(df_i)

#Plot
ggplot(df_i, aes(x = Longitude, y = Latitude, fill = value)) +
geom_tile() +
#add US map layer with state boundaries
geom_map(data = map_data("state"), map = map_data("state"),
     aes(long, lat, map_id = region), 
     color = "black", fill = NA,size = 0.5) +
expand_limits(x = map_data("state")$long, y = map_data("state")$lat)

Result: enter image description here

But there are more grids that overlapped by polygons: enter image description here

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
Arwen
  • 77
  • 4

1 Answers1

1

You can obtain the desired output using mask from the terra package, since this function has the touches argument (which is not available in the raster package mask). Since you did not add a minimal reproducible example, here is one.

Example data

library(terra)
library(raster)
library(sf)

# Create example data
rast1 <- rast(matrix(rnorm(225,50,20), 15))

vect1 <- vect("POLYGON ((6.5 6.5, 6.2 7.8, 7.3 7.3, 7.5 6.2, 6.5 6.5))")

# Show example data
plot(rast1)
plot(vect1, add = T)

example

Mask with terra package

# touches = T is the default, but adding it to make it explicit
plot(terra::mask(rast1, vect1, touches = T))

maskw/terra

Mask with raster package

rast2 <- raster::raster(rast1)
# Convert vect1 to sf to use raster::mask
plot(raster::mask(rast2, st_as_sf(vect1)))

maskw/raster