Counting how many point are in the polygon of region [R]

Question

I would like to ask how to calculate number of point that are in some region when we have longtitue and latitude variables of point and polygon of country and its regions.

I provided example below: I would like to calculate how many point are in what regions (including zero when there is no point) and than add this variables to data2@data object so one can use count measures to fill each regions polygons.

library(leaflet)
library(tidyverse)
set.seed(101)

URL2 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_FRA_2_sp.rds"
data2 <- readRDS(url(URL2))


URL3 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_ESP_2_sp.rds"
data3 <- readRDS(url(URL3))

URL4 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_PRT_2_sp.rds"
data4 <- readRDS(url(URL4))

URL5 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_GBR_2_sp.rds"
data5 <- readRDS(url(URL5))

random_long_lat <- 
  data.frame(
    long = sample(runif(300, min = -2.5, max = 15.99), replace = F),
    lat = sample(runif(300, min = 42.69, max = 59.75), replace = F)
  )

all <- rbind(data2, data3, data4, data5)

leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data=all, stroke = TRUE, color = "black", weight="", smoothFactor = 0.95, 
              fill = F) %>% 
  addCircles(lng = random_long_lat$long, lat = random_long_lat$lat)


# Here add new variable called count, that would count how many point are in the region
all@data  

I would like the result so one can calculate something like this:

all@data <- 
  all@data %>% 
  mutate(counts = rnorm(nrow(all), 100, sd = 20))

cuts <- c(0, 5, 20, 40, 80, 150, max(all@data$counts))
cuts <- colorBin("Greens", domain = all$counts, bins = cuts)

leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data=all, stroke = TRUE, color = "white", weight="", smoothFactor = 0.95, 
              fillOpacity = 0.65, fillColor = ~cuts(all$counts)) %>% 
  addLegend(pal = cuts, 
            values = all$counts,
            labFormat = labelFormat(suffix = " "),
            opacity = 0.85, title = "How many point were counted in each region", position = "topright")

however I dont know is it posible to calculate number of point in each regions having just coordinances?

Answer 1

If you transform the points and France polygons to sf objects, you can use st_intersects() to count the number of points in each polygon.

Note that I updated your sample points so that each quintile break in cuts is unique. With your original data, the first three quintiles were just 0 so the coloring in the leaflet map didn't work.

new sample data

library(leaflet)
library(tidyverse)
set.seed(101)

URL2 <- "https://biogeo.ucdavis.edu/data/gadm3.6/Rsp/gadm36_FRA_2_sp.rds"
data2 <- readRDS(url(URL2))

random_long_lat <- 
  data.frame(
    long = sample(runif(1000, min = -2.5, max = 5.99), replace = F),
    lat = sample(runif(1000, min = 42.69, max = 49.75), replace = F)
  )

make sf objects and count points in polygons

library(sf)

data_sf <- data2 %>% st_as_sf()
random_long_lat_sf <- random_long_lat %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

data_sf_summary <- data_sf %>% 
  mutate(counts = lengths(st_intersects(., random_long_lat_sf)))

define breaks for legend and draw map

cuts <- quantile(data_sf_summary$counts, probs = seq(0, 1, 0.2))
cuts <- colorBin("Greens", domain = data_sf_summary$counts, bins = cuts)

leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data=data_sf_summary, stroke = TRUE, color = "white", weight="", smoothFactor = 0.95, 
              fillOpacity = 0.65, fillColor = ~cuts(data_sf_summary$counts)) %>% 
  addLegend(pal = cuts, 
            values = data_sf_summary$hdp,
            labFormat = labelFormat(suffix = " "),
            opacity = 0.85, title = "How many point were counted in each region", position = "topright")

Also note that tmap package, which lets you switch between static and interactive maps using the same syntax (which resembles ggplot syntax).

same map with tmap:

library(tmap)

tmap_mode("view") # make map interactive
tm_shape(data_sf_summary) + 
  tm_polygons(col = "counts",
              n = 5,
              palette = "Greens",
              title = "How many point were counted in each region")

static map with tmap:

library(tmap)

tmap_mode("plot") # make map static
tm_shape(data_sf_summary) + 
  tm_polygons(col = "counts",
              n = 5,
              palette = "Greens",
              title = "How many point were counted in each region") +
  tm_layout(legend.position = c("right","top"))

---

For multiple countries

First create new sample points that cover Europe:

random_long_lat <- 
  data.frame(
    long = sample(runif(1000, min = -7.5, max = 19.99), replace = F),
    lat = sample(runif(1000, min = 38.69, max = 60.75), replace = F)
  )

all <- rbind(data2, data3, data4, data5)

Then make the sf objects and find the counts of points in every polygon:

all_sf <- all %>% st_as_sf()

random_long_lat_sf <- random_long_lat %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

all_sf_summary <- all_sf %>% 
  mutate(counts = lengths(st_intersects(., random_long_lat_sf)))

qtm(random_long_lat_sf)

Option 1: Choose maps from a list object by name using the NAME_0 column.

tmap_mode("view") # make maps interactive

list_of_maps <- map(unique(all_sf_summary$NAME_0),
                    ~ tm_shape(all_sf_summary %>% 
                                 filter(NAME_0 == .x)) + # filter the data for your country of interest 
                      tm_polygons(col = "counts",
                                  n = 5,
                                  palette = "Greens",
                                  alpha = 0.85,
                                  border.col = NA,
                                  title = "How many point were counted in each region") +
                      tm_layout(legend.position = c("right","top"))) %>% 
  set_names(unique(all_sf_summary$NAME_0))

list_of_maps[['France']]

list_of_maps[['Portugal']]

Option 2: Show all the maps at once

### all maps at once
tm_shape(all_sf_summary) + # filter the data for your country of interest 
  tm_polygons(col = "counts",
              n = 5,
              palette = "Greens",
              alpha = 0.85,
              border.col = NA,
              title = "How many point were counted in each region") +
  tm_layout(legend.position = c("right","top")) +
  tm_facets(by = c("NAME_0"), ncol = 2, showNA = FALSE)