使用 R 中的 shapefile 提取栅格值

Question

我正在尝试使用 R 提取 shapefile 中栅格的像素单元。

我正在使用包 raster 中的函数 extract。

栅格是格式为 .tiff 的栅格，形状文件是具有四个多边形的形状文件。

r.vals <- extract(raster,sdata[1,])

我在 r.vals 中得到一个空列表，

我认为这是投影问题。

这是栅格的投影：

> show(raster)
class      : RasterLayer 
dimensions : 1801, 3600, 6483600  (nrow, ncol, ncell)
resolution : 0.1, 0.1  (x, y)
extent     : -0.05, 359.95, -90.05, 90.05  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : 1950_01.tif 
names      : X1950_01 
values     : -1.387779e-17, 0.0448938  (min, max)

这是形状的投影：

> show(sdata)
Simple feature collection with 4 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -113.2154 ymin: 33.26956 xmax: -109.1728 ymax: 35.87181
Geodetic CRS:  WGS 84
# A tibble: 4 x 3
  Station      Area_m2                                                                                geometry
    <dbl>        <dbl>                                                                      <MULTIPOLYGON [°]>
1 9502000 16087018587. (((-109.653 34.10926, -109.6533 34.10926, -109.6533 34.10871, -109.653 34.10871, -10...
2 9510000 16004291893. (((-112.0876 34.57458, -112.0876 34.57431, -112.0879 34.57431, -112.0879 34.57458, -...
3 9498500 11107332300. (((-109.653 34.10927, -109.6533 34.10927, -109.6533 34.10871, -109.653 34.10871, -10...
4 9508500 15216249842. (((-112.0876 34.57459, -112.0876 34.57431, -112.0879 34.57431, -112.0879 34.57459, -...

我尝试使用 rgdal 中的 st_transform 函数

st_transform(sdata, CRS(projection( raster)))

具有以下结果

Simple feature collection with 4 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -110 ymin: 33 xmax: -110 ymax: 36
CRS:           GEOGCRS["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]],
    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]]]]
# A tibble: 4 x 3
  Station      Area_m2                                                                                geometry
*   <dbl>        <dbl>                                                                      <MULTIPOLYGON [°]>
1 9502000 16087018587. (((-110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, ...
2 9510000 16004291893. (((-112 35, -112 35, -112 35, -112 35, -112 35)), ((-112 35, -112 35, -112 35, -112 ...
3 9498500 11107332300. (((-110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, -110 34, ...
4 9508500 15216249842. (((-112 35, -112 35, -112 35, -112 35, -112 35)), ((-112 35, -112 35, -112 35, -112 ...

，但 r.vals 仍然为空。

知道我做错了什么吗？