Как настроить проекцию растров и шейп-файлов в r, чтобы сделать их пригодными для обрезки?

Я новичок в геопространственных данных и пытаюсь обрезать объект tif file 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, ...

Ссылка на шейп-файл: 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) 

Ошибка в st_crop.stars (ind_region_stars, ind_outline): для обрезки CRS обоих объектов должны быть идентичны

ОБНОВЛЕНИЕ:

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]]

person ViSa    schedule 29.06.2021    source источник
comment
вам просто нужно перепроецировать их в одну и ту же желаемую CRS.   -  person dww    schedule 29.06.2021
comment
да, но я не смог понять lcc .. Думаю, это сработает: stackoverflow.com/questions/30287065/   -  person ViSa    schedule 29.06.2021
comment
если вы хотите установить для полигонов тот же crs, что и растр, сделайте это явно - т.е. перепроецируйте на ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))   -  person dww    schedule 29.06.2021
comment
@dww Я пробовал код ниже, но получаю ошибку; crs <- CRS("+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=126 +datum=WGS84") ind_outline_crs <- SpatialPoints(ind_outline, proj4string=crs)   -  person ViSa    schedule 29.06.2021
comment
получение ошибки не очень полезно. Прочтите https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example   -  person dww    schedule 29.06.2021
comment
спасибо большое @dww, ваш код работал. Если вы добавите это в ответ, я приму это и закрою этот пост.   -  person ViSa    schedule 29.06.2021


Ответы (1)


Чтобы установить для полигонов тот же crs, что и для растра, вы должны сделать это явно, используя

ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))
person dww    schedule 29.06.2021