Memilih fitur "di atas" atau "di bawah" garis menggunakan R

8

Diberi garis dan sekumpulan poin, saya tidak tahu cara menggunakannya sfuntuk mengidentifikasi sisi mana dari setiap titik yang jatuh.

Contoh kecil yang bisa direproduksi berikut, diadaptasi dari pertanyaan yang berbeda

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

output plot

Dalam contoh ini, mudah untuk memverifikasi bahwa titik 2 dan 6 berada di selatan garis, dan sisanya di utara. Bagaimana saya bisa mengotomatiskan pelabelan?

Non sfjawaban berdasarkan menyambut juga.

r sf HAVB
sumber

Jawaban:

8

Jawaban yang diberikan terkait dengan pertanyaan ini Bagaimana cara mengelompokkan objek SpatialPoints untuk mendapatkan titik yang terletak di setiap sisi objek SpatialLines menggunakan R? tetapi menggunakan sfperpustakaan alih-alih sp.

Periksa kode yang dikomentari di bawah ini. 

# 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
sumber
1
Mengubah ke CRS yang berbeda mungkin bukan hal yang tepat untuk dilakukan - mungkin saja si penanya lebih memilih "utara" dan "selatan" merujuk pada garis lintang. Saya tidak tahu mengapa Anda memilih yang itu karena agak diputar dibandingkan dengan lat-long pada saat itu. Saya akan tetap dengan lat-panjang atau mercator. Saya akan menghapus jawaban saya karena ini adalah implementasi bagus dari kerangka saya.
Spacedman
Hai @Spacedman! Apakah maksud Anda mentransformasikannya ke crs 32721? Anda benar, itu tidak perlu. Ini zona UTM 21 Selatan; Saya pikir poin jatuh di zona itu tapi saya tidak yakin. Saya akan mengubah jawaban saya menggunakan crs 4326. Terima kasih!
Guzmán
OP bertanya di atas / bawah dalam konteks 32721, jika bukan itu yang dimaksudkan maka Q harus diperbarui, bukan jawaban untuk Q saat ini (Q mencakup di atas / di bawah dan utara / selatan sehingga benar-benar ambigu atmosfer) .
mdsumner
1
CRS 32721 dalam Q asli saya sampai di sana hanya karena digunakan dalam cuplikan kode yang saya bertujuan ulang ... buruk saya! Saya telah memperbarui pertanyaan saya untuk menghilangkan transformasi yang berpotensi membingungkan dari 4326 menjadi 3721.
HAVB
6

Algoritme garis besar, yang juga memberikan definisi yang lebih kuat tentang "utara atau selatan" garis:

Ubah garis menjadi poligon dengan menambahkan dua segmen garis tambahan dari titik akhir ke Y = -Infinity, atau setidaknya lebih jauh ke selatan daripada titik paling selatan. Kemudian lakukan tes point-in-polygon. Poin dalam poligon berada di sebelah selatan garis.

Ulangi untuk membuat poligon dengan segmen ekstra positif tanpa batas (atau besar). Itu memberi Anda poin di utara garis.

Poin dalam poligon tidak terdefinisi ke utara-selatan sifat garis - mereka timur atau barat dari garis.

Spacedman
sumber