Suite

Comment trouver la distance d'un point à un polygone, via un autre point, dans R

Comment trouver la distance d'un point à un polygone, via un autre point, dans R


J'essaie de calculer les distances entre un point central et le bord d'un polygone, le long de lignes droites tracées entre le point central et un autre ensemble de points (le cas échéant prolongés au-delà du point extérieur jusqu'au polygone).

Mes données ressemblent à ceci :

library(sp) library(raster) # Center point Center <- data.frame(x= 38.35419, y= 4.483533) Center <- SpatialPointsDataFrame(coords=Centre, data=Centre) projection(Centre) <- CRS('+proj =longlat') crs(Centre) <- "+proj=longlat +datum=WGS84" # Sites Site <- c(seq(1,10,1)) x <- c(38.29013,37.90987,38.29999,38.66006,38.15978 ,38.25978,38.15997,38.24062,38.25005,38.00712) y <- c(5.33006,4.96017,5.23995,4.74010,5.09011,5.09015,5.06016,4.99723,5.06002,4.94634) Sites <- data.frame(Site, x, y) xy <- Sites[2:3] Sites2 <- SpatialPointsDataFrame(coords=xy, data=Sites) projection(Sites2) <- CRS('+proj=longlat') crs(Sites2) <- "+proj=longlat +datum= WGS84" # Polygone Sr1 = Polygone(cbind(c(37.95169, 37.95169, 38.65000, 38.65000), c(5.17153, 3.85601, 3.85601, 5.17153))) Srs1 = Polygones(list(Sr1), "s1") Coque = SpatialPolygons( liste(Srs1), proj4string=CRS("+proj=longlat +datum=WGS84"))

J'ai réussi à calculer la distance du point central à chaque point périphérique en utilisant :

spDistsN1(Sites2, Centre, longlat = TRUE)

Ce qui renvoie la distance en kilomètres.

Pour visualiser les données, j'utilise le code :

library(rgeos) library(maptools) data(wrld_simpl) plot(wrld_simpl, xlim=c(37.5.39.5), ylim=c(4,5), axes=T) plot(Hull, col="green", axes= T, add=T) points(Sites2$x, Sites2$y, pch=20) plot(Centre, add=T, col="yellow", pch="+", cex=1,5)

et ont réussi à tracer des lignes du point central à chaque site :

C_P_lines <- list() for (i in 1:length(Sites$x)) { C_P_lines[i] <- Lines(list( Line(rbind([email protected], Sites[i,2:3]))) , ID=Site.coords[i,1]) } C_P_spLines <- SpatialLines(C_P_lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) plot(C_P_spLines , ajouter=T)

Ce que je dois faire, en fait, c'est a) étendre ces lignes jusqu'au bord du polygone (pour les sites à l'intérieur du polygone), puis b) extraire les points d'intersection des lignes avec le polygone.

a) Jusqu'à présent, je n'ai pas réussi à trouver une solution pour allonger les lignes.

b) Je dois obtenir les résultats indiqués ici : http://uk.mathworks.com/help/map/ref/polyxpoly.html mais dans R. J'ai essayé d'utilisergIntersectiondans legéospaquet,

Edge_points <- gIntersection(C_P_spLines, Hull, byid=T) plot(Edge_points, col="red", add=T)

mais cela ne fait que recadrer les lignes jusqu'à l'étendue du polygone.

Je suis vraiment nouveau dans l'utilisation de R pour le travail SIG.


D'autres utilisateurs pourraient être intéressés de savoir que j'ai réussi à résoudre ce problème.

Afin de prolonger les lignes à l'extérieur du polygone (problème a) :

# Enregistrer les coordonnées du point central séparément A.x <- [email protected][1] A.y <- [email protected][2]

En utilisant l'excellente réponse publiée ici : https://stackoverflow.com/a/7741655/4312169 j'ai appliqué le code :

# Générer des points à l'extérieur du polygone le long de la ligne 'C_P_spLines' (j'ai joué avec le multiplicateur à la fin pour m'assurer que tous les points se retrouvent à l'extérieur du polygone, mais pas trop loin) Out.site.coords <- Sites[, 1:3] pour (i en 1:longueur(Sites$x)) { Out.site.coords$x[i] <- Sites[i,2] + (Sites[i,2] - Ax) / sqrt( (Ax - Sites[i,2])^2 + (Ay - Sites[i,3])^2) * 1 Out.site.coords$y[i] <- Sites[i,3] + (Sites[ i,3] - Ay) / sqrt((Ax - Sites[i,2])^2 + (Ay - Sites[i,3])^2) * 1 } # En faire un SpatialPointsDataFrame xy <- Out.site .coords[2:3] Out.sites <- SpatialPointsDataFrame(coords=xy, data=Out.site.coords) projection(Out.sites) <- CRS('+proj=longlat') crs(Out.sites) < - "+proj=longlat +datum=WGS84" class(Out.sites) Out.sites # Re-tracer avec une étendue plus large (basé sur les coordonnées min et max de "Out.sites") plot(wrld_simpl, xlim=c (37.40.5), ylim=c(2.5,6.6), axes=T) plot(Hull, col="green", add=T) points(Sites$x, Sites$y, pch=20) plot(Centre , add=T, col="yellow", pch="+", cex=1.5) points(Out.sites $x, Out.sites$y, pch=20, col="blue") # Créer un objet SpatialLines pour "centre to out.sites" Radial_Lines <- list() for (i in 1:length(Out.site.coords $x)) { Radial_Lines[i] <- Lines(list( Line(rbind([email protected], Out.site.coords[i,2:3]))) , ID=Out.site.coords[i,1 ]) } # Faites-en un objet SpatialLines Radial_spLines <- SpatialLines(Radial_Lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

j'ai ensuite utiliségIntersectionet quelques sous-ensembles pour extraire les points extérieurs des lignes - voir la réponse ici : /a/154864/55218