r - 从坐标点创建分区统计图

Question

我有一个由具有特定地理坐标(纬度和经度)的多个数据点组成的数据框。我希望创建一个等值区域样式的世界地图，其中地理区域根据落在该区域边界内的数据点数量进行着色。

有没有一种简单的方法来完成我在 R 中尝试做的事情，最好使用“maps”包的世界地图和“ggplot2” map 绘制函数？

这是我所拥有的最低限度可重现的结果:

library(ggplot2)
library(maps)

data <- data.frame(lat = 40.730610, lon = -73.935242)

ggplot() + 
  geom_polygon(data = map_data("world"), aes(x = long, y = lat, group = group, fill = group)) +
  coord_fixed(1.3)

我注意到绘图项函数上的 fill 参数可用于创建分区统计图效果。此处，geom_polygon() 函数的 aes() 函数上的 fill 参数用于创建一个分区统计图，其中每个组都采用颜色编码不同。

Answer 1

有很多方法可以完成此任务。总体思路是将点数据和多边形数据转换为空间对象。之后，计算有多少个点落在该多边形内。我知道我们可以使用 sp 包来完成此操作，该包在 R 社区中广泛且众所周知，但我决定使用 sf 包，因为 sf 将成为 R 中空间对象的下一代标准 ( https://cran.r-project.org/web/packages/sf/index.html )。了解 sf 的用法和功能可能会有所帮助。

首先，OP 提供了一个示例点，但我决定添加更多点，以便我们可以了解如何计算点并聚合数据。为此，我使用 ggmap pakcage 对我选择作为示例的一些城市进行地理编码。

# Load package
library(tidyverse)
library(ggmap)
library(maps)
library(maptools)
library(sf)

# Create point data as a data frame
point_data <- data.frame(lat = 40.730610, lon = -73.935242)

# Geocode a series of cities
city <- c("Detroit", "Seattle", "Toranto", "Denver", "Mexico City", "Paris", "New Orleans",
          "Tokyo", "Osaka", "Beijing", "Canberra", "New York", "Istanbul", "New Delhi",
          "London", "Taipei", "Seoul", "Manila", "Bangkok", "Lagos", "Chicago", "Shanghai")
point_data2 <- geocode(city)

# Combine OP's example and the geocoding result
point_data3 <- bind_rows(point_data, point_data2)

接下来，我将 point_data3 数据帧转换为 sf 对象。我还将使用 maps 包获取世界的多边形数据，并将其转换为 sf 对象。

# Convert to simple feature object
point_sf <- st_as_sf(point_data3, coords = c("lon", "lat"), crs = 4326)

# Get world map data
worldmap <- maps::map("world", fill = TRUE, plot = FALSE)

# Convert world to sp class
IDs <- sapply(strsplit(worldmap$names, ":"), "[", 1L)
world_sp <- map2SpatialPolygons(worldmap, IDs = IDs, 
                                proj4string = CRS("+proj=longlat +datum=WGS84"))

# Convert world_sp to simple feature object
world_sf <- st_as_sf(world_sp)

# Add country ID
world_sf <- world_sf %>%
  mutate(region = map_chr(1:length(world_sp@polygons), function(i){
    world_sp@polygons[[i]]@ID
  }))

现在 point_sf 和 world_sf 都是 sf 对象。我们可以使用 st_within 函数来检查哪些点位于哪些多边形内。

# Use st_within
result <- st_within(point_sf, world_sf, sparse = FALSE)

# Calculate the total count of each polygon
# Store the result as a new column "Count" in world_sf
world_sf <- world_sf %>%
  mutate(Count = apply(result, 2, sum))

总计数信息位于world_sf的Count列中。我们可以像OP一样使用map_data函数获取世界数据框。然后我们可以合并 world_data 和 world_df。

# Convert world_sf to a data frame world_df 
world_df <- world_sf
st_geometry(world_df) <- NULL

# Get world data frame
world_data <- map_data("world")

# Merge world_data and world_df
world_data2 <- world_data %>%
  left_join(world_df, by = c("region"))

现在我们准备绘制数据了。以下代码与 OP 的 ggplot 代码相同，只是输入数据现在是 world_data2 和 fill = Count。

ggplot() + 
  geom_polygon(data = world_data2, aes(x = long, y = lat, group = group, fill = Count)) +
  coord_fixed(1.3)