rprojectionrgeo-shapefile

How to buffer/expand a census shapefile by 1 mile using rgeos::gBuffer?


I'm trying to expand/enlarge a shapefile of the Cherokee Nation that I downloaded from the Census dataset here: https://www2.census.gov/geo/tiger/TIGER2019/AIANNH/ using the rgeos package. The code I have for this currently is:

library(rgeos)
library(dplyr)

tribe_shp <- readOGR("/location of file/", layer = "tl_2019_us_aiannh")
tribe_shp <- tribe_shp %>% filter(GEOID == "5550R") #filter to cherokee nation

expand_obj <- gBuffer(tribe_shp, byid = F, width = 1000)

>Warning message:
In gBuffer(tribe_shp, byid = F, width = 1000) :
  Spatial object is not projected; GEOS expects planar coordinates

plot(expand_obj)

The resulting object loses the dataframe of the original SPDF and is pretty much just a circle that doesn't resemble the original shape at all. Is there something I am missing?


Solution

  • gBuffer uses the units of the data. In this case, the data are in degrees lat-long, and so the buffer is 1000 degrees in width. To make a buffer in metres transform to another coordinate system in metres.

    There are lots of coordinate systems, and you should really find one appropriate for your location. I think the US has a number of systems designed on a per-state basis, so that would probably be best. But for now I'll use EPSG:3857 which is what Google use for maps and isn't really that accurate.

    Read the data:

    tribe_shp <- readOGR("./", layer = "tl_2019_us_aiannh")
    

    Subset using selection - dplyr::filter doesn't work for me here but this will:

    tribe_shp = tribe_shp[tribe_shp$GEOID=="5550R",]
    

    Now transform to another coordinate system:

    tribe_shp_trans = spTransform(tribe_shp, "+init=epsg:3857")
    

    and make a 1km buffer. If you want a 1 mile buffer then use however many metres in a mile - one thousand six hundred and something?

    tribe_shp_buf = gBuffer(tribe_shp_trans, width=1000)
    

    If you plot them you can just about see the buffer is bigger than the original region:

    plot(tribe_shp_trans)
    plot(tribe_shp_buf,lty=3,add=TRUE)
    

    Detail of plot:

    enter image description here

    If you need the buffer in lat-long then do another transform of the buffer to "+init=epsg:4326".

    You can also do this with the more modern sf package, use st_read to read the data, st_transform to reproject, and st_buffer to do the buffer. It should even be faster.