从坐标点创建分区统计图

Question

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

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

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

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绘图项函数上的参数可用于创建分区统计图效果。在这里，函数的fill参数用于创建分区统计图，其中每个组的颜色编码不同。aes()geom_polygon()

Answer 1

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

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

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

总计数信息在Count列中world_sf。我们可以像 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代码相同，只是输入数据是 nowworld_data2和fill = Count。

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