Question

I'm trying to loop through folders in a directory while reading and assigning files to variables in R.

The files I want to assign to variables are shapefiles so I use the function readOGR from the rgdal package. The purpose is to later merge all shapefiles belonging to a particular species. The hierarchy/structure of the directory is type1 > species > ids. The shapefiles look like id.shp etc.

Example shapefiles can be downloaded here and here

#code
setwd("~/type1/")

#Extract ids belonging to $species. Later use in readOCR function

sp_id <- function(species){
 wd = "~/type1/"
 list_shp <- list.files(path=paste(wd,species,sep='/'), full.names = F, recursive = F, include.dirs = F)
 vec <- character() 

  for (shp in list_shp){
   y <- unlist(strsplit(shp, '\\.', perl=T))
   vec <- unique(c(vec,y[1]))
  }

 #"1905" "4279"

#Extract dirs for where to perform readOCR function

  list_dir <- list.dirs(path=paste(wd,species,sep='/'), full.names = F, recursive = F)

   for (id in list_dir){
    setwd(id)
    print(getwd())
   }

#"~/type1/speciesX1/1905"
#"~/type1/speciesX1/4279"

    for (i in vec){
     assign(paste("", i, sep=""), readOGR(".", i))
    break
    }
 }

sp_id('speciesX1')

[1] "~/type1/speciesX1/1905"
OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "1905"
with 10 features and 3 fields
Feature type: wkbPolygon with 2 dimensions
[1] "~/type1/speciesX1/4279"
Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv) : 
Cannot open layer

The problem is that the code only executes readOGR for one shapefile in one of the dirs, seems to change dir again but doesn't execute the last readOGR.

Was it helpful?

Solution

Your function seems to be broken in many different sections, but you could just do:

library(rgdal)

setwd("~/type1/")
species <- 'speciesX1'
list_shp <- list.files(path=species, pattern="*.shp", full.names = TRUE,
                       recursive = TRUE, include.dirs = FALSE)
shp_objects <- lapply(list_shp, function(x) {readOGR(dsn=x, 
                                                     layer=ogrListLayers(x))})

This will give you a list of SpatialPolygonDataframe objects that you can then merge.

OTHER TIPS

I could not get ogrInfo() to read those files getting a missing layer error). This provided a method to read and get some of the attributes. (A Mac will take duplicate names of directories and append "(n)" to them, so the identically named but different files were iho.zip and iho.zip:

library(sp)

ca3 = readShapeSpatial("~/Downloads/iho/iho.shp")
ca3 = readShapeSpatial("~/Downloads/iho(2)/iho.shp")

> attributes(ca3)$data
                               name   id mrgid
0 Mediterranean Sea - Western Basin 28Aa  4279
1               Strait of Gibraltar  28a  3346
2                       Alboran Sea  28b  3324
3                      Balearic Sea  28c  3322
4                      Ligurian Sea  28d  3363
5                    Tyrrhenian Sea  28e  3386
> attributes(ca2)$data
                               name   id mrgid
0 Mediterranean Sea - Western Basin 28Aa  4279
1 Mediterranean Sea - Eastern Basin 28Bb  4280
2               Strait of Gibraltar  28a  3346
3                       Alboran Sea  28b  3324
4                      Balearic Sea  28c  3322
5                      Ligurian Sea  28d  3363
6                    Tyrrhenian Sea  28e  3386
7                      Adriatic Sea  28g  3314
8                        Ionian Sea  28f  3351
9                        Aegean Sea  28h  3315
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top