r - 如何通过点是否在多边形内来标记点

Question

我有新西兰鸟类观测的经度和纬度，存储在 Count.df 中。 , 在变量 count 下, longitude , 和 latitude .然而，其中一些出现在近海/海洋中(这是一个公民科学数据集)。

我想根据这些点是否超出 maps::map("nz") 中提供的信息来对这些点进行子集化。 ，有两个原因:

我想在 map 上绘制这些离岸点，也许使用不同的颜色

为了进一步分析，我可能会将它们从数据集中删除。

如何在 Count.df 中创建新变量基于它们是否属于 map("nz") 所描述的组 (1:3) 之内/之外，保存为“nzmap.dat”？

谢谢，请尽量保持编码术语简单。
Count.df 的子集, 我有超过 17000 次观察，这里有 10 次。请注意“计数”的值对这个问题没有意义。

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数据到空间过滤点，因为数据实际上是用于绘图而不是用于进行空间操作。幸运的是，使用 sf 相当容易。包来做这种空间工作。

首先，我从 rnaturalearth 获得新西兰边界。包，这是一个方便的国界来源。我做了一些转换以仅保留 shapefile 中最大的三个多边形，否则我们将绘制许多可能与此 map 无关的遥远岛屿。

然后我在新西兰周围生成一些随机点。您可以看到 tibble只是一列经度和一列纬度。我们使用 st_as_sf 将其转换为点几何。并绘制它，以显示它的样子。

最后，我们可以使用 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 由 reprex package (v0.2.0)。