Círculos de 1 km em torno de pontos longos em muitos lugares do mundo

11

Eu tenho centenas de pontos de latitude e longitude espalhados pelo mundo e preciso criar polígonos circulares em torno de cada um deles, com raio exatamente de 1000 metros. Entendo que os pontos devem primeiro ser projetados de graus (lat long) para algo com unidades de medidor, mas como isso pode ser feito sem procurar e definir manualmente as zonas UTM para cada ponto?

Aqui está um grande destaque para o primeiro ponto na Finlândia.

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

Mas com o segundo ponto (Canadá) não funciona (porque a zona UTM está errada).

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

Como isso pode ser feito sem obter e especificar manualmente o ponto da zona UTM por ponto? Não tenho mais informações por ponto do que o lat.

Atualizar:

Usando e combinando as ótimas respostas de AndreJ e Mike T, aqui está o código para versões e plotagens. Eles são diferentes na quarta casa decimal, aproximadamente, mas as duas respostas são muito boas!

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

(Não sabe ao certo como inserir o enredo na atualização).

Chris
fonte
1
Uma possível solução para a parte de pesquisa manual: E se você obtiver uma grade da Zona UTM e cruzar com seus pontos, para adicionar a zona apropriada como um atributo? O atributo pode ser o nome da zona ou o código EPSG, mas algo que pode ser alimentado como uma variável para selecionar automaticamente o CRS correto para cada ponto.
Chris W
1
Eu tenho um problema com "exatamente 1000m" e a frase "polígonos de círculo". Seus polígonos circulares precisam de segmentos infinitos com exatamente 1000m e a conversão para UTM (ou qualquer outro sistema planar) introduzirá ainda mais erros. Tenha cuidado com o uso de "exato".
Spacedman
Sim, eu não deveria ter expressado de maneira diferente. Eu quis dizer que 1100m ou 900m estaria muito desligado e que cerca de 20 segmentos no círculo estão ok.
22414 Chris

Respostas:

9

Semelhante ao @AndreJ, mas use uma projeção gnômica dinâmica , quero dizer uma projeção equidistante azimutal dinâmica para obter ainda mais precisão. Uma projeção AEQ centralizada em cada ponto projetará distâncias iguais em todas as direções, como um círculo em buffer. (Uma projeção de Mercator terá algumas distorções nas direções norte e leste, uma vez que envolve a lateral de um cilindro.)

Portanto, para o seu primeiro ponto na Finlândia, a string PROJ.4 será assim :

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

Então você pode criar uma função R para fazer esta projeção dinâmica:

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

Em seguida, faça algo parecido com isto no Canadá (item 2):

aeqd.buffer(the.points.sp[2,], 1000)
Mike T
fonte
1
Na página da wikipedia: "Nenhuma distorção ocorre no ponto tangente, mas a distorção aumenta rapidamente para longe dele". Você fez um cálculo de deslocamento de amostra? Talvez en.wikipedia.org/wiki/Azimuthal_equidistant_projection seja mais adequado.
Andrej
2
Qualquer projeção que tenha a escala correta na origem do círculo e seja conforme, será bem, simplesmente porque 1000m é muito pequeno. Para raios muito maiores, no entanto, uma projeção gnomônica será terrível. Você provavelmente pretendia estipular uma projeção Equidistante .
whuber
2
Ótimo feedback, uma projeção AEQ está obviamente tendo um desempenho muito melhor para essa técnica, então mudei a gnômica. O AEQP também aguenta distâncias muito maiores, como na faixa de 10.000 km.
Mike T
1
Talvez eu esteja entendendo mal o código, mas você só precisa construir o polígono do buffer uma vez, em qualquer projeção AEQD (o centro é sempre zero, o mínimo de coordenação é sempre -1k, o máximo é sempre + 1k. AEQD centrada em cada um dos pontos que você precisa para obter os valores de lat / lon ...
mkennedy
2
@mkennedy você tem um bom argumento. projectedé sempre sempre em (0, 0) e bufferedtem os pontos ± 1000 m nas direções x e y . Se for crítico otimizar isso, basta transformar uma versão cartesiana simples do buffer do dinâmico AEQD para WGS84.
Mike T
4

Em vez de procurar a zona UTM correta, você pode criar uma projeção mercator transversal personalizada para cada ponto com

+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

Desenhe o círculo nessa projeção. As coordenadas projetadas do vértice do círculo sempre serão as mesmas, portanto, você deve criá-las apenas uma vez. Para o seguinte, basta atribuir o novo CRS personalizado a eles.

Reprojete o círculo no EPSG: 4326 para uso posterior.

Dentro do intervalo de 1000m, o círculo será quase exato. Caso contrário (ou para círculos maiores), use em aeqdvez de tmerc.

AndreJ
fonte
0

E se você adotar a abordagem de criar um medidor de 1000 metros no EPSG: 4326 em torno de cada um dos seus pontos. Em seguida, converta o EPSG: 4326 no outro sistema de coordenadas? A vantagem de projetar o argumento é que você não precisa se preocupar com a curvatura da terra com o EPSG: 4326.

Greg
fonte
1
Como exatamente você criaria buffers de 1000 m a partir do EPSG: 4326, que possui unidades de comprimento em graus?
Mike T
Uma maneira de abordar isso é criar um buffer de 1000 metros no EPSG: 32635. Converta isso em EPSG: 4326 e agora você terá o número que precisa.
Greg
1
Essa é a mesma abordagem descrita na pergunta, juntamente com as limitações dessa técnica.
Mike T