R: looping to save multiple shapefiles

Question

When I try to use writeOGR, plus a loop, to save my shapefiles, it doesn't do anything but give me an error message:

Error in writeOGR(plot.locationsSP_DROUGHT, dsn, layer1, driver = "ESRI Shapefile") : layer exists, use a new layer name

Essentially, I am converting each object into a CSV file, then into a shapefile, and want to save both the CSV files and the shapefiles. Here is my code fragment:

for (m in 1:500){
#First I want to save my CSV files:
drought.slice <- rotate(drought.array[m,,])
drought.vec <- as.vector(drought.slice)
length(drought.vec)
drought.df01 <- data.frame(cbind(lonlat, drought.vec))
names(drought.df01) <- c("lon", "lat", paste(dname, as.character(m), sep = "_"))
head(na.omit(drought.df01))
csvfile<-paste0("cru_drought_",m,".csv")

#Next I want to create shapefiles from the CSV files:
plot.locations_DROUGHT <- read.csv(paste0("cru_drought_",m,".csv"), stringsAsFactors = FALSE)
plot.locationsSP_DROUGHT <- SpatialPointsDataFrame(plot.locations_DROUGHT[,1:2], plot.locations_DROUGHT)
proj4string(plot.locationsSP_DROUGHT) <- CRS("+init=epsg:4326")
dsn <- layer1 <- gsub(".csv","cru_drought_",m)
writeOGR(plot.locationsSP_DROUGHT, dsn, layer1, driver="ESRI Shapefile")
}

Here is the full code I'm using:

    #Open and read the NCDF file, along with longitude and latitude
rm(list=ls())
library(lattice)
library(ncdf4)
library(chron)
library(rgdal)
library(sp)
library(raster)
library(RColorBrewer)
setwd('/Users/Neil/Dropbox/Drought Maps')
ncname <- "owda-orig"
ncfname <- paste(ncname,".nc",sep="")
dname <- "pdsi"
ncin <- nc_open(ncfname)
print(ncin)

lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(lat)

print(c(nlon, nlat))

t <- ncvar_get(ncin, "time")
nt <- dim(t)
head(t)

drought.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
#fillvalue <- ncatt_get(ncin, dname, "_FillValue")
dim(drought.array)

creation_date <- ncatt_get(ncin, 0, "creation_date")
Description <- ncatt_get(ncin, 0, "Description")

nc_close(ncin)



rotate <- function(x) t(apply(x, c(1, 2), rev))

m <- 333
drought.slice <- rotate(drought.array[m,,])
image(lon, lat, drought.slice, col = brewer.pal(10, "BrBG"))

lonlat <- expand.grid(lon, lat)
drought.vec <- as.vector(drought.slice)
length(drought.vec)

drought.df01 <- data.frame(cbind(lonlat, drought.vec))
names(drought.df01) <- c("lon", "lat", paste(dname, as.character(m), sep = "_"))
head(na.omit(drought.df01))


for (m in 1:500){
    drought.slice <- rotate(drought.array[m,,])
    drought.vec <- as.vector(drought.slice)
    length(drought.vec)
    drought.df01 <- data.frame(cbind(lonlat, drought.vec))
    names(drought.df01) <- c("lon", "lat", paste(dname,         as.character(m), sep = "_"))
    head(na.omit(drought.df01))
    csvfile<-paste0("cru_drought_",m,".csv")
    plot.locations_DROUGHT <- read.csv(paste0("cru_drought_",m,".csv"), stringsAsFactors = FALSE)
    plot.locationsSP_DROUGHT <- SpatialPointsDataFrame(plot.locations_DROUGHT[,1:2], plot.locations_DROUGHT)
    proj4string(plot.locationsSP_DROUGHT) <- CRS("+init=epsg:4326")
    dsn <- layer1 <- gsub(".csv","cru_drought_",m)
    writeOGR(plot.locationsSP_DROUGHT, dsn, layer1, driver="ESRI Shapefile")
}

Help would be most appreciated. I'm probably doing something very silly.

score 0 · Answer 1 · answered Apr 15 '18 at 00:08

I'm pretty sure your problem is in the line where you make the shapefile name. Using gsub(".csv","cru_drought_",m) would replace ".csv" with "cru_drought_" in the string m, which is just the integer you're looping over. That string isn't found, so it just assigns the integer m to dsn, so you end up trying to write shapefiles with names like "1", "2", etc. I think you just have the gsub arguments kinda scrambled.

I commented out the parts of your code that are specific to your files, and tried to just put together the file names that seem like what you're after.

for (m in 1:5){
    # drought.slice <- rotate(drought.array[m,,])
    # drought.vec <- as.vector(drought.slice)
    # length(drought.vec)
    # drought.df01 <- data.frame(cbind(lonlat, drought.vec))
    # names(drought.df01) <- c("lon", "lat", paste(dname,         as.character(m), sep = "_"))
    # head(na.omit(drought.df01))
    csvfile<-paste0("cru_drought_",m,".csv")
    # plot.locations_DROUGHT <- read.csv(paste0("cru_drought_",m,".csv"), stringsAsFactors = FALSE)
    # plot.locationsSP_DROUGHT <- SpatialPointsDataFrame(plot.locations_DROUGHT[,1:2], plot.locations_DROUGHT)
    # proj4string(plot.locationsSP_DROUGHT) <- CRS("+init=epsg:4326")
    dsn <- layer1 <- gsub(".csv","_shape",csvfile)
    # writeOGR(plot.locationsSP_DROUGHT, dsn, layer1, driver="ESRI Shapefile")
    print(dsn)
}
#> [1] "cru_drought_1_shape"
#> [1] "cru_drought_2_shape"
#> [1] "cru_drought_3_shape"
#> [1] "cru_drought_4_shape"
#> [1] "cru_drought_5_shape"

Created on 2018-04-14 by the reprex package (v0.2.0).

Thank you ever so much, this works! I have another question: it seems that this process creates separate folders for each shapefile. Is there a way to simply do this and just save the shapefiles in the working directory folder? — GIS_newb, Apr 15 '18 at 00:53
The `dsn` argument to `writeOGR` is the directory where you're saving the shapefile. Try changing to something like `writeOGR(shapefile, "folder_of_shapes", layer1, driver="ESRI Shapefile")` — camille, Apr 15 '18 at 01:51

score 0 · Answer 2 · answered Apr 15 '18 at 18:16

0

It is a bit easier to use raster::shapefile instead of writeOGR

shapefile(plot.locationsSP_DROUGHT, dsn)

answered Apr 15 '18 at 18:16

Robert Hijmans

40,301
4
55
63

R: looping to save multiple shapefiles

2 Answers2