Подсчет, сколько точек в многоугольнике области [R]

Question

Я хотел бы спросить, как рассчитать количество точек в некотором регионе, когда у нас есть переменные долготы и широты точки и многоугольника страны и ее регионов.

Я привел пример ниже: Я хотел бы вычислить, сколько точек находится в каких регионах (включая ноль, когда точки нет) и затем добавить эти переменные в объект data2@data, чтобы можно было использовать count меры для заполнения каждого полигона областей.

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  

Я хотел бы получить результат, чтобы можно было вычислить что-то вроде этого:

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

однако я не знаю, можно ли подсчитать количество точек в каждой области, имеющей только координаты?

Answer 1

Если вы преобразуете точки и многоугольники Франции в sf объекты, вы можете использовать st_intersects() для подсчета количества точек в каждом многоугольнике.

Обратите внимание, что я обновил ваши выборочные точки так, чтобы каждый квинтиль перерыв в cuts уникален. Исходя из ваших исходных данных, первые три квинтиля были равны 0, поэтому раскраска на листовой карте не сработала.

новые образцы данных

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

сделать sf объекты и считать точки в многоугольниках

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

определить разрывы для легенды и нарисовать карту

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

Также обратите внимание, что tmap пакет, который позволяет переключаться между stati c и интерактивными картами, используя один и тот же синтаксис (который напоминает синтаксис ggplot).

та же карта с 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")

stati c карта с 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"))

---

Для нескольких стран

Сначала создайте новые точки отбора проб, которые охватывают Европу:

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)

Затем создайте объекты sf и найдите количество точек в каждом многоугольнике:

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)

Вариант 1: Выберите карты из объекта списка по имени, используя NAME_0 столбец.

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']]

Опция 2: Показать все карты одновременно

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