提取与多边形重叠网格单元失败，crop和mask未识别全部网格

问题根源 #

• raster::mask()默认仅保留网格中心落在多边形范围内的单元，无法捕获「仅部分区域与多边形重叠」的网格。
• 代码混用sp与sf对象，易引发隐性的坐标参考系(CRS)匹配错误。
• 未验证栅格数据与矢量多边形的CRS是否完全一致，空间计算会因此出现偏差。

修正方案与代码 #

# 加载所需包
library(sf)
library(raster)
library(ggplot2)
library(maps)

# 读取美国边界（用sf统一处理矢量）
Us <- st_read(dsn = "C:/shapefile", layer = "boundaries")

# 读取臭氧栅格数据
data <- brick(path = "C:/", pattern = "nc", varname = "o3")

# 确保栅格与US边界CRS一致，转换栅格CRS（如果需要）
if (!compareCRS(data, Us)) {
  data <- projectRaster(data, crs = st_crs(Us)$wkt)
}

# 裁剪+保留美国境内的栅格
data_us <- crop(data, Us)
data_us2 <- mask(data_us, Us)

# 读取urban多边形并统一CRS
urban <- st_read("C:/shapefile_urban")
# 转换为与栅格一致的CRS（直接用栅格的CRS，不要硬编码）
urban_T <- st_transform(urban, crs = st_crs(data_us2))

# 生成栅格的网格多边形（每个网格单元转为sf多边形）
grid_poly <- st_as_sf(data_us2, as_points = FALSE, merge = FALSE)

# 判断哪些网格与urban多边形有重叠
intersect_idx <- st_intersects(grid_poly, urban_T, sparse = FALSE)
# 保留至少与一个urban多边形重叠的网格
data_us4 <- data_us2[rowSums(intersect_idx) > 0]

# 转换为数据框用于绘图
df_i <- as.data.frame(data_us4, xy = TRUE)
names(df_i) <- c("Longitude", "Latitude", "value")
df_i <- na.omit(df_i)

# 绘图
ggplot(df_i, aes(x = Longitude, y = Latitude, fill = value)) +
  geom_tile() +
  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)

关键调整说明 #

• 用sf::st_read替代readOGR，统一矢量数据处理流程，避免sp/sf兼容性问题。
• 动态匹配CRS：不再硬编码WGS84，而是直接对齐栅格数据的CRS，确保空间计算准确。
• 生成网格多边形并通过st_intersects判断重叠：捕获所有与urban多边形有接触的网格单元，而非仅中心落在多边形内的单元。

内容的提问来源于stack exchange，提问作者Arwen