R – Create 1km Circles Around Latitude-Longitude Points Globally

coordinate systemr

I have hundreds of lat-long-points spread around the world, and must create circle-polygons around each of them, with radius exactly 1000 meters. I understand the points must first be projected from degrees (lat long) to something with meter-units, but how can this be done without manually searching and defining the UTM-zones for each point?

Here is a mwe for first point in Finland.

library(sp)
library(rgdal)
library(rgeos)
the.points.latlong <- data.frame(
  Country=c("Finland", "Canada", "Tanzania", "Bolivia", "France"),
  lat=c(63.293001, 54.239631, -2.855123, -13.795272, 48.603949),
  long=c(27.472918, -90.476303, 34.679950, -65.691146, 4.533465))
the.points.sp <- SpatialPointsDataFrame(the.points.latlong[, c("long", "lat")], data.frame(ID=seq(1:nrow(the.points.latlong))), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

the.points.projected <- spTransform(the.points.sp[1, ], CRS( "+init=epsg:32635" ))  # Only first point (Finland)
the.circles.projected <- gBuffer(the.points.projected, width=1000, byid=TRUE)
plot(the.circles.projected)
points(the.points.projected)

the.circles.sp <- spTransform(the.circles.projected, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

But with second point (Canada) it doesn't work (because wrong UTM-zone).

the.points.projected <- spTransform(the.points.sp[2, ], CRS( "+init=epsg:32635" ))

How can this be done without manually getting and specifying UTM-zone point per point? I don't have any more info per point than lat long.

Update:

Using and combining the great answers from both AndreJ and Mike T, here is the code for both versions and plots. They are different on 4th decimal or so, but both very good answers!

gnomic.buffer <- function(p, r) {
  stopifnot(length(p) == 1)
  gnom <- sprintf("+proj=gnom +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
                  p@coords[[2]], p@coords[[1]])
  projected <- spTransform(p, CRS(gnom))
  buffered <- gBuffer(projected, width=r, byid=TRUE)
  spTransform(buffered, p@proj4string)
}

custom.buffer <- function(p, r) {
  stopifnot(length(p) == 1)
  cust <- sprintf("+proj=tmerc +lat_0=%s +lon_0=%s +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs", 
                  p@coords[[2]], p@coords[[1]])
  projected <- spTransform(p, CRS(cust))
  buffered <- gBuffer(projected, width=r, byid=TRUE)
  spTransform(buffered, p@proj4string)
}

test.1 <- gnomic.buffer(the.points.sp[2,], 1000)
test.2 <- custom.buffer(the.points.sp[2,], 1000)

library(ggplot2)
test.1.f <- fortify(test.1)
test.2.f <- fortify(test.2)
test.1.f$transf <- "gnomic"
test.2.f$transf <- "custom"
test.3.f <- rbind(test.1.f, test.2.f)

p <- ggplot(test.3.f, aes(x=long, y=lat, group=transf))
p <- p + geom_path()
p <- p + facet_wrap(~transf)
p

(Not sure how to get the plot into the update).

Best Answer

Similar to @AndreJ, but use a dynamic gnomic projection, I mean a dynamic azimuthal equidistant projection for even more accuracy. An AEQ projection centred on each point will project equal distances in all directions, such as a buffered circle. (A Mercator projection will have some distortions in north and eastern directions, since it wraps around the side of a cylinder.)

So for your first point around Finland, the PROJ.4 string will look like this:

+proj=aeqd +lat_0=63.293001 +lon_0=27.472918 +x_0=0 +y_0=0

So you can make an R function to make this dynamic projection:

aeqd.buffer <- function(p, r)
{
    stopifnot(length(p) == 1)
    aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
                    p@coords[[2]], p@coords[[1]])
    projected <- spTransform(p, CRS(aeqd))
    buffered <- gBuffer(projected, width=r, byid=TRUE)
    spTransform(buffered, p@proj4string)
}

Then do something like this for Canada (item 2):

aeqd.buffer(the.points.sp[2,], 1000)
Related Question