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

Improve this question
3
  • 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
    Commented Nov 9, 2014 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
    Commented Nov 9, 2014 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
    Commented Nov 10, 2014 at 14:49

3 Answers 3

Reset to default
8

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)
Improve this answer
5
  • 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
    Commented Nov 10, 2014 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
    Commented Nov 10, 2014 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
    Commented Nov 10, 2014 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
    Commented Nov 10, 2014 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
    Commented Nov 10, 2014 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.

Improve this answer
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.

Improve this answer
3
  • 1
    How exactly would you create 1000 m buffers from EPSG:4326, which has length units in degrees?
    – Mike T
    Commented Nov 10, 2014 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
    Commented Nov 10, 2014 at 3:50
  • 1
    That's the same approach described in the question, along with the limitations of this technique.
    – Mike T
    Commented Nov 10, 2014 at 11:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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