Geometri c операций над sf с двумя геометриями в R - PullRequest
Новая
       34

Geometri c операций над sf с двумя геометриями в R

0 голосов
/ 07 марта 2020

У меня есть набор данных sf с двумя столбцами геометрии. Вот как это выглядит: 

> trip_geo

   dstid sourceid                    dest_geom                    source_geom
1    1        1   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.607985 5...
2    1        2   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.57022 51...
3    1        3   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.593213 5...
4    1        4   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.608686 5...
5    1        5   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.512852 5...

Активная геометрия - dest_geom.

Каждая строка соответствует поездке между окрестностями. Для каждой поездки я хочу знать, какие районы l ie на пути поездки. То есть, если бы вы нарисовали прямую линию между source_geom и dest_geom для каждой строки, какая геометрия коснулась бы этой прямой линии? Я хочу получить все соприкасающиеся геометрии для строки, а затем объединить их.

У меня есть другой набор данных с геометрией, соответствующей каждому идентификатору: 

> id_geo

   dstid                       geometry
1      1 MULTIPOLYGON (((-2.607985 5...
2      2 MULTIPOLYGON (((-2.57022 51...
3      3 MULTIPOLYGON (((-2.593213 5...
4      4 MULTIPOLYGON (((-2.608686 5...
5      5 MULTIPOLYGON (((-2.512852 5...

Я представляю, что первый шаг было бы определить линию между центроидами source_geom и dest_geom для каждой поездки. Затем создайте sf, где столбец геометрии содержит список полигонов, которые касаются линии (я не знаю, возможно ли иметь несколько геометрий в одном столбце в sf). Затем объедините геометрии, содержащиеся в одной строке / списке.

Я не думаю, что подход к решению проблемы является правильным, поскольку, насколько мне известно, нельзя выполнять операции над двумя геометриями SF, например, определение линии. Кроме того, я не знаю, как бы интегрировать список в фрейм данных / sf.

Если бы вы могли предложить более реалистичный способ решения проблемы, я был бы очень признателен.

Ответы [ 2 ]

3 голосов
/ 08 марта 2020

Нет ничего плохого в наличии двух геометрий во фрейме данных sf, но только одна из них может быть геометрией, если она используется с функциями, которые принимают неявную геометрию, например, st_centroid(foo), которая получает центроиды заданной геометрии.

Для других геометрий вы можете работать с этим именованным столбцом, например st_centroid(foo$source_geom).

Для фрейма данных nc с двумя геометриями многоугольника вы можете сначала вычислить линию между центроидами объединение точек, чтобы сделать MULTIPOINT, а затем приведение их к LINESTRING. Например, для первой строки: 

> st_cast(st_union(c(st_centroid(nc$source_geom[1]), st_centroid(nc$dest_geom[1]))),"LINESTRING")
Geometry set for 1 feature 
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -81.49826 ymin: 36.42228 xmax: -77.41056 ymax: 36.4314
epsg (SRID):    NA
proj4string:    NA
LINESTRING (-81.49826 36.4314, -77.41056 36.42228)

Вы должны сделать это построчно, иначе вы в конечном итоге будете работать со всем вектором геометрии.

Полный пример. Сделайте library(spdep) и example(poly2nb), чтобы получить nc.sids.

Сначала разбейте его на два столбца и 5 случайных строк: 

> nc = nc.sids[,c("NAME","FIPS")]
> nc = nc[sample(1:nrow(nc.sids),5),]
> nc
Simple feature collection with 5 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -83.73952 ymin: 34.36442 xmax: -78.16968 ymax: 36.54607
epsg (SRID):    NA
proj4string:    NA
         NAME  FIPS                       geometry
82 Cumberland 37051 MULTIPOLYGON (((-78.49929 3...
96     Bladen 37017 MULTIPOLYGON (((-78.2615 34...
13  Granville 37077 MULTIPOLYGON (((-78.74912 3...
78      Macon 37113 MULTIPOLYGON (((-83.10629 3...
14     Person 37145 MULTIPOLYGON (((-78.8068 36...

Допустим, функция 1 переходит к функции 2, 2 к 3 и т. Д. c. Создайте новый столбец геометрии: 

> nc$dest_geom = nc$geometry[c(2,3,4,5,1)]

Теперь сделайте линию, соединяющую центроиды: 

> nc$join_geom = st_sfc(sapply(1:nrow(nc),function(i)st_cast(st_union(c(st_centroid(nc$geometry[i]), st_centroid(nc$dest_geom[i]))),"LINESTRING")))

Участок:

centroids joined

> plot(nc$geom)
> plot(nc$join_geom,add=TRUE,lty=2)
0 голосов
/ 16 марта 2020

Ответ @Spacedman был идеальным для рисования линий между интересующими полигонами. Ниже приведен код для (1) рисования линий и (2) выбора полигонов, которые пересекают каждую линию и объединяют выборку. 

#Draw a line between each source and destination
library(sf)
trip_geo$join_geom = st_sfc(sapply(1:nrow(trip_geo),function(j)
+ st_cast(st_union(c(st_centroid(trip_geo$dest_geom[j]),
+ st_centroid(trip_geo$source_geom[j]))),"LINESTRING")), crs=4326)

#For each trip, select the polygons that intersect the line, then union them together
df=data.frame()
for (i in 1:nrow(trip_geo)){
    id_geo$touch=apply(st_intersects(id_geo, trip_geo$join_geom[i]), 1, any)
    touch=subset(id_geo, touch==T)
    touch=st_union(touch)
    df[i,1]=touch #append the unioned polygons to an empty dataframe
  }

#Add the unioned polygons to the original dataset
sf=st_sf(df, crs=4326) 
bound=cbind(trip_geo, sf)

Вычисления занимают некоторое время (около 1 минуты для рисования линий и 1 минуты для l oop, с 3500 полигонами), и я уверен, что мой код можно улучшить.

...