如何在 shapefile (sf) 列表中设置 crs？答案

【问题标题】：How to set the crs on a list of shapefiles (sf)?如何在 shapefile (sf) 列表中设置 crs？
【发布时间】：2021-09-17 09:03:37
【问题描述】：

我正在阅读 shapefile 列表（如 sf），我想设置 crs。我知道这可以很容易地为一个 shapefile 完成（例如，here），但我没有找到将它应用于多个或一系列 shapefile 的任何帮助。

我正在阅读带有 list.files 和 lapply 的 shapefile。

library(sf)

# Example (Not Run)
temp = list.files(path = "/data",
                  pattern = "*WEST.shp$",
                  full.names = TRUE)

# Then, apply the st_read function to the list of shapefiles (Not Run).
myfiles = lapply(temp, st_read)

如果您正在阅读 1 个 shapefile，那么您可以通过以下方式简单地设置 crs：

boundary <- st_read(
         "/data/tsugWEST.shp", crs = 4326)

但是，如果我尝试使用lapply 执行此操作，则会返回以下错误：

st_read.default(crs = 4326) 中的错误：dsn 应指定数据来源或文件名

样本数据（即读入 R 后的 shapefile 列表）

shapefiles <- list(abieWEST.shp = structure(list(AREA = 1, PERIMETER = 0, ABIELASI_ = 0, 
    ABIELASI_I = 0, CODE = 0L, geometry = structure(list(structure(list(
        structure(c(-126.623274277836, -122.269490649612, -117.311014850802, 
        -114.529430866103, -113.803800261399, -112.836292788461, 
        -108.603447594355, -103.765910229662, -102.193710586137, 
        -101.105264679081, -101.58901841555, -107.515001687299, 
        -115.255061470807, -117.552891719036, -125.897643673132, 
        -130.49330416959, -141.740578542501, -146.940931209546, 
        -144.643100961316, -136.056472138987, -126.623274277836, 
        63.8736898397207, 61.9386748938436, 55.5289378856255, 
        55.4079994515082, 54.0776766762177, 52.9892307691618, 
        48.9982624432902, 45.3701094197706, 42.8304023033068, 
        36.7834805974407, 28.3177902092282, 23.8430681468874, 
        24.3268218833566, 27.834036472759, 35.5740962562675, 
        49.2401393115248, 55.6498763197429, 63.1480592350168, 
        65.4458894832459, 65.4458894832459, 63.8736898397207), .Dim = c(21L, 
        2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
        input = NA_character_, wkt = NA_character_), class = "crs"), class = c("sfc_POLYGON", 
    "sfc"), precision = 0, bbox = structure(c(xmin = -146.940931209546, 
    ymin = 23.8430681468874, xmax = -101.105264679081, ymax = 65.4458894832459
    ), class = "bbox"))), row.names = c(NA, -1L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(AREA = NA_integer_, 
PERIMETER = NA_integer_, ABIELASI_ = NA_integer_, ABIELASI_I = NA_integer_, 
CODE = NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"))), acerWEST.shp = structure(list(AREA = 1, 
    PERIMETER = 0, ABIELASI_ = 0, ABIELASI_I = 0, CODE = 0L, 
    geometry = structure(list(structure(list(structure(c(-108.146180374448, 
    -125.155521463775, -140.514702596675, -134.929545821075, 
    -128.963582901684, -126.805681420203, -124.647779938721, 
    -123.886167651139, -120.458912357021, -120.458912357021, 
    -114.746820200157, -111.573435668566, -108.526986518239, 
    -107.003761943076, -107.003761943076, -104.591989699066, 
    -104.084248174012, -102.561023598848, -104.845860461594, 
    -102.053282073794, -98.1182852546211, -94.310223816712, -108.146180374448, 
    15.5862714176048, 37.6730277574772, 57.9826887596588, 60.013654859877, 
    55.9517226594406, 56.9672057095497, 55.5709165156497, 55.6978518969134, 
    56.0786580407043, 54.3015627030134, 51.2551135526862, 47.827858258568, 
    47.9547936398317, 46.5585044459317, 45.416086014559, 42.2427014829681, 
    37.6730277574772, 37.1652862324227, 34.1188370820955, 30.3107756441864, 
    30.5646464067137, 26.2488434437501, 15.5862714176048), .Dim = c(23L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
        input = NA_character_, wkt = NA_character_), class = "crs"), class = c("sfc_POLYGON", 
    "sfc"), precision = 0, bbox = structure(c(xmin = -140.514702596675, 
    ymin = 15.5862714176048, xmax = -94.310223816712, ymax = 60.013654859877
    ), class = "bbox"))), row.names = c(NA, -1L), class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(AREA = NA_integer_, 
PERIMETER = NA_integer_, ABIELASI_ = NA_integer_, ABIELASI_I = NA_integer_, 
CODE = NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"))))

对于列表中的 1 个 shapefile，我可以通过以下方式轻松设置 crs：

st_crs(shapefiles[[1]]) <- 4326

如果我尝试使用此代码使用lapply 或purrr::map，那么它不起作用（虽然我没想到它会这样工作）。

lapply(myfiles, st_crs() <- 4326)

错误

st_crs() 中的错误

我也尝试过以类似的方式使用st_transform，但也没有用。

我想做的是使用st_crs(shapefiles[[1]]) <- 4326 之类的东西，但将其应用于列表中的每个元素（即 shapefile）。另外，我想避免使用 for 循环。我宁愿使用apply 或purrr 函数。

【问题讨论】：

我目前无权访问sf，但您能否尝试将crs 作为“…FUN 的可选参数”（其中FUN 是@ 987654343@) 在lapply，比如lapply(temp, st_read, crs = 4326)
@Henrik 是的，这实际上也有效！我得到与下面相同的结果。我没想到把它放在st_read之外。
太棒了。然后可以在阅读文件时一口气完成。无需第二次lapply。

标签： r gis sf

【解决方案1】：

我们只需要lapply 中的 lambda（匿名）函数来进行赋值，然后返回原始对象

shapefiles2 <- lapply(shapefiles, function(x) {st_crs(x) <- 4326; x})

-输出

shapefiles2
$abieWEST.shp
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -146.9409 ymin: 23.84307 xmax: -101.1053 ymax: 65.44589
Geodetic CRS:  WGS 84
  AREA PERIMETER ABIELASI_ ABIELASI_I CODE                       geometry
1    1         0         0          0    0 POLYGON ((-126.6233 63.8736...

$acerWEST.shp
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -140.5147 ymin: 15.58627 xmax: -94.31022 ymax: 60.01365
Geodetic CRS:  WGS 84
  AREA PERIMETER ABIELASI_ ABIELASI_I CODE                       geometry
1    1         0         0          0    0 POLYGON ((-108.1462 15.5862...

或者使用map

library(purrr)
map(shapefiles, ~ {
              st_crs(.x) <- 4326
              .x
        })

【讨论】：

我尝试过与lapply 类似的东西，但在我的函数中缺少; x，这一切都不同了。那究竟是在做什么？我可以看到它保留了所有原始信息，而没有它它只为列表中的每个元素返回 4326。
@AndyBrown 在您的帖子中，您只是在调用st_crs()，内部没有任何参数。这里的x 是您正在执行分配的列表元素
@AndyBrown 与在任何其他列表操作中进行赋值相同。例如。假设我们要将 data.frame 列表的第一列更改为整数。 lst1 <- list(iris, iris); lapply(lst1, function(x) {x[[1]] <- as.integer(x[[1]]); x})

【解决方案2】：

另一个选项（如上面@Henrik 所建议的）是将crs 作为 FUN 的可选参数（即st_read）。

library(sf)

# Example (Not Run)
temp = list.files(path = "/data",
                  pattern = "*WEST.shp$",
                  full.names = TRUE)

myfiles = lapply(temp, st_read, crs = 4326)

输出

Reading layer `abieWEST' from data source 
  `/data/abieWEST.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -146.9409 ymin: 23.84307 xmax: -101.1053 ymax: 65.44589
Geodetic CRS:  WGS 84

Reading layer `acerWEST' from data source 
  `/data/acerWEST.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -140.5147 ymin: 15.58627 xmax: -94.31022 ymax: 60.01365
Geodetic CRS:  WGS 84

【讨论】：