从坐标点创建Choropleth贴图 [英] Create choropleth map from coordinate points

问题描述

我有一个数据框，其中包含多个具有特定地理坐标(纬度和经度)的数据点.我正在尝试创建一个choroplethth风格的世界地图，其中根据有多少数据点落在该区域的边界内来对地理区域进行着色.

I have a data frame consisting of multiple data points with specific geocoordinates (latitude and longitude). I'm looking to create a choropleth-style world map where geographical regions are shaded according to how many data points fall within the boundaries of the region.

是否有一种简单的方法来完成我要在R中尝试做的事情，最好使用"maps"程序包的世界地图和"ggplot2"地图绘图功能?

Is there a simple way to accomplish what I'm trying to do in R, preferably using the "maps" package's world map and the "ggplot2" map-plotting functions?

以下是我所拥有的内容的可复制性最低的结果:

Here is a minimally reproducible result of what I have:

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参数可用于创建一个Choropleth效果.在这里，geom_polygon()函数的aes()函数上的fill参数用于创建一个Choropleth，其中每个组的颜色编码都不同.

I've noticed that the fill parameter on plot item functions can be used to create a choropleth effect. Here, the fill parameter on the aes() function of the geom_polygon() function is used to create a choropleth where each group is color coded differently.

推荐答案

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

There are many ways to achieve this task. The general idea is to convert both the point data and polygon data to spatial objects. After that, count how many points fall within that polygon. I know we can do this using the sp package, which is widespread and well-known in the R community, but I decided to use the sf package because sf would be the next generation standard of spatial objects in R (https://cran.r-project.org/web/packages/sf/index.html). Knowing the usage and functionality of sf will probably be beneficial.

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

First, the OP provided an example point, but I decided to add more points so that we can see how to count the points and aggregate the data. To do so, I used the ggmap pakcage to geocode some cities that I selected as an example.

# 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对象.

Next, I converted the point_data3 data frame to the sf object. I will also get the polygon data of the world using the maps package and convert it to an sf object.

# 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函数检查哪些点在哪些多边形内.

Now both point_sf and world_sf are sf objects. We can use the st_within function to examine which points are within which polygons.

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

The total count information is in the Count column of world_sf. We can get the world data frame as the OP did using the map_data function. We can then merge world_data and 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"))

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

Now we are ready to plot the data. The following code is the same as the OP's ggplot code except that the input data is now world_data2 and fill = Count.

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

这篇关于从坐标点创建Choropleth贴图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！