如何在R中从CORINE提取葡萄园土地覆盖并解决投影冲突问题

一、确认CORINE提取葡萄园的步骤是否正确 #

CORINE Land Cover的二级分类里，值15确实对应葡萄园（属于农业区域下的子类），你的提取逻辑没问题，可以用以下代码验证并确认结果：

library(terra)
# 导入CORINE的TIFF文件
clc <- rast("你的CORINE文件路径.tif")
# 提取值为15的葡萄园区域
vineyard <- clc == 15
# 查看分类统计，确认15类的数量是否合理
freq(vineyard)

小提示：最好核对下你下载的CORINE数据元数据，确认分类编码对应关系，不同年份的CORINE可能有细微调整，但15作为葡萄园的编码是比较稳定的。

二、解决[crop] extents do not overlap错误 #

这个报错的核心原因就是两个栅格的坐标系（CRS）不统一，或者空间范围完全没重叠，按以下步骤排查解决：

1. 统一所有数据的CRS #

你的意大利行政区划shapefile用的是WGS 84 / UTM zone 32N（EPSG:32632），但CORINE通常用ETRS89 LAEA投影（EPSG:3035），ERA5的NetCDF一般是WGS84地理坐标系（EPSG:4326），三者CRS必须一致才能操作。

先查看各数据的CRS：

# 查看葡萄园栅格的CRS
crs(vineyard)
# 查看掩膜（你的mask）的CRS，如果是RasterLayer用crs(mask)，如果是矢量用sf包的st_crs(mask)

然后统一到同一个CRS，比如统一到shapefile的EPSG:32632：

# 如果vineyard的CRS和mask不一样，转换vineyard的坐标系
vineyard <- project(vineyard, crs(mask))
# 如果你的mask是矢量shapefile，先转成栅格掩膜（确保分辨率和vineyard一致）
mask_rast <- rast(mask, resolution = res(vineyard))
mask_rast <- rasterize(mask, mask_rast)

2. 检查空间范围是否真的重叠 #

如果CRS统一后还是报错，就检查两者的空间范围：

# 查看vineyard的范围
ext(vineyard)
# 查看mask的范围
ext(mask)

要是确实没重叠，大概率是你下载的CORINE数据不是意大利区域的，或者mask的范围超出了CORINE的覆盖范围。这时候先把vineyard裁剪到意大利的大致范围再操作：

# 从mask（意大利行政区划）获取范围
italy_ext <- ext(mask)
# 裁剪vineyard到意大利范围
vineyard_italy <- crop(vineyard, italy_ext)
# 再用mask做掩膜
vineyard_masked <- mask(vineyard_italy, mask)

三、把葡萄园栅格重采样到ERA5的0.1度分辨率 #

ERA5的0.1度分辨率是基于WGS84地理坐标系（EPSG:4326）的，重采样步骤如下：

# 导入ERA5温度数据，获取它的栅格模板（分辨率、范围、CRS）
era5_temp <- rast("你的ERA5温度NetCDF路径.nc")
# 先把掩膜后的葡萄园栅格转成和ERA5一致的CRS
vineyard_wgs84 <- project(vineyard_masked, crs(era5_temp))
# 重采样到ERA5的分辨率和范围，用sum方法是为了后续计算每个网格的葡萄园面积
vineyard_resampled <- resample(vineyard_wgs84, era5_temp, method = "sum")

四、按葡萄园占比加权计算市政温度 #

要实现这个需求，核心是先计算每个市政内的葡萄园占比，再用这个占比加权温度，具体步骤：

library(sf)
# 导入意大利市政矢量文件
municipalities <- st_read("你的市政shapefile路径.shp")
# 转成和ERA5一致的CRS
municipalities_wgs84 <- st_transform(municipalities, crs(era5_temp))

# 1. 把市政矢量栅格化到ERA5的分辨率，用市政ID作为字段
muni_rast <- rasterize(municipalities_wgs84, era5_temp, field = "你的市政ID字段名")

# 2. 计算每个市政的葡萄园面积和总面积
# 先计算每个ERA5网格的实际面积
cell_area <- cellSize(vineyard_resampled)
# 计算每个网格的葡萄园面积
vineyard_area <- vineyard_resampled * cell_area
# 聚合到市政级别，分别计算葡萄园总面积和市政总面积
muni_vineyard_area <- zonal(vineyard_area, muni_rast, fun = "sum")
muni_total_area <- zonal(cell_area, muni_rast, fun = "sum")
# 合并数据并计算葡萄园占比
muni_vineyard_ratio <- merge(muni_vineyard_area, muni_total_area, by = "你的市政ID字段名")
muni_vineyard_ratio$vineyard_ratio <- muni_vineyard_ratio$sum.x / muni_vineyard_ratio$sum.y

# 3. 计算加权温度
# 方式一：先算市政平均温度，再用占比加权（简单版）
muni_avg_temp <- zonal(era5_temp, muni_rast, fun = "mean")
muni_weighted_temp <- merge(muni_avg_temp, muni_vineyard_ratio, by = "你的市政ID字段名")
# 这里的加权逻辑可以根据需求调整，比如只针对葡萄园区域的温度，或者整个市政的加权
muni_weighted_temp$weighted_temp <- muni_weighted_temp$mean * muni_weighted_temp$vineyard_ratio + muni_weighted_temp$mean * (1 - muni_weighted_temp$vineyard_ratio)

# 方式二：更精确的加权——用每个网格的葡萄园占比加权该网格的温度，再聚合到市政（推荐）
# 先计算每个网格的葡萄园占比（0-1）
grid_vineyard_ratio <- vineyard_resampled / max(vineyard_resampled, na.rm = TRUE)
# 用占比加权每个网格的温度
weighted_temp_grid <- era5_temp * grid_vineyard_ratio
# 聚合到市政，得到加权后的温度总和（如果是只考虑葡萄园区域的温度，直接用sum即可）
muni_precise_weighted_temp <- zonal(weighted_temp_grid, muni_rast, fun = "sum")

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