1

I would like to calculate a buffer of 5km around country polygons. However, I somehow always get a warning and an error. Here is a reproducible example code.

library(sp)
library(sf)
library(tmap) 
library(rnaturalearth)
countries_sp <- ne_countries(scale = 10)

countries_sp <- countries_sp[(countries_sp@data$geounit %in% "France"), ]

proj4string(countries_sp) <- CRS("+init=epsg:4326") 

countries_for_buff <- st_as_sf(countries_sp, crs ="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

countries_buff_5km <- st_buffer( countries_for_buff, 5000)

here the warning:

dist is assumed to be in decimal degrees (arc_degrees).
Warning message:
In st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle,  :
st_buffer does not correctly buffer longitude/latitude data

Then

tm_shape(countries_sp) +
  tm_polygons() +
tm_shape(countries_buff_5km) +
  tm_polygons()

gives the error

 Error in if (all(tiny) && show.messages) warning("all relative area sizes are below thres.poly",  : 
   missing value where TRUE/FALSE needed

What am I doing wrong? I guess there is something about the projection that does not quite work out the way it should.

  • Which query language? – jlSta Feb 8 at 11:51
  • The language is R, sorry – S Front Feb 8 at 11:53
1

This is unnecessary because the object is already in that coordinate system:

proj4string(countries_sp) <- CRS("+init=epsg:4326") 

This doesn't do a coordinate transformation - the crs argument is ignored - what made you think it did something?

countries_for_buff <- st_as_sf(countries_sp, crs ="+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

If you don't believe its ignored, try:

> countries_for_buff <- st_as_sf(countries_sp, crs ="bogus")
> 

So then when you do:

countries_buff_5km <- st_buffer( countries_for_buff, 5000)

then countries_for_buff is still in lat-long coordinates.

To convert from one to another you need either spTransform for sp objects or st_transform for sf objects.

so something like this on the sp object:

  csp2 = spTransform(countries_sp,"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")
  • great thank you! Because the plotting st_buffer does not accept sp objects, I used st_as_sf and the warning let me believe that something about the projection was wrong so I just simply tried some other projections – S Front Feb 8 at 12:07
  • If you want to buffer sp objects then use gBuffer from the rgeos package. But sf objects are mostly better and I'd advise converting to sf as early as possible. – Spacedman Feb 8 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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