如何使用空间多边形(shapefile)提取NetCDF中的O3数据替代经纬度范围筛选

以下为R语言环境下的完整操作流程：

依赖安装 #

首先安装所需的空间处理包：

install.packages(c("ncdf4", "terra", "sf"))

操作步骤 #

• 读取目标数据并对齐坐标系

library(sf)
library(terra)
library(ncdf4)

# 读取多边形shapefile，替换为你本地的文件路径
shp <- st_read("your_region.shp")
# 转换为WGS84坐标系，和NetCDF常用的经纬度坐标系对齐
shp_wgs84 <- st_transform(shp, crs = "EPSG:4326")

# 读取NetCDF文件
nc_obj <- nc_open("your_o3_data.nc")
# 提取O3变量、经纬度数组，变量名可根据你的NetCDF实际字段调整
o3_array <- ncvar_get(nc_obj, "O3")
lon_arr <- ncvar_get(nc_obj, "lon")
lat_arr <- ncvar_get(nc_obj, "lat")
nc_close(nc_obj)

• 将O3数组转换为空间栅格对象

# 构造栅格对象，多时间维度数据会自动识别为多层栅格
o3_rast <- rast(
  x = o3_array,
  extent = ext(range(lon_arr), range(lat_arr)),
  crs = "EPSG:4326"
)

• 多边形掩膜提取

# 转换shapefile格式适配栅格操作
shp_terra <- vect(shp_wgs84)
# 裁剪并掩膜，仅保留多边形覆盖范围内的O3数据
o3_extracted <- crop(o3_rast, shp_terra, mask = TRUE)

• 结果导出

# 导出为带经纬度的表格
o3_df <- as.data.frame(o3_extracted, xy = TRUE)
write.csv(o3_df, "o3_polygon_extracted.csv", row.names = FALSE)

# 导出为新的NetCDF文件
writeCDF(o3_extracted, "o3_polygon_extracted.nc", 
         varname = "O3", unit = "μg/m³", # 单位替换为你数据的实际单位
         overwrite = TRUE)

注意事项 #

若你的NetCDF使用的是投影坐标系而非WGS84经纬度，需将st_transform的crs参数改为NetCDF对应的坐标系编码，保证两者坐标系完全一致，避免提取位置偏移。多时间维度的O3数据无需额外循环处理，上述流程会自动对每个时间层执行掩膜操作。

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