2

The use case I have is this: I want to take points (in my case a data.frame of lat/long coordinates) and calculate a polygon around each point of a given radius (in km, miles, whatever). To be clear, not a polygon around all points at once, but a polygon around each point separately. This seems like it would be straightforward, but how does one do this in R?

The next thing I want to do is take those polygons and get well-known text (WKT) polygons because the web API I am querying against wants WKT. I think I can get a WKT polygon using rgeos like

library(rgeos)
g3 <- readWKT("POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))")
writeWKT(g3)

[1] "POLYGON ((1.0000000000000000 1.0000000000000000, 1.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 5.0000000000000000 1.0000000000000000, 1.0000000000000000 1.0000000000000000), (2.0000000000000000 2.0000000000000000, 3.0000000000000000 2.0000000000000000, 3.0000000000000000 3.0000000000000000, 2.0000000000000000 3.0000000000000000, 2.0000000000000000 2.0000000000000000))"
3

Start with a data frame of lat, lon, and some ID variable. Convert to a SpatialPointsDataFrame, and give it the usual lat-long coordinate system code:

require(sp)
require(rgeos)
d=data.frame(lat=c(33.95,34.95,34.70), lon=c(-118.40,-118.22,-118.43),ID=1:3)
coordinates(d)=~lon+lat
proj4string(d)=CRS("+init=epsg:4326")

Now do the buffering:

 buf = gBuffer(d,width=0.2)
 plot(buf)

Extract the parts and convert to WKT:

writeWKT(buf[1])
[1] "POLYGON ((-117.4000000000000057 33.9500000000000028, -117.4489434799999970 33.6409830099999994, -117.5909830100000022 33.3622147499999997, -117.8122147499999954 33.1409830099999994, .....

Note this has done the buffering for all three points in the data frame, but you can only do that if you want the same buffer size for each point. I think the rgeos guys might be working on a new version where the buffer size can be a vector, for example another column in the data frame, and so different for each point.

I'll repeat the warning about working with lat-long as if it was cartesian - use spTransform to convert to metres using a locally cartesian coordinate system such as a UTM zone. If you really need to do this on great circle distances then rgeos, and GEOS in particular, probably won't help you.

  • I just started playing with gBuffer. I am running your code, but I receive the following error message. Error: is.logical(byid) is not TRUE. Would you mind if I ask you what causes this error? Thank you for your help. – jazzurro Dec 18 '14 at 7:52
  • Create a reproducible example and post it as a new question. – Spacedman Dec 18 '14 at 9:47
  • Thank you for your reply. I am happy to create an example. But, your code above does not work for me although I uploaded sp and rgeos. Would you mind if I ask you to run the code on your machine and see if it works or not? – jazzurro Dec 18 '14 at 13:51
  • 1
    agh! they changed the args for gBuffer, you now have to spell it out unless you also want byid=TRUE – Spacedman Dec 18 '14 at 14:03
0

Okay, this does what I want:

library(rgeos)
foo <- function(lat, lon, ...)
{
  latlon <- paste(lon, lat)
  poly <- readWKT(sprintf("POLYGON((%s,%s,%s,%s))", latlon,latlon,latlon,latlon))
  g <- gBuffer(poly, ...)
  writeWKT(g)
}

foo(lat=33.95, lon=-118.40, width=0.2)

[1] "POLYGON ((-118.2000000000000028 33.9500000000000028, -118.2097887000000043 33.8881966000000006, -118.2381965999999949 33.8324429500000008, -118.2824429500000036 33.7881965999999991, -118.3381966000000034 33.7597887000000014, -118.4000000000000057 33.7500000000000000, -118.4618033999999938 33.7597887000000014, -118.5175570499999935 33.7881965999999991, -118.5618034000000023 33.8324429500000008, -118.5902113000000071 33.8881966000000006, -118.5999999999999943 33.9500000000000028, -118.5902113000000071 34.0118033999999980, -118.5618034000000023 34.0675570499999978, -118.5175570499999935 34.1118033999999994, -118.4618033999999938 34.1402112999999972, -118.4000000000000057 34.1499999999999986, -118.3381966000000034 34.1402112999999972, -118.2824429500000036 34.1118033999999994, -118.2381965999999949 34.0675570499999978, -118.2097887000000043 34.0118033999999980, -118.2000000000000028 33.9500000000000028))"

It's sort of cheating I guess since I'm creating polygon with all four corners the same exact point, then use gBuffer to make a wider polygon of width width, then write a WKT polygon

Thanks to Simon Goring for the help (https://twitter.com/sjGoring/status/408458101252030464)

  • Can't you use gBuffer on SpatialPoint objects to get circular polygon buffers? Also, you need to be very careful with mixing lat-long coordinates with anything you think is a 'distance'. For a small region, convert to a cartesian coordinate reference system first. – Spacedman Dec 5 '13 at 8:18
  • Sorry for the late reply on this, Thanks for pointing that out – sckott Mar 11 '14 at 3:16

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.