1 км кругов вокруг лат-длинных точек во многих местах мира

11

У меня есть сотни широчайших точек по всему миру, и я должен создать многоугольники вокруг каждого из них с радиусом ровно 1000 метров. Я понимаю, что точки должны сначала проецироваться из градусов (широта в длину) во что-то с единицами измерения, но как это можно сделать без ручного поиска и определения UTM-зон для каждой точки?

Вот mwe для первой точки в Финляндии.

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

Но со вторым пунктом (Канада) он не работает (из-за неправильной UTM-зоны).

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

Как это можно сделать без ручного получения и указания точки зоны UTM для каждой точки? У меня нет больше информации по каждому пункту, чем последний.

Обновить:

Используя и комбинируя отличные ответы от AndreJ и Mike T, вот код как для версий, так и для графиков. Они различаются на четвертом десятичном числе или около того, но оба очень хорошие ответы!

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

(Не уверен, как получить сюжет в обновлении).

Крис
источник
1
Возможное решение для поиска вручную: что если вы получите сетку UTM Zone и пересечете ее с вашими точками, чтобы добавить соответствующую зону в качестве атрибута? Атрибутом может быть либо имя зоны, либо код EPSG, но что-то, что может быть введено как переменная для автоматического выбора правильного CRS для каждой точки.
Крис W
1
У меня проблема с "ровно 1000м" и фразой "круг-полигоны". Вашим круговым полигонам нужны бесконечные отрезки, чтобы быть ровно 1000 м, и преобразование в UTM (или любую другую планарную систему) приведет к еще большему количеству ошибок. Будьте осторожны с использованием «точных».
Spacedman
Да, я не должен был выражать это по-другому. Я имел в виду, что 1100 м или 900 м будет слишком далеко, и что около 20 сегментов на круге в порядке.
Крис

Ответы:

9

Подобно @AndreJ, но использую динамическую проекцию гномов, я имею в виду динамическую азимутальную эквидистантную проекцию для еще большей точности. Проекция AEQ с центром в каждой точке будет проецировать одинаковые расстояния во всех направлениях, например, в буферную окружность. (Проекция Меркатора будет иметь некоторые искажения в северном и восточном направлениях, поскольку она оборачивается вокруг стороны цилиндра.)

Итак, для вашей первой точки вокруг Финляндии строка PROJ.4 будет выглядеть так :

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

Таким образом, вы можете сделать функцию R для создания этой динамической проекции:

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

Затем сделайте что-то подобное для Канады (пункт 2):

aeqd.buffer(the.points.sp[2,], 1000)
Майк Т
источник
1
Со страницы википедии: «Никаких искажений не возникает в точке касания, но искажение быстро увеличивается от нее». Вы сделали расчет смещения образца? Может быть, en.wikipedia.org/wiki/Azimuthal_equidistant_projection лучше подходит.
AndreJ
2
Любая проекция, которая имеет правильную шкалу в начале круга и является конформной там, подойдет, просто потому, что 1000 м очень мала. Однако для гораздо больших радиусов гномическая проекция будет ужасной. Вы, вероятно, хотели указать проекцию Equidistant .
whuber
2
Отличная обратная связь, AEQ-проекция, очевидно, работает намного лучше для этой техники, поэтому я выбрал gnomic. AEQP также выдержит гораздо большие расстояния, например, в диапазоне более 10 000 км.
Майк Т
1
Возможно, я неправильно понимаю код, но вам нужно построить многоугольник буфера только один раз, в любой проекции AEQD (центр всегда равен нулю, минимальная координата всегда равна -1k, максимальная всегда равна + 1k. Затем отмените проецирование на широту / долготу, используя AEQD сосредоточен на каждой из точек, которые вам нужны, чтобы получить значения
широты / долготы
2
@mkennedy у тебя есть хорошая мысль. projectedдействительно всегда в (0, 0), и bufferedимеет точки ± 1000 м в x- и y- направлениях. Если это было критически важно для оптимизации, то просто преобразуйте простую декартову версию буфера из динамического AEQD в WGS84.
Майк Т
4

Вместо поиска нужной UTM-зоны, вы можете создать собственную поперечную проекцию Меркатора для каждой точки с

+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

Нарисуйте круг в этой проекции. Координаты вершин проекции круга всегда будут одинаковыми, поэтому их нужно создавать только один раз. Для следующего, просто назначьте им новый пользовательский CRS.

Перепроектируйте круг в EPSG: 4326 для дальнейшего использования.

В пределах 1000 м круг будет почти точным. Если нет (или для больших кругов), используйте aeqdвместо tmerc.

Andrej
источник
0

Что, если вы выберете 1000 метров в EPSG: 4326 вокруг каждой из ваших точек. Затем преобразовать EPSG: 4326 в другую систему координат? Преимущество проецирования точки в том, что вам не нужно беспокоиться о искривлении земли с помощью EPSG: 4326.

Greg
источник
1
Как именно вы создадите 1000-метровые буферы из EPSG: 4326, у которого есть единицы измерения длины в градусах?
Майк Т
Один из подходов, который я бы использовал, - создать 1000-метровый буфер в EPSG: 32635. Переведите это в EPSG: 4326, и теперь у вас будет нужный номер.
Грег
1
Это тот же подход, который описан в вопросе, наряду с ограничениями этой техники.
Майк Т