This is probably one of the simplest GIS tasks to do, but I am struggling. I have tried to adapt this answer but it appears not to work Create buffers around points and merge with SpatialPolygonsDataFrame to create a list of matched polygons and the area of intersection

Here is my data with long and lat in WGS84

S10.df<-structure(list(id = c(273, 1028, 3635, 6693, 13430), long_wgs = c(-1.564636192, 
-1.56677146, -1.554945908, -1.56773222, -1.564188582), lat_wgs = c(53.82838137, 
53.82924493, 53.83130225, 53.82848279, 53.82857442)), class = "data.frame", row.names = c("273", 
"1023", "3357", "6326", "12301"))

and this code:

stores.sf <- st_as_sf(S10.df, coords = c("long_wgs", "lat_wgs"), crs = 27700, dim = "XY")
stores.buffer <- st_buffer(stores.sf, 500)

to create 500m polygon buffers around each location. When I try to plot I just get one large circle. Further inspection shows:

> stores.buffer$geometry
Geometry set for 5 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -501.5677 ymin: -446.1716 xmax: 498.4451 ymax: 553.8313
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
POLYGON ((498.4354 53.82838, 497.7501 27.6604, ...
POLYGON ((498.4332 53.82924, 497.748 27.66127, ...
POLYGON ((498.4451 53.8313, 497.7598 27.66332, ...
POLYGON ((498.4323 53.82848, 497.747 27.6605, 4...
POLYGON ((498.4358 53.82857, 497.7506 27.6606, ...

Which looks an odd set of POLYGON co-ordinates. Has the "transformation" worked so that st_buffer knows the CRS 27700 is in units of meters? Ideally I'd like the "answer" back in WGS84.

  • You are defining your buffer distance assuming a planar projection. Since your data is unprojected the function is probably assuming a great circle distance, which would translate to 500 kilometers. Alternately, it could be using the projection units which would be 500 degrees, not distance units at all. – Jeffrey Evans Oct 1 at 18:41
up vote 0 down vote accepted

The problem is that the data the you load is in WGS84 coordinates, however when you use the function st_as_sf() you specify the CRS as EPSG:27700, thus the buffer has strange coordinates.

You should first load your data with argument crs = 4326 then transform them to crs = 27700:

library("sf")

S10.df <- structure(list(id = c(273, 1028, 3635, 6693, 13430),
                         long_wgs = c(-1.564636192, -1.56677146, -1.554945908, -1.56773222, -1.564188582),
                         lat_wgs = c(53.82838137, 53.82924493, 53.83130225, 53.82848279, 53.82857442)),
                    class = "data.frame", row.names = c("273", "1023", "3357", "6326", "12301"))

stores_wgs <- st_as_sf(S10.df, coords = c("long_wgs", "lat_wgs"),
                       crs = 4326, dim = "XY")
stores_proj <- st_transform(stores_wgs, 27700)

st_buffer(stores_proj, 500)
#> Simple feature collection with 5 features and 1 field
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 428049.6 ymin: 436489.6 xmax: 429889.3 ymax: 437818.5
#> epsg (SRID):    27700
#> proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
#>          id                       geometry
#> 273     273 POLYGON ((429253.5 436989.6...
#> 1023   1028 POLYGON ((429112.4 437084.8...
#> 3357   3635 POLYGON ((429889.3 437318.5...
#> 6326   6693 POLYGON ((429049.6 436999.7...
#> 12301 13430 POLYGON ((429282.8 437011.3...

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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