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.

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)

2 Answers 2

2

Did you read the warning messages generated by the code you say didnt work?

> test_coords <- st_geometry(breweries[1:2,])
> buff_test_coords <- st_buffer(test_coords, dist = 10000)
dist is assumed to be in decimal degrees (arc_degrees).
Warning message:
In st_buffer.sfc(test_coords, dist = 10000) :
  st_buffer does not correctly buffer longitude/latitude data

Specifically dist is assumed to be in decimal degrees (arc_degrees).. You are doing a buffer that is 10000 degrees (and a sphere has only 360 all the way round.

The reason CRS 3488 works is because it is a system with metres as its coordinates, not degrees lat-long like breweries has:

> st_crs(sf_test_coords)
Coordinate Reference System:
  EPSG: 3488 
  proj4string: "+proj=aea +lat_1=34 +lat_2=40.5
                +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000
                +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 
                +units=m +no_defs"

I don't know where you picked 3488 out of, but because of the problems of flattening a sphere onto a flat map, you need to make sure that the coordinate system you choose for a buffer operation is appropriate to the geographic location of your data.

0

Copied from https://stackoverflow.com/a/60899504/7877917

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)

Your Answer

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

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