方法简介
我们以加拿大的13个省级行政单位的GDP增长率为例，绘制GDP的空间分布填色图。大致思路如下：我们先把省份 – 地理坐标对应关系找到。接着，准备好2015-16年的各省GDP年增长率，随后将其与地图中的每个省一一匹配，做成GDP – 省份的对应关系。这样GDP就可以经由省份对应到地理位置上。最后根据增速快慢赋予每个省不同的颜色。
准备工作
-
我们需要先下载加拿大的省级行政区的底图。Global Administrative Areas (http://www.gadm.org/country) 提供了中国的行政区划图，我们可以选择“Canada”和“Shapefile”，下载底图 (CAN_adm_shp.zip)。
网站上也提供了专门针对R的底图——R (SpatialPolygonsDataFrame)，但是这一文件格式只适用于R，而Shapefile也可以利用ArcGIS处理，使用范围更广，因此本推送以Shapefile格式为例进行说明。 -
下载各个省份的GDP年度数据，计算增长率。这一数据可以从加拿大国家统计局网站 (http://www.statcan.gc.ca/pub/13-016-x/13-016-x2017001-eng.htm) 上下载，保存为excel格式后计算增长率。
读取底图 (省份 – 地理坐标)
首次操作，需要先安装rgdal和ggplot2库。
> install.packages("rgdal")
> install.packages("ggplot2")
我们使用readORG命令读取刚刚下载的Shapefile底图。我们把CAN_adm_shp.zip解压至R的工作目录下，读取CAN_adm1.shp。
若不清楚工作目录，可以使用getwd()命令获取。
> library(rgdal)
> shp <- readOGR("CAN_adm1.shp")
操作成功后，R会输出该底图的一些基础信息
OGR data source with driver: ESRI Shapefile Source: "USA_adm1.shp", layer: "CAN_adm1"with 15 features It has 12 fields Integer64 fields read as strings: ID_0 ID_1 CCN_1
随后，我们将读入的shp中的data另存出来，用以加入GDP增长率数据。
> data <- shp@data
同样，我们把shp中的几何映射地图 (即坐标 – 地点对应关系) 读出
> library(ggplot2) > map <- fortify(shp)
将GDP增速加入到数据中 (GDP – 省份)
读取加拿大13个省2015到16年的GDP增长率。同样，这份GDP增速的excel也需要提前放在工作目录下。
> install.packages("readxl") #如果已经安装过该包，无需重复安装
> library(readxl)
> province <- read_excel("growth.xlsx")
然后，把growth里的growth列增加到先前的data里，创建GDP – 省份对应关系。
> growthnew <- rbind(growth[1:8,], growth[8,], growth[8:nrow(state),]) > data$growth <- growthnew$growth > data$ID_1 <- c(0:13) #因为map中，省份ID从0而开始，而data中的ID从1开始，因此我们必须在这里将两者的ID统一从0开始
其中第一步是为了处理努纳武特省的问题，因为其行政区划并不是一个连续的封闭图形，因此在shp文件中，同一个省对应了3个ID。
我们第一步的目的就是使得同一个省的GDP对应到这3个ID上。类似问题有很多，比如美国阿拉斯加州、中国的河北省，都可能需要进行类似的操作。
将数据和底图结合 (GDP – 省份 – 地理坐标)
> new <- merge(map[,c(-4,-5)], data[,c("ID_1","NAME_1","growth")],by.x="id",by.y="ID_1")
这样，我们就得到了new，其中GDP增速已经对应到地理坐标上去了。
最后我们用ggplot2绘图：
> library(ggthemes)
> ggplot() +
+ geom_polygon(data = new, aes(x = long, y = lat, group = group, fill = growth), col= "grey95") +
+ scale_fill_gradient(low = "#995A7B", mid = "#9ED0E6", high = "#EAA1AC") +
+ coord_map("polyconic") +
+ theme_map()
