计算区域 [R] 的多边形中有多少个点

Question

我想问当我们有国家及其地区的点和多边形的经度和纬度变量时，如何计算某个地区的点数。

我在下面提供了示例：我想计算哪些区域中有多少个点（当没有点时包括零），然后将此变量添加到data2@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))


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包，它允许您使用相同的语法（类似于 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")

静态地图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)