Sélection d'entités «au-dessus» ou «en dessous» d'une ligne à l'aide de R

8

Étant donné une ligne et un ensemble de points, je ne peux pas comprendre comment utiliser sfpour identifier de quel côté de la ligne chaque point tombe.

Un petit exemple reproductible suit, adapté d'une question différente

# Load Libraries ----------------------------------------------------------

library('sf')

# Test data ---------------------------------------------------------------

points.df <- data.frame(
    'x' = c(-53.50000, -54.15489, -54.48560, -52.00000, -52.57810, -49.22097, -48.00000),
    'y' = c(-38.54859, -41.00000, -38.80000, -38.49485, -38.00000, -40.50000, -37.74859)
)


line.df <- data.frame(
    'x' = c(-54.53557, -52.00000, -50.00000, -48.00000, -46.40190),
    'y' = c(-39.00000, -38.60742, -38.08149, -38.82503, -37.00000)
)

# Create 'sf' objects -----------------------------------------------------

points.sf <- st_as_sf(points.df, coords = c("x", "y"))

st_crs(points.sf) <- st_crs(4326) # assign crs

line.sf <- st_sf(id = 'L1', st_sfc(st_linestring(as.matrix(line.df), dim = "XY")))
st_crs(line.sf) <- st_crs(4326) # assign crs


# Plots -------------------------------------------------------------------

xmin <- min(st_bbox(points.sf)[1], st_bbox(line.sf)[1])
ymin <- min(st_bbox(points.sf)[2], st_bbox(line.sf)[2])
xmax <- max(st_bbox(points.sf)[3], st_bbox(line.sf)[3])
ymax <- max(st_bbox(points.sf)[4], st_bbox(line.sf)[4])

plot(points.sf, pch = 19, xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin,xmax), ylim = c(ymin,ymax), graticule = st_crs(4326), axes = TRUE)

plot(line.sf, col = "#C72259", add = TRUE)
text(st_coordinates(points.sf), as.character(1:7), pos = 3)

sortie de tracé

Dans cet exemple, il est facile de vérifier que les points 2 et 6 tombent au sud de la ligne et le reste au nord. Comment automatiser l'étiquetage?

Les sfréponses non fondées sont également les bienvenues.

r sf HAVB
la source

Réponses:

8

La réponse fournie est liée à cette question Comment sous-ensemble un objet SpatialPoints pour obtenir les points situés de chaque côté d'un objet SpatialLines en utilisant R? mais en utilisant la sfbibliothèque au lieu de sp.

Vérifiez le code commenté ci-dessous. 

# Load Libraries ----------------------------------------------------------

library('sf')

# Test data ---------------------------------------------------------------

points.df <- data.frame(
  'x' = c(-53.50000, -54.15489, -54.48560, -52.00000, -52.57810, -49.22097, -48.00000),
  'y' = c(-38.54859, -41.00000, -38.80000, -38.49485, -38.00000, -40.50000, -37.74859),
  'id' = as.character(c(1:7))
)


line.df <- data.frame(
  'x' = c(-54.53557, -52.00000, -50.00000, -48.00000, -46.40190),
  'y' = c(-39.00000, -38.60742, -38.08149, -38.82503, -37.00000)
)

# Create 'sf' objects -----------------------------------------------------

points.sf <- st_as_sf(points.df, coords = c("x", "y"))

st_crs(points.sf) <- st_crs(4326) # assign crs

line.sf <- st_sf(id = 'L1', st_sfc(st_linestring(as.matrix(line.df), dim = "XY")))
st_crs(line.sf) <- st_crs(4326) # assign crs    

# Plots -------------------------------------------------------------------

xmin <- min(st_bbox(points.sf)[1], st_bbox(line.sf)[1])
ymin <- min(st_bbox(points.sf)[2], st_bbox(line.sf)[2])
xmax <- max(st_bbox(points.sf)[3], st_bbox(line.sf)[3])
ymax <- max(st_bbox(points.sf)[4], st_bbox(line.sf)[4])

plot(points.sf, pch = 19, xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin,xmax), ylim = c(ymin,ymax), graticule = st_crs(4326), axes = TRUE)

plot(line.sf, col = "#272822", lwd = 2, add = TRUE)
text(st_coordinates(points.sf), as.character(points.sf$id), pos = 3)

map1

# Create Polygons from line -----------------------------------------------

# Add x and y offsets (in degrees units)
offsetX <- 0
offsetY <- 3

polySideUp <- rbind(c(st_bbox(line.sf)['xmax'] + offsetX, 
                       st_bbox(line.sf)['ymax'] + offsetY),
                     c(st_bbox(line.sf)['xmin'] - offsetX, 
                       st_bbox(line.sf)['ymax'] + offsetY),
                     as.data.frame(st_coordinates(line.sf))[,c(1,2)],
                     c(st_bbox(line.sf)['xmax'] + offsetX, 
                       st_bbox(line.sf)['ymax'] + offsetY))

