J'ai des centaines de points lat-long répartis dans le monde et je dois créer des cercles-polygones autour de chacun d'eux, avec un rayon d'exactement 1000 mètres. Je comprends que les points doivent d'abord être projetés de degrés (lat long) à quelque chose avec des unités de mètre, mais comment cela peut-il être fait sans rechercher et définir manuellement les zones UTM pour chaque point?
Voici un mwe pour le premier point en Finlande.
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"))
Mais avec le deuxième point (Canada) ça ne marche pas (car mauvaise zone UTM).
the.points.projected <- spTransform(the.points.sp[2, ], CRS( "+init=epsg:32635" ))
Comment cela peut-il se faire sans obtenir et spécifier manuellement point par zone UTM? Je n'ai pas plus d'informations par point que lat long.
Mettre à jour:
En utilisant et en combinant les bonnes réponses d'AndreJ et de Mike T, voici le code pour les versions et les graphiques. Ils sont différents à la 4ème décimale environ, mais les deux très bonnes réponses!
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
(Je ne sais pas comment obtenir l'intrigue dans la mise à jour).
la source
Réponses:
Semblable à @AndreJ, mais utilisez une
projection gnomiquedynamique , je veux dire une projection équidistante azimutale dynamique pour encore plus de précision. Une projection AEQ centrée sur chaque point projettera des distances égales dans toutes les directions, comme un cercle tamponné. (Une projection Mercator aura quelques distorsions dans les directions nord et est, car elle s'enroule autour du côté d'un cylindre.)Donc, pour votre premier point autour de la Finlande, la chaîne PROJ.4 ressemblera à ceci:
Vous pouvez donc créer une fonction R pour réaliser cette projection dynamique:
Ensuite, faites quelque chose comme ça pour le Canada (point 2):
la source
projected
est en effet toujours à (0, 0), etbuffered
a les points ± 1000 m dans les directions x et y . S'il était essentiel d'optimiser cela, il suffit de transformer une simple version cartésienne du tampon de l'AEQD dynamique en WGS84.Au lieu de rechercher la bonne zone UTM, vous pouvez créer une projection Mercator transversale personnalisée pour chaque point avec
Tracez le cercle dans cette projection. Les coordonnées du sommet du cercle projeté seront toujours les mêmes, vous devez donc les créer une seule fois. Pour les éléments suivants, attribuez-leur simplement le nouveau CRS personnalisé.
Reprojetez le cercle à EPSG: 4326 pour une utilisation ultérieure.
Dans un rayon de 1000 m, le cercle sera presque exact. Sinon (ou pour les cercles plus grands), utilisez
aeqd
plutôt quetmerc
.la source
Que se passe-t-il si vous adoptez l'approche de la création d'un 1000 mètres en EPSG: 4326 autour de chacun de vos points. Puis convertir l'EPSG: 4326 à votre autre système de coordonnées? L'avantage de projeter le point, c'est que vous n'avez pas à vous soucier de la courbure de la terre avec EPSG: 4326.
la source