12 votos

1km círculos alrededor de lat-long-puntos de que en muchos lugares del mundo

Tengo cientos de lat-long-puntos repartidos por todo el mundo, y debe crear círculo-polígonos alrededor de cada uno de ellos, con radio exactamente 1000 metros. Entiendo los puntos primero debe ser proyectada de grados (lat long) a algo con medidor de unidades, pero, ¿cómo se puede hacer esto sin la búsqueda manual y la definición de la UTM-zonas para cada punto?

Aquí está una mwe para el primer punto en Finlandia.

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

Pero con el segundo punto (Canadá) no funciona (ya que un mal UTM-zona de).

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

Cómo se puede hacer esto sin manualmente llegar y especificando UTM-zona de punto por punto? No tengo más info por punto de lat long.

Actualización:

El uso y la combinación de las grandes respuestas de ambos AndreJ y Mike T, aquí está el código para ambas versiones y parcelas. Son diferentes en la 4ª decimal o así, pero ambas muy buenas respuestas!

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

(No estoy seguro de cómo obtener la parcela en la actualización).

9voto

hernan43 Puntos 566

Similares a @AndreJ, pero el uso de una dinámica gnomónicos proyección, me refiero a una dinámica azimutal equidistante de proyección para obtener aún más la precisión. Un AEQ proyección centrada en cada punto del proyecto distancias iguales en todas las direcciones, como un tampón de círculo. (Una proyección de Mercator tendrá algunas distorsiones en el norte y este de las direcciones, ya que se envuelve alrededor de la parte lateral de un cilindro.)

Así que para el primer punto alrededor de Finlandia, el PROJ.4 cadena de aspecto como este:

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

Así que usted puede hacer una función de R para hacer esta dinámica de proyección:

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

Luego de hacer algo como esto para Canadá (punto 2):

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

4voto

SpliFF Puntos 214

En lugar de buscar el derecho UTM zona, se podría crear una personalizada transversal de mercator la proyección de cada punto de

+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

Dibujar el círculo en que la proyección. La proyectada círculo de vértices de coordenadas siempre será el mismo, por lo que usted tiene que crear ellos sólo una vez. Para la siguiente, simplemente asignar a la nueva costumbre de CRS para ellos.

Reproyectar el círculo a EPSG:4326 para su uso posterior.

Dentro de la gama de 1000m, el círculo será casi exacta. Si no (o para círculos más grandes), el uso de aeqd en lugar de tmerc.

0voto

ane Puntos 116

¿Qué pasa si usted toma el enfoque de la creación de un 1000 metros en EPSG:4326 alrededor de cada uno de sus puntos. A continuación, convertir el EPSG:4326 a su otro sistema de coordenadas? La ventaja de proyectar el punto, es que usted no tiene que preocuparse acerca de la curvatura de la tierra con EPSG:4326.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X