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).
sumber
Jawaban:
Mirip dengan @AndreJ, tetapi menggunakan
proyeksi gnomic yangdinamis , maksud sayaproyeksiazimutal 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:
Jadi, Anda dapat membuat fungsi R untuk membuat proyeksi dinamis ini:
Kemudian lakukan sesuatu seperti ini untuk Kanada (item 2):
sumber
projected
memang selalu di (0, 0), danbuffered
memiliki 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.Alih-alih mencari zona UTM yang tepat, Anda dapat membuat proyeksi mercator melintang khusus untuk setiap titik dengan
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
aeqd
sebagai gantitmerc
.sumber
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.
sumber