11

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).

  • 1
    A possible solution to the manual searching part: What if you get a UTM Zone grid and intersect that with your points, so that you add the appropriate zone as an attribute? The attribute could either be zone name or EPSG code, but something that could be fed in as a variable to automatically select the right CRS for each point. – Chris W Nov 9 '14 at 20:32
  • 1
    I have a problem with "exactly 1000m" and the phrase "circle-polygons". Your circle-polygons need infinite segments to be exactly 1000m, and converting to UTM (or any other planar system) is going to introduce even more errors. Be careful with the use of "exact". – Spacedman Nov 9 '14 at 23:23
  • Yes, I shouldn't have expressed it differently. I meant that 1100m or 900m would be too off, and that about 20 segments on the circle is ok. – Chris Nov 10 '14 at 14:49
9

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)
  • 1
    From the wikipedia page: "No distortion occurs at the tangent point, but distortion increases rapidly away from it". Have you made a sample offset computation? Maybe en.wikipedia.org/wiki/Azimuthal_equidistant_projection is better suited. – AndreJ Nov 10 '14 at 12:01
  • 2
    Any projection that has the correct scale at the circle's origin and is conformal there will do fine, simply because 1000m is so small. For much larger radii, however, a Gnomonic projection will be awful. You probably meant to stipulate an Equidistant projection. – whuber Nov 10 '14 at 14:53
  • 2
    Great feedback, an AEQ projection is obviously performing much better for this technique, so I've switched out gnomic. AEQP will also hold up for much greater distances, like in the 10,000+ km range. – Mike T Nov 10 '14 at 16:58
  • 1
    I may be misunderstanding the code, but you only need to build the buffer polygon once, in any AEQD projection (Center is always zero, min coord is always -1k, max is always +1k. Then unproject it to lat/lon using an AEQD centered on each of the points that you need to get the lat/lon values... – mkennedy Nov 10 '14 at 21:01
  • 2
    @mkennedy you have a good point. projected is indeed always at (0, 0), and buffered has the points ±1000 m in x- and y-directions. If it were critical to optimise this, then just transform a simple Cartesian version of the buffer from the dynamic AEQD to WGS84. – Mike T Nov 10 '14 at 21:17
4

Instead of searching for the right UTM zone, you could create a custom transverse mercator projection for every point with

+proj=tmerc +lat_0=.... +lon_0=... +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Draw the circle in that projection. The projected circle vertex coordinates will always be the same, so you have to create them only once. For the following, just assign the new custom CRS to them.

Reproject the circle to EPSG:4326 for further use.

Within the range of 1000m, the circle will be almost exact. If not (or for larger circles), use aeqd instead of tmerc.

0

What if you take the approach of creating a 1000 meter in EPSG:4326 around each of your points. Then convert the EPSG:4326 to your other coordinate system? The advantage of projecting the point, is that you don't have to worry about curvature of the earth with EPSG:4326.

  • 1
    How exactly would you create 1000 m buffers from EPSG:4326, which has length units in degrees? – Mike T Nov 10 '14 at 1:54
  • One way I'd approach this is to create a 1000 meter buffer in EPSG:32635. Convert that to EPSG:4326 and now you would have the number you need. – Greg Nov 10 '14 at 3:50
  • 1
    That's the same approach described in the question, along with the limitations of this technique. – Mike T Nov 10 '14 at 11:28

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.