1 km okrąża długie punkty w wielu miejscach na świecie

Mam setki lat-długich punktów rozmieszczonych na całym świecie i muszę tworzyć wielokąty kołowe wokół każdego z nich, o promieniu dokładnie 1000 metrów. Rozumiem, że punkty należy najpierw rzutować ze stopni (długości lat) na coś z jednostkami metrowymi, ale jak można to zrobić bez ręcznego wyszukiwania i definiowania stref UTM dla każdego punktu?

Oto mwe za pierwszy punkt w Finlandii.

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

Ale z drugim punktem (Kanada) to nie działa (ponieważ niewłaściwa strefa UTM).

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

Jak można to zrobić bez ręcznego uzyskiwania i określania punktu strefy UTM na punkt? Nie mam więcej informacji na punkt niż lat long.

Aktualizacja:

Używając i łącząc świetne odpowiedzi zarówno AndreJ, jak i Mike T, oto kod dla obu wersji i fabuł. Różnią się mniej więcej na czwartym miejscu po przecinku, ale obie bardzo dobre odpowiedzi!

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

(Nie jestem pewien, jak wprowadzić fabułę do aktualizacji).

Podobnie jak @AndreJ, ale używaj dynamicznej projekcji gnomicznej , mam na myśli dynamiczną, azymutalną równoodległą projekcję dla jeszcze większej dokładności. Projekcja AEQ wyśrodkowana na każdym punkcie spowoduje rzutowanie równych odległości we wszystkich kierunkach, takich jak buforowany okrąg. (Projekcja Mercatora będzie miała pewne zniekształcenia w kierunku północnym i wschodnim, ponieważ owija się wokół boku cylindra.)

Tak więc dla twojego pierwszego punktu w Finlandii łańcuch PROJ.4 będzie wyglądał następująco:

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

Możesz więc wykonać funkcję R, aby wykonać tę dynamiczną projekcję:

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

Następnie zrób coś takiego dla Kanady (pozycja 2):

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

Zamiast szukać właściwej strefy UTM, możesz utworzyć niestandardową projekcję poprzecznego mercatora dla każdego punktu, za pomocą którego

+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

Narysuj okrąg w tej projekcji. Współrzędne rzutowanego wierzchołka okręgu zawsze będą takie same, więc musisz je utworzyć tylko raz. W tym celu wystarczy przypisać im nowy niestandardowy CRS.

Ponownie skieruj okrąg na EPSG: 4326 do dalszego użycia.

W zakresie 1000 m okrąg będzie prawie dokładny. Jeśli nie (lub w przypadku większych kręgów), użyj aeqdzamiast tmerc.

Co się stanie, jeśli podejmiesz podejście polegające na utworzeniu 1000 metrów w EPSG: 4326 wokół każdego z punktów? Następnie przekonwertować EPSG: 4326 na inny układ współrzędnych? Zaletą rzutowania punktu jest to, że dzięki EPSG: 4326 nie musisz się martwić o krzywiznę ziemi.