R投影问题：栅格与多边形实际重叠却无法匹配裁剪

我想用USGS HUC多边形裁剪Sentinel-2栅格影像。栅格的CRS是弃用的Proj.4格式：+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs；多边形的CRS是NAD83地理坐标系。我试过把多边形转成栅格的CRS，也试过把栅格投影到多边形的CRS，二者实际空间范围是重叠的，但执行裁剪时一直报错extents do not overlap。我的代码如下：

UpCarson <- read_sf("UpperCarson.shp")
crs(UpCarson)
> crs(Mar17_project)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=NAD83 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 

Mar17 <- raster("S2A2A_20230317_070_sen2r_NDSI_10.tif")
crs(Mar17)
> crs(Mar17)
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]],
        ID["EPSG",16010]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 

Mar17_project <- projectRaster(Mar17, crs = crs(UpCarson))
crs(Mar17_project)

> crs(Mar17_project)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=NAD83 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 

##Crop the raster by the polygon
Mar17_crop <- crop(Mar17_project, UpCarson)
> Mar17_crop <- crop(Mar17_project, UpCarson)
Error in .local(x, y, ...) : extents do not overlap

问题根源 #

报错核心是WGS84与NAD83基准面的微小偏移，即使完成投影转换，坐标的细微差异会让程序误判范围不重叠；另外用projectRaster转换栅格时，若未指定正确参数，也可能导致栅格范围计算偏差。

正确处理步骤 #

方案1：将多边形转换到栅格的UTM坐标系（推荐，避免栅格重采样精度损失） #

优先转换矢量数据，因为栅格重采样会损失精度：

library(sf)
library(raster)

# 读取多边形与栅格
UpCarson <- read_sf("UpperCarson.shp")
Mar17 <- raster("S2A2A_20230317_070_sen2r_NDSI_10.tif")

# 将多边形转换到栅格的CRS（直接调用栅格的crs对象，避免手动输入Proj4字符串）
UpCarson_utm <- st_transform(UpCarson, crs = crs(Mar17))

# 先裁剪到多边形范围，再掩膜去除外部像素
Mar17_crop <- crop(Mar17, extent(UpCarson_utm))
Mar17_masked <- mask(Mar17_crop, UpCarson_utm)

方案2：修复栅格投影到NAD83的问题 #

若必须转换栅格，需明确指定基准面转换参数，使用官方EPSG代码更可靠：

library(sf)
library(raster)

# 读取数据
Mar17 <- raster("S2A2A_20230317_070_sen2r_NDSI_10.tif")
UpCarson <- read_sf("UpperCarson.shp")

# 用EPSG:4269指定NAD83地理坐标系（官方标准，比手动Proj4更准确）
nad83_crs <- st_crs(4269)

# 投影栅格时使用WKT格式CRS，指定重采样方法
Mar17_project <- projectRaster(Mar17, crs = nad83_crs$wkt, method = "bilinear")

# 检查范围是否重叠，确认数值匹配
print(extent(Mar17_project))
print(st_bbox(UpCarson))

# 执行裁剪和掩膜
Mar17_crop <- crop(Mar17_project, extent(UpCarson))
Mar17_masked <- mask(Mar17_crop, UpCarson)

关键注意点 #

• 用sf包的st_transform处理矢量数据，比旧的sp包方法更稳定，支持现代WKT格式CRS。
• 优先使用EPSG代码指定CRS（如NAD83用4269，UTM10N WGS84用32610），避免手动编写Proj4字符串减少错误。
• 裁剪建议分两步：先crop缩小栅格范围，再mask保留多边形内像素，便于排查问题。
• 转换后务必检查extent（栅格）和st_bbox（矢量）的数值，确认范围确实重叠。

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