Plot spatial polygon on worldmap and extract the coordinates within it (R)

Question

Given a set of coordinates with worldwide distribution, I would like to plot a polygon (preferably a rectangle) on a worldmap and then extract all the coordinates that fall within these plots. The location of the polygon would be placed manually based on point density but their size has to remain constant.

Here, I provide a minimal reproducible example of points randomly distributed across Europe. I also add an orientative image that hopefully helps to understand the desired output. First, I would like to add the polygon on the map and then extract all points within that area. Thanks for any possible help in advance.

#Load libraries
library(rgdal) #v1.5-28
library(rgeos) #v.0.5-9
library(ggplot2) # 3.3.5
library(rworldmap) #plot worldmap v.1.3-6
library(dplyr) #v.1.0.7

#Create dataframe of coordinates that fall in Europe
coord <- data.frame(cbind(runif(1000,-15,45),runif(1000,30,75)))
colnames(coord) <- c("long","lat")

#Exlude ocean points following this post
#https://gis.stackexchange.com/questions/181586/remove-points-which-are-out-of-shapefile-or-raster-extent
URL <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_ocean.zip"
fil <- basename(URL)
if (!file.exists(fil)) download.file(URL, fil)
fils <- unzip(fil)
oceans <- readOGR(grep("shp$", fils, value=TRUE), "ne_110m_ocean",
                  stringsAsFactors=FALSE, verbose=FALSE)


europe_coord <- data.frame(long = coord$long,
                          lat = coord$lat)

coordinates(europe_coord) <- ~long+lat
proj4string(europe_coord) <- CRS(proj4string(oceans))

ocean_points <- over(europe_coord, oceans)

#Add ocean points to dataset
coord$ocean <- ocean_points$featurecla
#Exlude ocean points
europe_land <- coord %>% filter(is.na(ocean))

#Load worldmap
world <- map_data("world")
#Plot europe spatial data
ggplot() + geom_map(data = world, map = world,
                    aes(long, lat, map_id = region), color = "white", 
                    fill = "lightgray", size = 0.1) +
    geom_point(data = europe_land,aes(long, lat),
               alpha = 0.7, size = 0.05) + ylim(0,70) +
    coord_sf(xlim = c(-15, 45), ylim = c(30, 75), expand = FALSE)

Answer 1

You can convert the objects to sf, then use st_intersection to extract the points for a given polygon. First, I create a bounding box (this is where you would enter your coordinates for the polygon extent that you want). Then, I convert the europe_land to an sf object. Then, I use st_intersection to return only the points that fall within the polygon.

library(sf)

bb <-
  c(
    "xmin" = 5.005,
    "xmax" = 12.005,
    "ymin" = 45.008,
    "ymax" = 51.005
  ) %>%
  sf::st_bbox() %>%
  sf::st_as_sfc() %>%
  sf::st_as_sf(crs = 4326) %>%
  sf::st_transform(crs = 4326)

europe_land_sf <- europe_land %>% 
  st_as_sf(coords = c("long", "lat"), dim = "XY") %>% 
  st_set_crs(4326)

bb_extraction <- st_intersection(bb, europe_land_sf)

Output

head(bb_extraction)

Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.136638 ymin: 45.69418 xmax: 11.78427 ymax: 49.51203
Geodetic CRS:  WGS 84
  ocean                         x
1  <NA> POINT (7.136638 46.90647)
2  <NA> POINT (9.035649 45.81661)
3  <NA> POINT (10.69865 45.91507)
4  <NA> POINT (11.78427 48.55559)
5  <NA> POINT (10.90349 45.69418)
6  <NA> POINT (9.417477 49.51203)

For plotting, once you have created the sf object, then you can use geom_sf in your existing code to add the polygon to the world map.

ggplot() +
  geom_map(
    data = world,
    map = world,
    aes(long, lat, map_id = region),
    color = "white",
    fill = "lightgray",
    size = 0.1
  ) +
  geom_point(data = europe_land,
             aes(long, lat),
             alpha = 0.7,
             size = 0.05) + ylim(0, 70) +
  geom_sf(data = bb, colour = "black", fill = NA) +
  coord_sf(xlim = c(-15, 45),
           ylim = c(30, 75),
           expand = FALSE)

Output