4

I have downloaded a US federal lands data layer (fedlanp010g.shp.tar.gz) from http://nationalmap.gov/small_scale/atlasftp.html?openChapters=chpbound#chpbound

In R, I am trying to create a buffer around each polygon, However gBuffer crashes when the code below is run. I added a pause (sys.sleep(1)) to add a delay, and I received the error message "Value of SET_STRING_ELT() must be "CHARSXP" not a 'NULL'" at row 14.

Here is the code I have been using (I am processing with a loop):

library(rgeos)
library(raster)

fed = readOGR("downloaded file from above here")
i = 1
temp = fed[i,]
buff = gBuffer(temp,byid=TRUE,width=25)
result = erase(buff,temp)

for(i in 2:nrow(fed)){
    temp = fed[i,]
    buff = gBuffer(temp,byid=TRUE,width=25)
    temp_buff = erase(buff,temp)
    result = rbind(result,temp_buff)
  } 
1
  • 1
    Are you reprojecting your data? The data in your post is in "longlat NAD83" which is not a planar coordinate system and as such, is not supported by gBuffer. Jan 21, 2016 at 18:42

1 Answer 1

1

Is R actually crashing or is the code just failing? These are very different outcomes. I believe that your problem was that your data is in a geographic and not planar projection. I reprojected the data, from your link, to the USGS definition of Albers and gBuffer ran with no problems. You will need rgdal along with sp for the projection transformation.

library(rgeos)
library(sp)
library(rgdal)

setwd("D:/TMP/test")
fed <- readOGR(getwd(), "fedlanp010g")

fed <- spTransform(fed, CRS( " +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0" ) )  

fed.buf <- gBuffer(fed, byid = TRUE, width = 25)

You may be getting into trouble with how the for loop is defined and the use of erase. Since you did not provide a reproducible example it is difficult to say what the exact issue is as it could be package, code or data related.

Try this code, where each buffered polygon is stored in a list and then combined using do.call with rbind. However, in this case, I am not sure what happens with intersecting geometry.

First, add required libraries and create some example polygon data.

library(rgeos)
library(sp)

p1=Polygons(list(Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))), "1")
p2=Polygons(list(Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))), "2")
p3=Polygons(list(Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))), "3") 
fed = SpatialPolygonsDataFrame(SpatialPolygons(list(p1,p2,p3)), data.frame(ID=c("1","2","3")))

Here is a for loop that buffers each polygon, adds to a list object and then uses do.call to rbind the buffered polygons together.

result <- list()
for(i in 1:length(fed)){result[[i]] <- gBuffer(fed[i,], byid=TRUE, width = -0.05)}
result <- do.call("rbind", result)

plot(result) 

I am wondering why you are applying a for loop. You get exactly the same results just using gBuffer with the entire spatial object and the byid=TRUE argument.

fed.buf <- gBuffer(fed, byid=TRUE, width = -0.05)
plot(fed.buf)
6
  • R was originally crashing, that is why I added a loop. gBuffer is the problem. When erase is removed from the loop, same problems occur. The code is reproducible; I added into the code the the object fed comes from the downloaded shapefile. I tried running your code, but it failed with the downloaded shapefile as well.
    – user44796
    Jan 21, 2016 at 16:57
  • I didn't specify that. Yes, I reprojected into Albers Equal Area. gBuffer works for the first 13 polygons, but fails on the 14th and returns the error message that I specified above
    – user44796
    Jan 21, 2016 at 18:53
  • I am having absolutely no issues with exactly this data. I even subset the first 50 polygons and ran the for loop and it worked. I would recommend downloading the data again, making sure that your R and package versions are up to date (R 3.2.3, sp 1.2-1, rgdal 1.1-1, rgeos 0.3-15) and using "readOGR" in the rgdal package, and not maptools, to read the shapefile. Jan 21, 2016 at 19:01
  • I updated R, sp, rgdal, and rgeos to the current versions. I redownloaded the shapefile and reprojected using spTransform(). Now I receive the error message "rgeos_buffer: geometry count/id count mismatch -id changed Non.unique values when setting row.names." This occurs at fed[30,]. duplicated(row.names(fed)) returns false for all row names.
    – user44796
    Jan 21, 2016 at 20:52
  • Are you you using rstudio? Jan 22, 2016 at 15:17

Your Answer

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

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