polySideDown <- rbind(c(st_bbox(line.sf)['xmax'] + offsetX, 
                       st_bbox(line.sf)['ymin'] - offsetY),
                     c(st_bbox(line.sf)['xmin'] - offsetX, 
                       st_bbox(line.sf)['ymin'] - offsetY),
                     as.data.frame(st_coordinates(line.sf))[,c(1,2)],
                     c(st_bbox(line.sf)['xmax'] + offsetX, 
                       st_bbox(line.sf)['ymin'] - offsetY))

# Create sf objects
polySideUp <- st_sf("id" = 'sideUp', st_sfc(st_polygon(list(as.matrix(polySideUp))), crs = 4326))
polySideDown <- st_sf("id" = 'sideDown', st_sfc(st_polygon(list(as.matrix(polySideDown))), crs = 4326))

# Plot
plot(polySideUp, xlab = "Longitude", ylab = "Latitude", col = "#C72259", 
     xlim = c(xmin - offsetX, xmax + offsetX), ylim = c(ymin - offsetY, ymax + offsetY), graticule = st_crs(4326), axes = TRUE)
plot(polySideDown, col = "#53A8BD", add = TRUE)
plot(points.sf$geometry, pch = 19, add = TRUE)
plot(line.sf, col = "#272822", lwd = 2, add = TRUE)
text(st_coordinates(points.sf), as.character(points.sf$id), pos = 3)

map2

# Select points in side up
pointsInSideUp <- st_intersection(points.sf, polySideUp)

print(pointsInSideUp)

print1

# Select points in side down
pointsInSideDown <- st_intersection(points.sf, polySideDown)

print(pointsInSideDown)

print2

# Plot intersection
plot(polySideUp, xlab = "Longitude", ylab = "Latitude", col = "#C72259", 
     xlim = c(xmin - offsetX, xmax + offsetX), ylim = c(ymin - offsetY, ymax + offsetY), graticule = st_crs(4326), axes = TRUE)
plot(polySideDown, col = "#53A8BD", add = TRUE)
plot(pointsInSideUp, pch = 19, col = "#53A8BD", add = TRUE)
plot(pointsInSideDown, pch = 19, col = "#C72259", add = TRUE)
plot(line.sf, lwd = 2, col = "#272822", add = TRUE)
text(st_coordinates(points.sf), as.character(points.sf$id), pos = 3)

map3

Guzmán
la source
1
La transformation en un CRS différent n'est peut-être pas la bonne chose à faire - il est possible que le questionneur préfère que "nord" et "sud" se réfèrent à la latitude. Je ne sais pas pourquoi vous avez choisi celui-ci en particulier, car il tourne un peu par rapport au lat-long à ce moment-là. Je resterais avec lat-long ou un mercator. Je vais supprimer ma réponse car c'est une belle mise en œuvre de mon plan.
Spacedman
Salut @Spacedman! Voulez-vous dire la transformation en crs 32721? Vous avez raison, ce n'est pas nécessaire. C'est la zone UTM 21 Sud; Je pense que les points tombent dans cette zone mais je ne suis pas sûr. Je modifierai ma réponse en utilisant crs 4326. Merci!
Guzmán
L'OP a demandé ci-dessus / ci-dessous dans le contexte de 32721, si ce n'est pas ce qui était prévu, le Q devrait être mis à jour, pas cette réponse au Q. actuel (le Q comprend ci-dessus / ci-dessous et nord / sud, donc c'est strictement ambigu) .
mdsumner
1
CRS 32721 dans mon Q d'origine est arrivé simplement parce qu'il a été utilisé dans l'extrait de code que j'ai redéfini ... mon mauvais! J'ai mis à jour ma question pour me débarrasser de la transformation potentiellement déroutante de 4326 à 3721.
HAVB
6

Algorithme de contour, qui donne également une définition plus forte du "nord ou sud" de la ligne:

Transformez la ligne en polygone en ajoutant deux segments de ligne supplémentaires depuis les points d'extrémité jusqu'à Y = -Infinity, ou au moins plus au sud que le point le plus au sud. Effectuez ensuite un test de point dans le polygone. Les points du polygone sont au sud de la ligne.

Répétez l'opération pour créer un polygone avec des segments supplémentaires positifs à l'infini (ou grands). Cela vous donne des points au nord de la ligne.

Les points dans aucun des polygones ne sont pas définis quant à leur nature nord-sud de la ligne - ils sont à l'est ou à l'ouest de la ligne.

Spacedman
la source