如何通过点是否在多边形内来标记点 [英] How to mark points by whether or not they are within a polygon

问题描述

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

I have the longitude and latitude of bird observations across New Zealand, stored in Count.df, under the variables count, longitude, and latitude. However some of these are appearing offshore/in the ocean (it's a Citizen Science data set).

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

I would like to subset these points based off whether they fall outside of the information given in maps::map("nz"), for two reasons:

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

如何根据它们是否落在由 map("nz")划定的组(1:3)的内部/外部，在 Count.df 中创建新变量

How can I create a new variable in Count.df based on whether they fall inside/outside of the groups (1:3) delineated by map("nz"), which is saved as "nzmap.dat"?

谢谢，请尽量简化编码术语.

Thanks, and please try to keep the coding lingo simple.

Count.df 的子集，我有17000多个观测值，这里有10个.请注意，"count"的值与这个问题无关.

Subset of Count.df, I have over 17000 observations, here are 10. Please note that the value of "count" is of no interest for this question.

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"
)

推荐答案

不幸的是，使用 maps 数据在空间上过滤点相当困难，因为该数据实际上是用于绘制而不是做很多事情与空间操作有关.幸运的是，使用 sf 包进行这种空间工作相当容易.

Unfortunately it's rather difficult to use maps data to spatially filter points, because the data is really for plotting and not so much for doing spatial operations with. Luckily, it's fairly easy to use the sf package to do this kind of spatial work.

1. 首先，我从 rnaturalearth 包中获得了一个新西兰边界，这是一个方便的国家边界来源.我进行了一些转换，仅将最大的三个多边形保留在shapefile中，因为否则我们将绘制很多可能与该地图无关的遥远岛屿.
2. 然后，我在新西兰周围产生一些随机点.您可以看到 tibble 只是一列经度和一列纬度.我们使用 st_as_sf 将其转换为点几何，并对其进行绘制，以显示其外观.
3. 最后，我们可以使用 st_within 检查每个点是否在 nz 边界之内. st_within 返回该点在每个点所在的多边形的索引列表，因此我们可以使用 lengths 来获得所需的结果.这里， 0 的任何内容都不在边界内，而 1 的任何内容都在边界内.使用此新的 on_land 属性进行绘制可显示离岸点已适当着色.

1. First I get a New Zealand boundary from the rnaturalearth package, which is a convenient source of country borders. I do a few transformations to keep only the largest three polygons in the shapefile, since otherwise we'd plot a lot of far flung islands that likely aren't relevant to this map.
2. Then I generate some random points around New Zealand. You can see that the tibble is just a column of longitudes and a column of latitudes. We turn this into point geometries with st_as_sf and plot it, to show what it looks like.
3. Last, we can use st_within to check for each point whether or not it is inside the nz boundary. st_within returns a list of indices of polygons that the point is within for each point, so we can use lengths to get the result we want. Here anything that is 0 is not within the boundary and anything that is 1 is within the boundary. Plotting with this new on_land attribute shows that offshore points are appropriately coloured.

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软件包(v0.2.0)创建于2018-05-02.

Created on 2018-05-02 by the reprex package (v0.2.0).

这篇关于如何通过点是否在多边形内来标记点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！