将不与多边形重叠的网格单元设置为NA

Question

我想提取那些与多边形重叠的网格单元。我使用了“crop”和“mask”，但它无法正常工作并且没有正确检测到所有网格单元。

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

结果：

但是有更多的网格被多边形重叠：

Answer 1

您可以使用

mask

包中的

terra

获得所需的输出，因为此函数具有

touches

参数（在

raster

包

mask

中不可用）。由于您没有添加最小的可重现示例，因此这里有一个。

示例数据

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)

带terra包的面具

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

带光栅包的掩码

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