如何在 shapefile (sf) 列表中设置 crs?

How to set the crs on a list of shapefiles (sf)?

我正在读取 shapefile 列表(如 sf),我想设置 crs。我知道这可以很容易地为一个 shapefile 完成(例如,),但我还没有找到针对多个或一系列 shapefile 的解决方案。

我正在阅读 list.fileslapply 的 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 来做,那么 returns 会出现以下错误:

Error in st_read.default(crs = 4326): dsn should specify a data source or filename

示例数据(即读入 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

如果我尝试使用此代码使用 lapplypurrr::map,则它不起作用。

lapply(myfiles, st_crs() <- 4326)

错误

Error in st_crs() <- 4326 : invalid (NULL) left side of assignment

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

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

我们只需要 lapply 中的 lambda(匿名)函数来进行赋值,然后 return 原始对象

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
        })

另一种选择(如上面@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