Почему мое пространственное объединение возвращает разные результаты, используя пакет sp и sf? - PullRequest
Новая
       94

Почему мое пространственное объединение возвращает разные результаты, используя пакет sp и sf?

0 голосов
/ 27 июня 2018

В следующем примере я запускаю пространственное соединение для некоторых точечных и полигональных данных, но неожиданно получаю другие результаты при использовании пакета sp от использования пакета sf. Почему это?

Я пытаюсь подсчитать acled баллов в prio квадратах сетки, но, как показано ниже, мои подсчеты различаются между пакетами, даже если я запускаю st_covers объединение из sf, если, насколько мне известно, функционально будет одинаковым используя метод over из sp. 

library(sp) # packageVersion("sp") #> [1] ‘1.2.7’
library(sf) # packageVersion("sf") #> [1] ‘0.6.3’
library(rgdal)
library(maptools)
library(dplyr); library(tibble)

Вот пример данных, с которыми я работаю: 

# prio (polygon squares) and acled (points); in both sp and sf objects:

# prio sf polygons object
priosf <- structure(list(
  CELL_ID = c(180365, 176783, 150830, 145866, 140055), 
  gwno = c(615L, 616L, 432L, 626L, 475L), 
  POP = c(111983.7, 107369.7, 12169.35, 23005.76, 527012.1), 
  prio_country = c("Algeria", "Tunisia", "Mali", "South Sudan", "Nigeria"), 
  geometry = structure(list(structure(list(structure(c(2, 2, 2.5, 2.5, 2, 35, 35.5, 35.5, 35, 35), 
                                                     .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
                            structure(list(structure(c(11, 11, 11.5, 11.5, 11, 32.5, 33, 33, 32.5, 32.5), 
                                                     .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
                            structure(list(structure(c(-5.5, -5.5, -5, -5, -5.5, 14.5, 15, 15, 14.5, 14.5), 
                                                     .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
                            structure(list(structure(c(32.5, 32.5, 33, 33, 32.5, 11, 11.5, 11.5, 11, 11),
                                                     .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
                            structure(list(structure(c(7, 7, 7.5, 7.5, 7, 7, 7.5, 7.5, 7, 7), 
                                                     .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), 
                       class = c("sfc_POLYGON", "sfc"), precision = 0, 
                       bbox = structure(c(-5.5, 7, 33, 35.5), 
                                        .Names = c("xmin", "ymin", "xmax", "ymax"), 
                                        class = "bbox"), 
                       crs = structure(list(epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), 
                                       .Names = c("epsg", "proj4string"), class = "crs"), n_empty = 0L)), 
              .Names = c("CELL_ID", "gwno", "POP", "prio_country", "geometry"),
  row.names = c(NA, -5L), class = c("sf", "tbl_df", "tbl", "data.frame"),
  sf_column = "geometry", agr = structure(c(NA_integer_, NA_integer_, NA_integer_, NA_integer_), 
                                          class = "factor", .Label = c("constant", "aggregate", "identity"), 
                                          .Names = c("CELL_ID", "gwno", "POP", "prio_country")))

# prio sp polygons object
priosp <- as(priosf, 'Spatial')

# acled data
acled <- structure(list(
  EVENT_ID_CNTY = c("ALG3195", "ALG3316", "ALG4228", 
      "ALG4824", "MLI1050", "MLI1144", "MLI1423", "MLI1672", "NIG4606", 
      "NIG4951", "NIG6196", "NIG7661", "NIG9100", "SSD1216", "SSD1504", 
      "SSD3232", "SSD3234", "SSD3231", "SSD3239", "TUN1376", "TUN2597", 
      "TUN3217", "TUN3633"), 
  COUNTRY = c("Algeria", "Algeria", "Algeria", 
              "Algeria", "Mali", "Mali", "Mali", "Mali", "Nigeria", "Nigeria", 
              "Nigeria", "Nigeria", "Nigeria", "South Sudan", "South Sudan", 
              "South Sudan", "South Sudan", "South Sudan", "South Sudan", "Tunisia", 
              "Tunisia", "Tunisia", "Tunisia"), 
  LATITUDE = c(35.2122, 35.4343, 35.2122, 35.2122, 14.8252, 14.8252, 14.7414, 14.8252, 7.3028, 
               7.3028, 7.3028, 7.3028, 7.3588, 11.05, 11.05, 11.05, 11.05, 11.05, 11.05, 32.8487, 32.7149, 32.7149, 32.7149), 
  LONGITUDE = c(2.3189, 2.2166, 2.3189, 2.3189, -5.2547, -5.2547, -5.3282, -5.2547, 7.0382, 7.0382, 7.0382, 7.0382, 7.0994, 32.7, 32.7, 32.7, 32.7, 32.7, 32.7, 11.4309, 11.012, 11.012, 11.012)), 
  row.names = c(NA, -23L), 
  class = c("tbl_df", "tbl", "data.frame"), 
  .Names = c("EVENT_ID_CNTY", "COUNTRY", "LATITUDE", "LONGITUDE"))

# acled sf points object
acledsf <- st_as_sf(
  acled,
  coords = c('LATITUDE', 'LONGITUDE'),
  crs = 4326
)

# acled sp points object
coordinates(acled) <- ~LONGITUDE+LATITUDE
  proj4string(acled) <- proj4string(priosp)
acledsp <- acled; rm(acled)

sp результат пространственного объединения пакетов. Я связал полигоны, которые пересекаются с каждой точкой, соединил результат с точками, а затем посчитал количество CELL_ID (полигонов): 

# sp spatial join:
addPolyDataToPts <- function (points, poly) {
  polysByPoint <- over(points, poly)
  points <- spCbind(points, polysByPoint)
}

acj <- addPolyDataToPts(acledsp, priosp)

(acled_count_sp <- acj@data %>% filter(!is.na(CELL_ID)) %>%
  group_by(CELL_ID, prio_country, POP) %>%
  summarize(acled_sp = n()) %>% arrange(CELL_ID) %>%
  rename(prio_country_sp = prio_country))
#> # A tibble: 5 x 4
#> # Groups:   CELL_ID, prio_country_sp [5]
#>   CELL_ID prio_country_sp     POP acled_sp
#>     <dbl> <chr>             <dbl>    <int>
#> 1 140055. Nigeria         527012.        5
#> 2 145866. South Sudan      23006.        6
#> 3 150830. Mali             12169.        4
#> 4 176783. Tunisia         107370.        4
#> 5 180365. Algeria         111984.        4

Аналогичный sf результат пакета пространственного соединения, где мой столбец подсчета acled_sf отличается от указанного выше столбца acled_sp для всех, кроме одного квадрата многоугольника. (140055; Нигерия): 

# sf spatial join:
(acled_count_sf <- 
  st_join(priosf, acledsf, join = st_covers) %>%
  group_by(CELL_ID, POP, prio_country) %>%
  summarize(acled_sf = n()) %>% ungroup %>% 
  arrange(CELL_ID) %>%
  rename(prio_country_sf = prio_country))
#> although coordinates are longitude/latitude, st_covers assumes that they are planar
#> Simple feature collection with 5 features and 4 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -5.5 ymin: 7 xmax: 33 ymax: 35.5
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 5 x 5
#>   CELL_ID     POP prio_country_sf acled_sf                        geometry
#>     <dbl>   <dbl> <chr>              <int>                   <POLYGON [°]>
#> 1 140055. 527012. Nigeria                5 ((7 7, 7 7.5, 7.5 7.5, 7.5 7, …
#> 2 145866.  23006. South Sudan            4 ((32.5 11, 32.5 11.5, 33 11.5,…
#> 3 150830.  12169. Mali                   1 ((-5.5 14.5, -5.5 15, -5 15, -…
#> 4 176783. 107370. Tunisia                6 ((11 32.5, 11 33, 11.5 33, 11.…
#> 5 180365. 111984. Algeria                1 ((2 35, 2 35.5, 2.5 35.5, 2.5 …

Моя теория работы заключается в том, что один метод связывает значения в неправильном порядке, но я не уверен, какой именно. В моем большем примере я получаю аналогичные значения, но привязываюсь к разным полигонам, т. Е. Точки «2706» сопоставляются с ячейкой 1 для соединения sf и с ячейкой 2 для соединения sp.

(И в некоторых случаях некоторые значения полностью отсутствуют в sf объединении)

Любое понимание того, как или почему мои результаты отличаются таким образом, было бы очень полезно.

1 Ответ

0 голосов
/ 29 июня 2018

Таким образом, мне потребовалось нанести данные в mapview, чтобы выяснить, что здесь происходит, но, по крайней мере, в вашем заданном представлении ваша проблема вызвана тем, что вы указали свою долготу и широту в обратном направлении при создании объекта acledsf. Создано в правильном порядке, и выходы соединения совпадают: 

# acled sf points object
acledsf <- st_as_sf(
  acled,
  coords = c('LONGITUDE', 'LATITUDE'),  ###notice the correct order here
  crs = 4326
) 

# acled sp points object
coordinates(acled) <- c("LONGITUDE", "LATITUDE")
proj4string(acled) <- proj4string(priosp)
acledsp <- acled; rm(acled)


addPolyDataToPts <- function (points, poly) {
  polysByPoint <- over(points, poly)
  points <- spCbind(points, polysByPoint)
}

acj <- addPolyDataToPts(acledsp, priosp)

(acled_count_sp <- acj@data %>% filter(!is.na(CELL_ID)) %>%
    group_by(CELL_ID, prio_country, POP) %>%
    summarize(acled_sp = n()) %>% arrange(CELL_ID) %>%
    rename(prio_country_sp = prio_country))
#> # A tibble: 5 x 4
#> # Groups:   CELL_ID, prio_country_sp [5]
#>   CELL_ID prio_country_sp     POP acled_sp
#>     <dbl> <chr>             <dbl>    <int>
#> 1  140055 Nigeria         527012.        5
#> 2  145866 South Sudan      23006.        6
#> 3  150830 Mali             12169.        4
#> 4  176783 Tunisia         107370.        4
#> 5  180365 Algeria         111984.        4


### sf
(acled_count_sf <- 
    st_join(priosf, acledsf, join = st_covers) %>%
    group_by(CELL_ID, prio_country,  POP) %>%
    summarize(acled_sf = n()) %>% ungroup %>% 
    arrange(CELL_ID) %>%
    rename(prio_country_sf = prio_country))
#> although coordinates are longitude/latitude, st_covers assumes that they are planar
#> Simple feature collection with 5 features and 4 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -5.5 ymin: 7 xmax: 33 ymax: 35.5
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 5 x 5
#>   CELL_ID prio_country_sf     POP acled_sf                        geometry
#>     <dbl> <chr>             <dbl>    <int>                   <POLYGON [°]>
#> 1  140055 Nigeria         527012.        5 ((7 7, 7 7.5, 7.5 7.5, 7.5 7, …
#> 2  145866 South Sudan      23006.        6 ((32.5 11, 32.5 11.5, 33 11.5,…
#> 3  150830 Mali             12169.        4 ((-5.5 14.5, -5.5 15, -5 15, -…
#> 4  176783 Tunisia         107370.        4 ((11 32.5, 11 33, 11.5 33, 11.…
#> 5  180365 Algeria         111984.        4 ((2 35, 2 35.5, 2.5 35.5, 2.5 …
...