如何通过判断点是否在多边形内来标记点

Question

我拥有存储在Count.df中的新西兰各地观鸟记录的经度和纬度数据，分别保存在变量count、longitude和latitude中。然而，其中一些出现在离岸/海洋区域(这是一个公民科学数据集)。

出于以下两个原因，我想根据maps::map("nz")提供的信息将这些点子集化，以确定它们是否在范围之外：

1. 我想在地图上绘制这些离岸点，也许使用不同的颜色
2. 为了进一步分析，我可能会将它们从数据集中删除。

如何基于map("nz")中的(1:3)组内/外的位置创建Count.df中的新变量？"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"
)

- jvan257

---

请编辑问题，将Count.df的至少一个样本作为纯文本包含在内：https://dev59.com/eG025IYBdhLWcg3whGSx?rq=1 - neilfws

Answer 1

很遗憾，使用地图数据来空间过滤点是相当困难的，因为这些数据更多用于绘图而不是执行空间操作。幸运的是，使用sf包来进行这种空间操作相当容易。

1. 首先，我从rnaturalearth包中获取新西兰边界，这是一个方便的国家边界数据源。我对shapefile进行了一些转换，只保留最大的三个多边形，因为否则我们将绘制许多远离主体的岛屿，这可能与此地图无关。
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)

这段内容是由reprex package（v0.2.0）于2018年5月2日创建的。

- Calum You

---

非常感谢您的帮助！我在运行您的代码时遇到了一些问题，因为需要一些额外的包“rnaturalearthhires”和“lwgeom”，但是在这之后一切都正常了。我将您的方法应用于自己的数据，没有遇到任何问题，现在我有了on_land变量，应该能够使用ggplot2和maps可视化我想要的内容。 - jvan257

Answer 2

你可以使用 maps::map.where 进行操作，如果需要的话，可以使用一个只包含三个主要岛屿的自定义地图：

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)

- Alex Deckmyn