R 栅格和多边形不重叠的投影问题，即使它们重叠

Question

我需要通过多边形（USGS HUC 数据）裁剪栅格（Sentinel 2 图像）。栅格的 crs 是：弃用的 Proj.4 表示：+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs。多边形的 crs 是 "GEOGCRS["NAD83", DATUM["北美基准 1983", ……”

我已经尝试将多边形空间转换为栅格的 crs，还尝试将栅格投影到多边形的 crs 中。这两个项目的范围重叠但单位不同（UTM 与纬度/经度）。

改变投影使两者重叠的技巧是什么？

代码如下：

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

Answer 1

多边形作为 sf 引入，但空间变换要求它是 sp。将 sf 转换为 sp 允许进行空间转换。

UpCarson_sp <- as(UpCarson, "Spatial")
##crs taken directly from the raster crs
UpCarson_project <- spTransform(UpCarson, CRSobj = "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs")