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

Question

У меня есть долгота и широта наблюдений за птицами по всей Новой Зеландии, сохраненные в Count.df, под переменными count, longitude и latitude. Однако некоторые из них появляются в море / в океане (это набор данных Citizen Science).

Я хотел бы установить эти точки на основании того, выходят ли они за пределы информации, указанной в maps::map("nz"), по двум причинам:

1. Я хотел бы нанести эти оффшорные точки на карту, возможно, используя другой цвет
2. Для дальнейшего анализа я, вероятно, удалю их из набора данных.

Как я могу создать новую переменную в Count.df в зависимости от того, попадают ли они в / вне групп (1: 3), обозначенных map("nz"), которые сохраняются как "nzmap.dat"?

Спасибо, и, пожалуйста, постарайтесь, чтобы код языка был простым.

Подмножество Count.df, у меня более 17000 наблюдений, здесь 10. Обратите внимание, что значение "count" не представляет интереса для этого вопроса.

df <- readr::read_table2(
  "count longitude latitude
3 174.7889 -41.24632
1 174.7764 -41.25923
4 176.8865 -39.67894
1 174.7656 -36.38182
2 175.5458 -37.13479
2 175.5458 -37.13479
1 176.8862 -39.67853
1 170.6101 -45.84626
5 174.9162 -41.25709
2 176.8506 -39.51831"
)

Answer 1

Вы можете сделать это с помощью maps::map.where и (если хотите) пользовательской карты, содержащей только 3 основных острова:

nz3 = map("nz",c("North","South","Stewart"), fill=TRUE, plot=FALSE)
points$onland = is.na(map.where(nz3, points$x, points$y))
plot(points, col = ifelse(points$onland, 1, 2), pch=19)

Answer 2

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

1. Сначала я получаю новозеландскую границу из пакета rnaturalearth, который является удобным источником границ страны. Я делаю несколько преобразований, чтобы сохранить только три самых больших полигона в шейп-файле, так как в противном случае мы бы построили множество удаленных островов, которые, вероятно, не имеют отношения к этой карте.
2. Затем я генерирую несколько случайных точек вокруг Новой Зеландии. Вы можете видеть, что tibble - это просто столбец долгот и столбец широт. Мы превращаем это в точечные геометрии с помощью st_as_sf и наносим на график, чтобы показать, как это выглядит.
3. Наконец, мы можем использовать st_within для проверки каждой точки, находится ли она внутри границы nz. st_within возвращает список индексов полигонов, в которых находится точка для каждой точки, поэтому мы можем использовать lengths, чтобы получить желаемый результат. Здесь все, что является 0, не находится в пределах границы, и все, что является 1, находится в пределах границы. Графики с этим новым атрибутом on_land показывают, что морские точки соответствующим образом окрашены.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(rnaturalearth)

nz <- ne_countries(country = "New Zealand", returnclass = "sf", scale = "large") %>%
  st_cast("POLYGON") %>%
  mutate(area = st_area(geometry)) %>%
  top_n(3, area)
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all sub-
#> geometries for which they may not be constant

points <- tibble(
  x = runif(1000, st_bbox(nz)[1], st_bbox(nz)[3]),
  y = runif(1000, st_bbox(nz)[2], st_bbox(nz)[4])
)
points
#> # A tibble: 1,000 x 2
#>        x     y
#>    <dbl> <dbl>
#>  1  167. -44.5
#>  2  175. -40.9
#>  3  177. -43.8
#>  4  167. -44.8
#>  5  173. -39.3
#>  6  173. -42.1
#>  7  176. -41.9
#>  8  171. -44.9
#>  9  173. -41.2
#> 10  174. -39.5
#> # ... with 990 more rows
points <- st_as_sf(points, coords = c("x", "y"), crs = 4326)
plot(nz$geometry, col = "red")
plot(points, pch = 19, cex = 1, add = TRUE)

points <- points %>% mutate(on_land = lengths(st_within(points, nz)))
#> although coordinates are longitude/latitude, st_within assumes that they are planar
plot(nz$geometry, col = "red")
plot(points, pch = 19, cex = 1, add = TRUE)

Создано в 2018-05-02 пакетом Представить (v0.2.0).