如何调整 r 中的光栅和 shapefile 投影以使其适合裁剪？

Question

我是地理空间数据的新手，正在尝试使用shapefile[=51=裁剪tif文件raster对象] 通过引用 https://www.youtube.com/watch?v=UP2Za1TizOc.

我已参考上述视频尝试了以下代码，但似乎存在投影问题：

library(tidyverse)
library(sf)
library(stars)

ind_outline <- sf::st_read("path\polymap15m_area.shp")
ind_outline


Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2815341 ymin: 2177499 xmax: 5678865 ymax: 5444567
Projected CRS: LCC_WGS84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((5547296 2230982, ...
2   0       1875 POLYGON ((5560180 2232030, ...
3   0       1875 POLYGON ((5549993 2253154, ...
4   0       1875 POLYGON ((5542651 2256150, ...
5   0       1875 POLYGON ((5533962 2260494, ...
6   0       1875 POLYGON ((5523175 2264240, ...
7   0       1875 POLYGON ((3223295 2294948, ...
8   0       1875 POLYGON ((5502051 2315325, ...
9   0       1875 POLYGON ((5522126 2328209, ...
10  0       1875 POLYGON ((5480027 2338995, ...

shapefile link: https://github.com/johnsnow09/covid19-df_stack-code/blob/main/polymap15m_area.shp

ind_outline %>% 
        st_as_sf() %>% 
        ggplot() +
        geom_sf()

ind_region_stars <- stars::read_stars("../gt30e060n40.tif")
ind_region_stars

stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                 Min. 1st Qu. Median     Mean 3rd Qu. Max.
gt30e060n40.tif   130     793   1052 1302.186    1648 4795
dimension(s):
  from   to offset       delta refsys point values x/y
x    1 4800     60  0.00833333 WGS 84 FALSE   NULL [x]
y    1 6000     40 -0.00833333 WGS 84 FALSE   NULL [y]

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c()

问题： 当我尝试将一个叠加在另一个上时，它不起作用

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c() +
        geom_sf(data = ind_outline, alpha = 0)

裁剪错误：

ind_region_stars_cropped <- st_crop(ind_region_stars, ind_outline) 

Error in st_crop.stars(ind_region_stars, ind_outline) : for cropping, the CRS of both objects have to be identical

更新：

st_crs(ind_outline)

Coordinate Reference System:
  User input: LCC_WGS84 
  wkt:
PROJCRS["LCC_WGS84",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",24,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",80,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",12.472944,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",35.172806,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    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]]]]

st_crs(ind_region_stars)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Answer 1

要将多边形设置为与栅格相同的 crs，则应使用

明确地执行此操作

ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))