1km lingkaran di sekitar titik-titik lat-panjang di banyak tempat di dunia

11

Saya memiliki ratusan titik lat-panjang yang tersebar di seluruh dunia, dan harus membuat lingkaran-poligon di sekitar masing-masing, dengan radius tepat 1.000 meter. Saya mengerti bahwa poin pertama harus diproyeksikan dari derajat (panjang lat) ke sesuatu dengan meter-unit, tetapi bagaimana hal ini dapat dilakukan tanpa secara manual mencari dan menentukan zona UTM untuk setiap titik?

Ini adalah mwe untuk poin pertama di 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"))

Tetapi dengan poin kedua (Kanada) itu tidak berfungsi (karena salah zona UTM).

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

Bagaimana ini bisa dilakukan tanpa secara manual mendapatkan dan menentukan titik zona-UTM per titik? Saya tidak punya info lebih banyak per poin dari lat panjang.

Memperbarui:

Menggunakan dan menggabungkan jawaban-jawaban hebat dari AndreJ dan Mike T, berikut adalah kode untuk kedua versi dan plot. Mereka berbeda pada desimal ke-4, tetapi keduanya adalah jawaban yang sangat bagus!

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

(Tidak yakin bagaimana memasukkan plot ke dalam pembaruan).

Chris
sumber
1
Solusi yang mungkin untuk bagian pencarian manual: Bagaimana jika Anda mendapatkan kisi UTM Zone dan memotongnya dengan poin Anda, sehingga Anda menambahkan zona yang sesuai sebagai atribut? Atribut dapat berupa nama zona atau kode EPSG, tetapi sesuatu yang dapat dimasukkan sebagai variabel untuk secara otomatis memilih CRS yang tepat untuk setiap titik.
Chris W
1
Saya memiliki masalah dengan "tepat 1000 m" dan frasa "lingkaran-poligon". Lingkaran-poligon Anda membutuhkan segmen tak hingga persis 1000m, dan mengonversi ke UTM (atau sistem planar lainnya) akan menimbulkan lebih banyak kesalahan. Hati-hati dengan penggunaan "tepat".
Spacedman
Ya, saya seharusnya tidak mengekspresikannya secara berbeda. Maksud saya, 1100m atau 900m akan terlalu mati, dan sekitar 20 segmen pada lingkaran itu ok.
Chris

Jawaban:

9

Mirip dengan @AndreJ, tetapi menggunakan proyeksi gnomic yang dinamis , maksud saya proyeksi azimutal yang sama untuk lebih akurat. Proyeksi AEQ yang dipusatkan pada setiap titik akan memproyeksikan jarak yang sama di semua arah, seperti lingkaran buffered. (Proyeksi Mercator akan memiliki beberapa distorsi di arah utara dan timur, karena membungkus sisi silinder.)

Jadi untuk poin pertama Anda di Finlandia, string PROJ.4 akan terlihat seperti ini:

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

Jadi, Anda dapat membuat fungsi R untuk membuat proyeksi dinamis ini:

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

Kemudian lakukan sesuatu seperti ini untuk Kanada (item 2):

aeqd.buffer(the.points.sp[2,], 1000)
Mike T
sumber
1
Dari halaman wikipedia: "Tidak ada distorsi yang terjadi pada titik singgung, tetapi distorsi meningkat dengan cepat darinya". Sudahkah Anda membuat perhitungan sampel offset? Mungkin en.wikipedia.org/wiki/Azimuthal_equidistant_projection lebih cocok.
AndreJ
2
Setiap proyeksi yang memiliki skala yang tepat pada asal lingkaran dan sesuai di sana akan baik-baik saja, hanya karena 1000m sangat kecil. Namun, untuk jari-jari yang jauh lebih besar, proyeksi Gnomonic akan mengerikan. Anda mungkin bermaksud untuk menetapkan proyeksi yang Sama .
whuber
2
Umpan balik yang bagus, proyeksi AEQ jelas berkinerja jauh lebih baik untuk teknik ini, jadi saya sudah mengganti gnomic. AEQP juga akan bertahan untuk jarak yang jauh lebih besar, seperti di kisaran 10.000+ km.
Mike T
1
Saya mungkin salah paham kode, tetapi Anda hanya perlu membangun buffer poligon sekali, dalam proyeksi AEQD (Center selalu nol, coord min selalu -1k, maks selalu +1k. Kemudian tidak memproyeksikannya ke lat / lon menggunakan AEQD berpusat pada masing-masing poin yang Anda butuhkan untuk mendapatkan nilai lat / lon ...
mkennedy
2
@kennedy Anda punya poin bagus. projectedmemang selalu di (0, 0), dan bufferedmemiliki poin ± 1000 m dalam x - dan y -directions. Jika sangat penting untuk mengoptimalkan ini, maka cukup ubah buffer Cartesian versi sederhana dari AEQD dinamis ke WGS84.
Mike T
4

Alih-alih mencari zona UTM yang tepat, Anda dapat membuat proyeksi mercator melintang khusus untuk setiap titik dengan

+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

Gambarlah lingkaran dalam proyeksi itu. Koordinat simpul lingkaran yang diproyeksikan akan selalu sama, jadi Anda harus membuatnya hanya sekali. Untuk yang berikut, tetapkan CRS khusus baru untuk mereka.

Proyeksi ulang lingkaran ke EPSG: 4326 untuk penggunaan lebih lanjut.

Dalam kisaran 1000m, lingkaran akan hampir tepat. Jika tidak (atau untuk lingkaran yang lebih besar), gunakan aeqdsebagai ganti tmerc.

AndreJ
sumber
0

Bagaimana jika Anda mengambil pendekatan menciptakan 1000 meter dalam EPSG: 4326 di sekitar setiap poin Anda. Kemudian konversi EPSG: 4326 ke sistem koordinat Anda yang lain? Keuntungan dari memproyeksikan intinya, adalah Anda tidak perlu khawatir tentang kelengkungan bumi dengan EPSG: 4326.

Greg
sumber
1
Bagaimana tepatnya Anda membuat buffer 1000 m dari EPSG: 4326, yang memiliki satuan panjang dalam derajat?
Mike T
Salah satu cara saya mendekati ini adalah dengan membuat buffer 1000 meter di EPSG: 32635. Ubah itu menjadi EPSG: 4326 dan sekarang Anda akan memiliki nomor yang Anda butuhkan.
Greg
1
Itu pendekatan yang sama dijelaskan dalam pertanyaan, bersama dengan keterbatasan teknik ini.
Mike T