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

Question

I am reading in a list of shapefiles (as sf) and I want to set the crs. I know this can easily be done for one shapefile (e.g., here), but I have not found a solution for multiple or a list of shapefiles.

I am reading in the shapefiles with list.files and lapply.

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)

If reading in 1 shapefile, then you can simply set the crs by:

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

However, if I do it with lapply, then it returns the following error:

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

Sample Data (i.e., list of shapefiles after being read into R)

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

For 1 shapefile in the list, I can easily set the crs by doing:

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

If I try to use lapply or purrr::map using this code, then it does not work.

lapply(myfiles, st_crs() <- 4326)

Error

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

I also tried using st_transform in similar ways, but it also didn't work.

What I want to do is to use something like st_crs(shapefiles[[1]]) <- 4326 but apply it to each element (i.e., shapefile) in the list. Also, I'd like to avoid using a for loop. I would rather use an apply or purrr function.

akrun · Accepted Answer

We just need a lambda (anonymous) function in lapply to do the assignment and then return the original object

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

-output

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

Or using map

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

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

Answers (2)

Related Questions