Как я могу выполнить пространственное соединение с пакетом sf с помощью st_join ()

Вот игрушечный пример, с которым я боролся

# Make points
point1 <- c(.5, .5)
point2 <- c(.6, .6)
point3 <- c(3, 3)
mpt <- st_multipoint(rbind(point1, point2, point3))  # create multipoint

# Make polygons
square1 <- rbind(c(0, 0), c(1, 0), c(1,1), c(0, 1), c(0, 0))
square2 <- rbind(c(0, 0), c(2, 0), c(2,2), c(0, 2), c(0, 0))
square3 <- rbind(c(0, 0), c(-1, 0), c(-1,-1), c(0, -1), c(0, 0))
mpol <- st_multipolygon(list(list(square1), list(square2), list(square2)))  # create multipolygon

# Convert to class' sf'
pts <- st_sf(st_sfc(mpt))
polys <- st_sf(st_sfc(mpol))

# Determine which points fall inside which polygons
st_join(pts, polys, join = st_contains)

Последняя строка производит

Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : 
  cannot coerce class "c("sfc_MULTIPOINT", "sfc")" to a data.frame

Как я могу выполнить пространственное соединение, чтобы определить, какие точки и какие полигоны попадают?


person Ben    schedule 05.05.2017    source источник
comment
Не могли бы вы пояснить, что вы подразумеваете под пространственным объединением? Каков был бы ожидаемый результат?   -  person lbusett    schedule 09.05.2017
comment
Учитывая набор многоугольников и набор точек, создайте сопоставление (PointId, PolygonId), которое указывает, какие точки и какие многоугольники содержатся.   -  person Ben    schedule 09.05.2017
comment
Недавно я написал это руководство для sf package, чтобы помочь себе и другим понять основные концепции. Понимание основ - ключ к решению конкретных проблем, подобных той, что была у меня здесь.   -  person Ben    schedule 02.08.2017


Ответы (3)


Я также работаю над особенностями пакета sf, поэтому извиняюсь, если это неверно или есть способы получше. Я думаю, что одна проблема здесь в том, что при построении геометрии, как в вашем примере, вы не получаете то, что думаете:

> pts
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
                     st_sfc.mpt.
1 MULTIPOINT(0.5 0.5, 0.6 0.6...

> polys
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
                    st_sfc.mpol.
1 MULTIPOLYGON(((0 0, 1 0, 1 ...

Вы можете видеть, что у вас есть только одна «функция» как в pts, так и в polys. Это означает, что вы строите один объект-мультиполигон (то есть, полигон, состоящий из 3 частей), а не три разных полигона. То же самое и с очками.

Немного покопавшись, я нашел другой (и, на мой взгляд, более простой) способ построения геометрии с использованием обозначения WKT:

polys <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))",
                     "POLYGON((0 0 , 0 2 , 2 2 , 2 0, 0 0 ))", 
                     "POLYGON((0 0 , 0 -1 , -1 -1 , -1 0, 0 0))")) %>% 
  st_sf(ID = paste0("poly", 1:3))    

pts <- st_as_sfc(c("POINT(0.5 0.5)",
                   "POINT(0.6 0.6)",
                   "POINT(3 3)")) %>%
  st_sf(ID = paste0("point", 1:3))

> polys
Simple feature collection with 3 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -1 ymin: -1 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
     ID                              .
1 poly1 POLYGON((0 0, 0 1, 1 1, 1 0...
2 poly2 POLYGON((0 0, 0 2, 2 2, 2 0...
3 poly3 POLYGON((0 0, 0 -1, -1 -1, ...

> pts
Simple feature collection with 3 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
      ID              .
1 point1 POINT(0.5 0.5)
2 point2 POINT(0.6 0.6)
3 point3     POINT(3 3)

вы можете видеть, что теперь и polys, и pts имеют три функции.

Теперь мы можем найти «матрицу пересечений», используя:

# Determine which points fall inside which polygons
pi <- st_contains(polys,pts, sparse = F) %>% 
  as.data.frame() %>% 
  mutate(polys = polys$ID) %>% 
  select(dim(pi)[2],1:dim(pi)[1])
colnames(pi)[2:dim(pi)[2]] = levels(pts$ID)

> pi
  polys point1 point2 point3
1 poly1   TRUE   TRUE  FALSE
2 poly2   TRUE   TRUE  FALSE
3 poly3  FALSE  FALSE  FALSE

это означает (как указано в @symbolixau в комментариях), что многоугольники 1 и 2 содержат точки 1 и 2, а многоугольник 3 не содержит никаких точек. Вместо этого точка 3 не содержится ни в одном многоугольнике.

HTH.

person lbusett    schedule 05.05.2017
comment
Очень интересно. Надеюсь, что и другие будут участвовать в этом. Я бы хотел, чтобы в sf документации были лучшие примеры st_join и создания простых функций с нуля. - person Ben; 05.05.2017
comment
Я тоже считаю, что построение функций в виде списков немного обременительно. Однако это правда, что мне редко нужно строить геометрию с нуля, и увеличение скорости по отношению к sp действительно велико! - person lbusett; 05.05.2017
comment
Я думаю, что результат говорит о том, что точки 1 и 2 находятся в обоих полигонах 1 и 2 (они перекрываются), а точка 3 не входит ни в один многоугольник, а многоугольник 3 не содержит никаких точек. (Я попробовал тот же вопрос, используя другой подход, отличный от sf, чтобы помочь интерпретировать результаты) - person SymbolixAU; 09.05.2017
comment
Если вам нужна более подробная документация о пространственных функциях, вы можете посмотреть здесь postgis.net/docs/reference .html - person troh; 13.05.2017

Я вижу другой результат:

> # Determine which points fall inside which polygons
> st_join(pts, polys, join = st_contains)
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
                        geometry
1 MULTIPOINT(0.5 0.5, 0.6 0.6...

это было с последней версией CRAN sf?

person Edzer Pebesma    schedule 02.06.2017
comment
Спасибо за ответ и отличный пакет. На момент написания этого вопроса это была самая последняя версия на CRAN. Теперь у меня 0.4-3, и я все еще получаю ошибку в том же месте (хотя и другое сообщение). Похоже, что виноват этот вызов pts <- st_sf(st_sfc(mpt)), и я предполагаю, что это потому, что простые функции не могут иметь просто геометрию - они также должны иметь данные. Вот почему что-то вроде pts <- st_sf(Dummy=1, st_sfc(mpt)) предотвращает ошибку. Недавно я подробно рассмотрел виньетки, которые во многом прояснили мою путаницу. - person Ben; 02.06.2017
comment
Спасибо - это было исправлено в версии для разработки некоторое время назад. - person Edzer Pebesma; 03.06.2017

Обратите внимание, что исходный набор мультиточечных и мультиполигональных объектов может быть преобразован в точку и многоугольник без создания новых объектов:

st_contains(polys %>% st_cast("POLYGON"), pts %>% st_cast("POINT"), sparse = F)
#      [,1]  [,2]  [,3]
#[1,]  TRUE  TRUE FALSE
#[2,]  TRUE  TRUE FALSE
#[3,] FALSE FALSE FALSE
person sebdalgarno    schedule 01.03.2018