3

I have three spatial dataframes (polygons), created via this tutorial. The three spatial dataframes must be exported as shapefiles via a for loop. I created a list polygons containing the three spatial dataframes. Each spatial dataframe should be treated as a separate item that is used in the for loop. However, when I run the for loop, the first item of the polygons is used three times, resulting in 3 shapefiles that have the same coordinates (i.e. the coordinates of spatial dataframe 1). In case I access the second item of the list polygons via polygons[[2]] the information corresponding to the second item is shown nevertheless.

Each iteration in the for loop should take the corresponding item (i.e. spatial dataframe) of the list, eventually resulting in three shapefiles that are situated at different locations - How can I do this?

code:

library(sp)
library(sf)

## create spatial dataframe 1 

FID_df1 <-  c(1)
x_coord <- c(16.48438,  17.49512,  24.74609, 22.59277, 16.48438)
y_coord <- c(59.736328125, 55.1220703125, 55.0341796875, 61.142578125, 59.736328125)
xym <- cbind(x_coord, y_coord)
xym

p = Polygon(xym)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
plot(sps)

proj4string(sps) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
sps$FID_df1 <- FID_df1

## create spatial dataframe 2

FID_df2 <-  c(1)
x_coord2 <- c(18.48438,  19.49512,  26.74609, 24.59277, 18.48438)
y_coord2 <- c(61.736328125, 57.1220703125, 57.0341796875, 63.142578125, 61.736328125)
xym2 <- cbind(x_coord2, y_coord2)
xym2

p2 = Polygon(xym2)
ps2 = Polygons(list(p2),1)
sps2 = SpatialPolygons(list(ps2))
plot(sps2)

proj4string(sps2) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
sps2$FID_df2 <- FID_df2
view(sps2)

## create spatial dataframe 3

FID_df3 <-  c(1)
x_coord3 <- c(-118.48438,  -119.49512,  -126.74609, -124.59277, -118.48438)
y_coord3 <- c(35.736328125, 33.1220703125, 33.0341796875, 37.142578125, 35.736328125)
xym3 <- cbind(x_coord3, y_coord3)
xym3

p3 = Polygon(xym3)
ps3 = Polygons(list(p3),1)
sps3 = SpatialPolygons(list(ps3))
plot(sps3)

proj4string(sps3) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
sps3$FID_df3 <- FID_df3
view(sps3)

## for loop - export every spatial dataframe as shapefile

#initialize list with different spatial dataframes that need to be exported as shapefile
#assign the spatial dataframes to variables with identical names. These variables can be used in the for loop
polygons = list(sps = sps, sps2 = sps2, sps3 = sps3)

#initialize list with different values to insert into the shapefile name, relating to the different spatial dataframes
numbers = list(1, 2, 3)

#for loop - write every spatial dataframe as shapefile
for(polygon in polygons)
{ for(num in numbers)
  {
  sf_polygon <- st_as_sf(polygon)
  layer <- paste("poly", num, sep = "")
  print(layer) #verify
  sf::st_write(sf_polygon, dsn='C:/Test',layer=layer, driver = "ESRI Shapefile")
  }
}

2 Answers 2

3

There is no need to use two for loop for this case. Actually, with two for() loops you're iterating 9 times in total. Just simply use a positional index (integer 1:n) for this purpose:

for(i in seq_along(polygons)){
  layer <- paste0('poly',i,'.shp') # less characters than paste(sep='')
  print(layer)
  write_sf(st_as_sf(polygons[[i]]),layer)
}
3

You could use lapply as your iterator. Please note that I added SpatialPolygonsDataFrame to your workflow as there are issues in exporting a shapefile with an empty attribute table.

library(sp)
library(sf)

polygons <- vector(mode = "list", length = 3)
polygons[[1]] <- SpatialPolygonsDataFrame(  
  SpatialPolygons((list(Polygons(list(Polygon(cbind(
         c(16.48438, 17.49512,24.74609, 22.59277, 16.48438), 
         c(59.736328125, 55.1220703125, 55.0341796875, 61.142578125, 59.736328125)))), 
         "1"))), proj4string = CRS("EPSG:4326")), data.frame(ID=1)) 
polygons[[2]] <- SpatialPolygonsDataFrame( 
  SpatialPolygons((list(Polygons(list(Polygon(cbind(
         c(18.48438,  19.49512,  26.74609, 24.59277, 18.48438), 
         c(61.736328125, 57.1220703125, 57.0341796875, 63.142578125, 61.736328125)))), 
         "1"))), proj4string = CRS("EPSG:4326")), data.frame(ID=2)) 
polygons[[3]] <- SpatialPolygonsDataFrame(   
  SpatialPolygons((list(Polygons(list(Polygon(cbind(
         c(-118.48438,  -119.49512,  -126.74609, -124.59277, -118.48438), 
         c(35.736328125, 33.1220703125, 33.0341796875, 37.142578125, 35.7363281255)))), 
         "1"))), proj4string = CRS("EPSG:4326")) , data.frame(ID=3))

Here we use an internal iterator using the "parent.frame" object call.

lapply(polygons, function(x){
  i = parent.frame()$i
  sf::write_sf(sf::st_as_sf(polygons[[i]]),
               paste0("poly_", i, ".shp"))})

Here, we specify the iterator using seq_along.

lapply(seq_along(polygons), function(i) {
  sf::write_sf(sf::st_as_sf(polygons[[i]]),
               paste0("poly", i,".shp"))})

You could save yourself quite a bit of headache by skipping sp entirely and just using sf. First, create a single data.frame with an ID field associated with the coordinates.

library(sf)
library(dplyr)

( xy <- rbind(
data.frame(ID=1, x=c(16.48438, 17.49512,24.74609, 22.59277, 16.48438), 
y=c(59.736328125, 55.1220703125, 55.0341796875, 61.142578125, 59.736328125)),
data.frame(ID=2, x=c(18.48438,  19.49512,  26.74609, 24.59277, 18.48438), 
y=c(61.736328125, 57.1220703125, 57.0341796875, 63.142578125, 61.736328125)),
data.frame(ID=3, x=c(-118.48438,  -119.49512,  -126.74609, -124.59277, -118.48438), 
y=c(35.736328125, 33.1220703125, 33.0341796875, 37.142578125, 35.7363281255))) )

Now, use st_as_sf to coerce to sf, group_by to create MULTIPOINT data for representing each polygon and then cast into a POLYGON.

xy <- xy %>%
  st_as_sf(coords = c("x", "y"), crs = 4326) %>%
    group_by(ID) %>%
      summarise(geometry = st_combine(geometry)) %>%
        st_cast("POLYGON") 

You can now iterate through the single POLYGON feature class object rather than having a list.

lapply(unique(xy$ID), function(i){
  sf::write_sf(xy[xy$ID == i,],
  paste0("poly_", i, ".shp"))}) 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.