1

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)
user2438292
  • 87
  • 9

1 Answers1

1

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)

dieghernan
  • 2,690
  • 8
  • 16
  • Thanks for feedback. I should clarify...I'm interested in drawing the buffers arbitrary distances from the points (e.g. 10 km, 100 km). Seems like my general approach (i.e. first transform to CRS that uses distance in meters, then plot) is the approved solution based on other posts. Also, seems like there are many known issues related to performing this calculation for many arbitrary points on the globe. – user2438292 Mar 30 '20 at 18:16
  • You’re welcome. If you need to buffer in units then you have no choice, you have to transform. Meters over lonlat doesn’t work, but if you really need to have the final result on lonlat you can try `project -> buffer -> unproject (4326)` – dieghernan Mar 30 '20 at 20:43