Блог

Главная - Вопросы по программированию - R, sf: пересекаем линии с границами мультиполигонов, извлекаем координаты этих пересечений

R, sf: пересекаем линии с границами мультиполигонов, извлекаем координаты этих пересечений

#r #sf

#r #sf

Вопрос:

Я новичок в SF и stack, надеюсь, мой вопрос достаточно ясен. Я смог создать набор линий, соединяющих 1 точку с набором точек по всей территории США. Я могу преобразовать округа США в мультиполигоны. Моя цель — найти и геолокировать все точки, где созданные мной линии пересекают границы округа.

До сих пор мне удавалось создавать линии из точек:

 points_to_lines <- dt %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  group_by(lineid) %>%
  summarize(do_union = FALSE) %>% lineid
  st_cast("LINESTRING")

Это начало строк

 Simple feature collection with 1628 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 30.1127 ymin: -91.32484 xmax: 37.23671 ymax: -82.31262
epsg (SRID):    4326
proj4string:     proj=longlat  datum=WGS84  no_defs
# A tibble: 1,628 x 2
   lineid                                 geometry
    <int>                         <LINESTRING [°]>
 1      1  (33.51859 -86.81036, 36.16266 -86.7816)
 2      2 (33.51859 -86.81036, 34.61845 -82.47791)

Это начало набора данных округа.

 Reading layer `US_county_1930_conflated' from data source `~/county_gis/1930' using driver `ESRI Shapefile'
Simple feature collection with 3110 features and 18 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -7115608 ymin: -1337505 xmax: 2258244 ymax: 4591848
epsg (SRID):    NA
proj4string:     proj=aea  lat_1=29.5  lat_2=45.5  lat_0=37.5  lon_0=-96  x_0=0  y_0=0  datum=NAD83  units=m  no_defs

Очень наивно я попытался присвоить им обоим одинаковый набор координат, а затем st_intersects их. Не разреженная матрица показывает, что все линии пересекают только один округ.

 gis1930_p <- st_set_crs(gis1930, 4326) %>% st_transform(4326)
st_intersects(points, gis1930_p, sparse=FALSE)

Я также наношу линии поверх округов, но отображается только карта округов США.

 plot(gis1930_p[0], reset = FALSE)
plot(points[0], add = TRUE)

Я был бы очень признателен за любую помощь и, пожалуйста, дайте мне знать, если я смогу предоставить какие-либо дополнительные подробности.

Ответ №1:

Вы не предоставили свои данные, поэтому я собираюсь использовать набор данных, предоставленный в:https://gis.stackexchange.com/a/236922/32531

Главное, что вам нужно, это st_intersection функция:

 library(sf)

line_1 <- st_as_sfc("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
poly_1 <- st_as_sfc("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")
pnts   <- st_intersection(line_1, 
                          st_cast(poly_1, "MULTILINESTRING", group_or_split = FALSE))

plot(poly_1)
plot(line_1, add = TRUE)
plot(pnts, add = TRUE, col = "red", pch = 21)

введите описание изображения здесь

Комментарии:

1. Привет, большое спасибо. Похоже, это действительно то, что я ищу. К сожалению, я не могу заставить это работать. Когда я выполняю st_interesection между округами и одной строкой, я получаю пустой набор данных. Что я замечаю, что отличается то, что в моей настройке у меня есть два объекта, которые привязаны по-разному, может быть? Я должен запустить st_set_crs, чтобы они были одинаковыми (например, 4326). Что странно, так это то, что если я проведу линию по округам, линия не появится.

2. Да, это звучит как проблема с проекцией. Вы должны использовать st_transform для изменения CRS, а не st_set_crs . Я использую mapview пакет для отладки ошибок проекции.

3. большое спасибо, это сработало после того, как я внес два изменения в свой код: 1. Я st_transform, а затем st_set_src(4326). 2. Я переключился на st_as_sf(coords = c("lat","lon")) вместо st_as_sf(coords = c("lon","lat")) ! Также подсказка в mapview была отличной для устранения второй ошибки!!