How to mark points by whether or not they are within a polygon

Question

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

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. I would like to plot these offshore points on a map, perhaps using a different colour
2. For further analysis I will probably remove them from the data set.

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.

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

Answer 1

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

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

Answer 2

You can do this with maps::map.where and (if you want) a custom map containing only the 3 main islands:

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)