Why is st_buffer function not creating an R object that correctly displays in mapview?

1.4k Views Asked by At

Trying to get st_buffer function to display a buffer within mapview. Got it to work, but I had to first perform a transform (I doubt this is necessary). Looking for a more straightforward way to do this.

EDIT: To clarify...the buffers (in this case, polygons) I'd like to draw around the points would be a distance (say, kilometers) around the points.

library(sf)
library(mapview)

data("breweries")

test_coords <- st_geometry(breweries[1:2,])

# This code doesn't work. Not sure why.
# buff_test_coords <- st_buffer(test_coords, dist = 10000)
# mapview(test_coords) + mapview(buff_test_coords)

# This code words. Not sure what's special about transforming to 3488
sf_test_coords <- test_coords %>% st_transform(3488)
sf_buff_test_coords <- st_buffer(sf_test_coords, 10000)
sf_buff_test_coords2 <- st_transform(sf_buff_test_coords, 4326)
mapview(test_coords) + mapview(sf_buff_test_coords2)
1

There are 1 best solutions below

2
On

Well, the warning is pretty clear, buffering doesn't work well for non-projected data.

#> Warning in st_buffer.sfc(test_coords, dist = 2): st_buffer does not correctly
#> buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).

Your data shows proj4string: "+proj=longlat +datum=WGS84 +no_defs" so you either project it or change the approach. On top of that, you tried to project lonlat (-180,180 degrees) by 10000, meaning 10000 degrees. So that is a non-sense buffer, buffering works on the same units of the projection.

You have here two approaches without projecting:

  • Buffering by 2 degrees, it works but the buffer is quite strange.
  • Another option is just to plot two times the same point, but not passing it to POLYGON, that is what the buffer does. Just plot it the second time with a bigger cex.
library(sf)
#> Warning: package 'sf' was built under R version 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(mapview)
#> Warning: package 'mapview' was built under R version 3.5.3

data("breweries")

test_coords <- st_geometry(breweries[1:2,])
st_crs(test_coords)
#> Coordinate Reference System:
#>   EPSG: 4326 
#>   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
buff_test_coords <- st_buffer(test_coords, dist = 2)
#> Warning in st_buffer.sfc(test_coords, dist = 2): st_buffer does not correctly
#> buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).

#Buffering non-projected coords
mapview(test_coords) + mapview(buff_test_coords)


#Plotting the same points with a bigger cex
mapview(test_coords) + mapview(test_coords, cex=100, col.regions ="red")

Created on 2020-03-28 by the reprex package (v0.3.0)