Selecionando recursos "acima" ou "abaixo" de uma linha usando R

8

Dada uma linha e um conjunto de pontos, não consigo descobrir como usar sfpara identificar qual lado da linha cada ponto cai.

Um pequeno exemplo reproduzível segue, adaptado de uma pergunta diferente

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

saída da trama

Neste exemplo, é fácil verificar se os pontos 2 e 6 ficam ao sul da linha e o restante ao norte. Como automatizar a rotulagem?

sfRespostas não baseadas também são bem-vindas.

r sf HAVB
fonte

Respostas:

8

A resposta fornecida está relacionada a esta pergunta Como subconjunto de um objeto SpatialPoints para obter os pontos localizados em cada lado de um objeto SpatialLines usando R? mas usando a sfbiblioteca em vez de sp.

Verifique o código comentado abaixo. 

# 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
fonte
1
Transformar em um CRS diferente pode não ser a coisa certa a fazer - é possível que o questionador prefira que "norte" e "sul" se refiram à latitude. Eu não sei por que você escolheu esse em particular porque é um pouco girado em comparação com o longo nesse momento. Eu ficaria com um comerciante de longa data. Excluirei minha resposta porque esta é uma boa implementação do meu esquema.
Spacedman
Olá @Spacedman! Você quer dizer transformar para crs 32721? Você está certo que não é necessário. É a zona 21 do sul da UTM; Eu acho que os pontos caem nessa zona, mas não tenho certeza. Modificarei minha resposta usando crs 4326. Obrigado!
Guzmán
O OP perguntou acima / abaixo no contexto de 32721, se não for esse o pretendido, o Q deve ser atualizado, não esta resposta para o Q atual. (O Q inclui acima / abaixo e norte / sul, por isso é atm estritamente ambíguo) .
Mdsumner 13/01/19
1
O CRS 32721 no meu Q original chegou lá simplesmente porque foi usado no snippet de código que eu propus novamente ... meu mal! Eu atualizei a minha pergunta para se livrar da transformação potencialmente confusa de 4326 a 3721
HAVB
6

Algoritmo de estrutura de tópicos, que também fornece uma definição mais forte de "norte ou sul" da linha:

Transforme a linha em um polígono adicionando dois segmentos de linha extras dos pontos finais até Y = -Infinity, ou pelo menos mais ao sul que o ponto mais ao sul. Em seguida, faça um teste point-in-polygon. Os pontos no polígono estão ao sul da linha.

Repita para criar um polígono com segmentos extras positivos infinitos (ou grandes). Isso fornece pontos ao norte da linha.

Os pontos em nenhum dos polígonos são indefinidos quanto ao norte-sul da natureza da linha - estão a leste ou oeste da linha.

Spacedman
fonte