Auswählen von Merkmalen "über" oder "unter" einer Linie mit R.

8

Angesichts einer Linie und einer Reihe von Punkten kann ich nicht herausfinden, wie ich sfidentifizieren kann, auf welche Seite der Linie jeder Punkt fällt.

Es folgt ein kleines reproduzierbares Beispiel, das aus einer anderen Frage stammt

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

Plotausgabe

In diesem Beispiel kann leicht überprüft werden, ob die Punkte 2 und 6 südlich der Linie und der Rest nördlich liegen. Wie kann ich die Kennzeichnung automatisieren?

Nicht sfbasierte Antworten sind ebenfalls willkommen.

r sf HAVB
quelle

Antworten:

8

Die bereitgestellte Antwort bezieht sich auf diese Frage. Wie kann ein SpatialPoints-Objekt untergeordnet werden, um die Punkte auf jeder Seite eines SpatialLines-Objekts mithilfe von R zu ermitteln? aber mit sfBibliothek anstelle von sp.

Überprüfen Sie den kommentierten Code unten. 

# 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
quelle
1
Die Umstellung auf ein anderes CRS ist möglicherweise nicht das Richtige. Möglicherweise möchte der Fragesteller lieber "Nord" und "Süd" als Breitengrad. Ich weiß nicht, warum Sie sich für dieses entschieden haben, weil es zu diesem Zeitpunkt im Vergleich zu lat-long etwas gedreht ist. Ich würde bei Lat-Long oder einem Mercator bleiben. Ich werde meine Antwort löschen, da dies eine schöne Implementierung meiner Gliederung ist.
Spacedman
Hallo @Spacedman! Meinen Sie die Umwandlung in crs 32721? Sie haben Recht, das ist nicht notwendig. Es ist UTM Zone 21 South; Ich denke, dass die Punkte in diese Zone fallen, aber ich bin nicht sicher. Ich werde meine Antwort mit crs 4326 ändern. Danke!
Guzmán
Das OP fragte oben / unten im Zusammenhang mit 32721: Wenn dies nicht beabsichtigt war, sollte das Q aktualisiert werden, nicht diese Antwort auf das aktuelle Q. (Das Q enthält oben / unten und Nord / Süd, daher ist es streng mehrdeutig atm) .
Mdsumner
1
CRS 32721 in meinem ursprünglichen Q ist einfach dort angekommen, weil es in dem Code-Snippet verwendet wurde, das ich neu geplant habe ... mein schlechtes! Ich habe meine Frage aktualisiert, um die möglicherweise verwirrende Umwandlung von 4326 auf 3721 zu beseitigen.
HAVB
6

Gliederungsalgorithmus, der auch eine stärkere Definition von "Norden oder Süden" der Linie gibt:

Verwandeln Sie die Linie in ein Polygon, indem Sie zwei zusätzliche Liniensegmente von den Endpunkten bis zu Y = -Infinity oder zumindest weiter südlich als der südlichste Punkt hinzufügen. Führen Sie dann einen Punkt-in-Polygon-Test durch. Punkte im Polygon befinden sich südlich der Linie.

Wiederholen Sie diesen Vorgang, um ein Polygon mit unendlich (oder großen) positiven zusätzlichen Segmenten zu erstellen. Das gibt Ihnen Punkte nördlich der Linie.

Punkte in keinem Polygon sind hinsichtlich ihrer Nord-Süd-Linie undefiniert - sie befinden sich östlich oder westlich der Linie.

Spacedman
quelle