左手用R右手Python系列12——空间数据可视化与数据地图
以前我一直觉得Python的绘图工具与R语言ggplot2比起来，不够优雅，这也是我一直坚定的选择使用R+ggplot2深入的学习数据可视化的原因，ggplot2在坐标系的整合与兼容性和扩展性上确实技高一筹，所以ggplot2成了可视化的巨无霸，成了可视化界的微信，不仅自身生态日趋完善，而且还有众多的开发者为其开发辅助功能包（你可以理解为依附于微信的小程序）。
最近偶然在学习Python可视化的过程中，了解到了geopandas，确实第一眼看着很眼熟，或许你第一眼就能把它与pandas联系起来。的确，它跟pandas有着千丝万缕的联系，并且继承了pandas的诸多高频函数。而geo是什么鬼呢？
geo是地理信息系统的简写，geopandas是Python中用于处理空间地理信息数据的后起之秀（为什么是后起之秀呢，因为有个叫basemap的包，据说很难用，我还没有深入了解过呢）。
今天要讲解的主角是R语言中的sf包和Python中的geopandas库。
为什么今天把geopandas和R语言空间数据可视化写在一起，因为他们很巧合的用到了相同的地理信息处理技术，无论是数据源的支持上、还是空间数据的结构存储还是投影设置上都是如此。
我以前写过大量ggplot2空间数据可视化的文章，但是那些大多是基于shp数据源，和ggplot2中的geom_ploygon或者gemo_map函数来制作的。
shp数据源导入所依赖的maptools包已经快要被遗弃了，而geom_ploygon、gemo_map函数所支持的数据结构复杂难懂，无论是对于初学者还是老手们都是一种挑战，需要做大量的数据合并、转换、匹配，在前期的数据处理上花费的时间和代码量已经远远超过了可视化的代码量。
好在新技术总是不断地出现，数据源上json格式的数据为我们提供了更为便捷、高效、低廉的空间数据信息，而sf包则可以使用直观易懂的Simple Features数据结构来从新规整地图数据源，使得过去需要分别准备地理边界属性信息和地理边界经纬点信息来呈现地理空间信息数据结构。
巧合的是，python中的geopandas用了同样的 技术来简化空间数据可视化的复杂度，其核心理念也是通过压缩单个地理多边形为一个Simple Features，使得所有的地理多边形与其属性信息严格对齐，行政一个呈现友好的带有地理信息数据的数据框。
也许以上描述过于抽象，因为涉及到到的内容比较深入，我实在是不知道该如何把这些内容将的通俗易懂，接下来会使用图片辅助演示。
R中制作地图传统的方式是使用geom_ploygon+maptools+shp数据
library(ggplot2)
library(plyr)
library(maptools)
#数据导入：
china_map<-readShapePoly("D:/R/rstudy/CHN_adm/bou2_4p.shp")
Warning message:
#已经开始提出警告了！
(替代方案，使用rgdal中的readORG函数或者sf包中的st_read函数)
use rgdal::readOGR or sf::st_read
china_map1<-fortify(china_map)
#从SP（空间数据对象）中剥离地理多边形边界点信息和多边形属性信息
x<-china_map@data
xs<-data.frame(id=row.names(x),x)
china_map_data <- join(china_map1, xs, type = "full")
#导入业务数据
province_city <- read.csv("D:/R/rstudy/Province/chinaprovincecity.csv")
mydata <- read.csv("D:/R/rstudy/Province/geshengzhibiao.csv")
china_data <- join(china_map_data, mydata, type="full")
#可视化代码：
ggplot(china_data,aes(long,lat))+
geom_polygon(aes(group=group),fill="white",colour="grey60")+
geom_point(data =province_city,aes(x = jd,y = wd),colour="red")+
coord_map("polyconic")+
theme_void()
过程何其辛苦！
为什么使用maptools+geom_ploygon技术组合这么辛苦呢，问题出在数据源上，如果你想要详细了解maptools导入的空间信息结果以及goem_ploygon根据什么规则映射地图信息，请看这一篇。
我能告诉你的是，geom_ploygon制作地图的时候，剥离了地理信息边界点数据和多边形属性信息，所以你需要同时兼顾、处理两个包含空间信息的数据框，如果是对不同区域进行等值线映射，你还需要对这两个数据框进行合并操作。
而sf包则使用了新的、更为优雅简洁的空间信息呈现技术——Simple Features
以上便是使用shp+maptools+geom_ploygon技术的核心数据结构概况，接下来我们会跟大家讲解新技术组合下所支持的空间数据结构。
sf包则也是同时支持shp数据源和json数据源
library("sf")
library("ggplot2")
china_map<-st_read("D:/R/rstudy/CHN_adm/bou2_4p.shp",stringsAsFactors=FALSE,quiet=TRUE)
Encoding(china_map$NAME)<-"GBK"
china_map<-st_read("D:/R/mapdata/State/china.geojson",stringsAsFactors=FALSE,quiet=TRUE)
####转换编码
st_crs(china_map)$epsg<-4267
st_crs(china_map)$proj4string<-"+proj=longlat +datum=NAD27 +no_defs"
china_map<-st_transform(china_map,3395)
#####合并数据
mydata<-data.frame(NAME=unique(china_map$NAME),value=runif(34,1,100))
china_map<-merge(china_map,mydata,by="NAME")
ggplot() +
geom_sf(data=china_map,aes(fill=value))
所以使用sf提供的新技术，制作数据地图通常仅需以上几步。任务量大大缩减。
而Python的geopandas包则也提供了相同的空间数据结构处理技术。
import pandas as pd
import numpy as np
import geopandas as gp
import matplotlib.pyplot as plt
province_city = pd.read_csv("D:/R/rstudy/Province/chinaprovincecity.csv", encoding = 'gb18030')
china_map=gp.GeoDataFrame.from_file("D:/R/rstudy/CHN_adm/bou2_4p.shp", encoding = 'gb18030')
china_map=gp.GeoDataFrame.from_file("D:/R/mapdata/State/china.geojson", encoding = 'gb18030')
geopandas包同时支持导入shp素材和json素材，导入之后得数据结构与R语言中得sf导入之后得结构是一致得，地理多边形边界点信息都被压缩成了一个非常整齐的列表存储，列表内每一个单独的子项目都代表着一个多边形。
实际上导入之后，你可以看到它的结构是一种特殊的带有地理信息列的数据框。
geopandas.geodataframe.GeoDataFrame
这种格式数据框继承了大多数pandas普通数据框的函数及属性，可以直接针对其使用plot函数绘图。
china_map.plot(figsize=(20,12))
其内部数据结构与sf包的数据结构如出一辙。
我们可以给其指定一个数值变量，使得映射出来的地图各区块根据数值单独填色。
china_map.plot(column="AREA",figsize=(20,12),cmap="Greens")
如果你想要在此图层上添加另外一个散点图层，则需设定两个具有同样投影信息的GeoDataFrame对象。
使用刚才导入的province_city数据，将其合并进china_map中去。
china_map=china_map.merge(province_city,left_on='NAME', right_on='province', how='left')
from shapely.geometry import Point
china_map["center"]=gp.GeoSeries([Point(x, y) for x, y in zip(china_map["jd"], china_map["wd"])])
china_map_ploygon=china_map[["AREA","NAME","geometry","class"]]
china_map_ploygon=china_map_ploygon
china_map_point=china_map[["AREA","NAME","center","class"]]
china_map_point=china_map_point.set_geometry('center')
###设置相同的投影：
china_map_ploygon.crs={'init': 'epsg:3395'}
china_map_point.crs ={'init': 'epsg:3395'}
###地图可视化过程：
base=china_map_ploygon.plot(column="AREA", edgecolor='black',figsize=(20,12),cmap="Greens")
china_map_point.plot(ax=base,marker='o',color='red',markersize=5)
最后让我们再次看一下R语言中的sf数据对象和Pyhton中的geodatafame对象的对比。
在线课程请点击文末阅读原文：
数据源请移步GitHub：
https://github.com/ljtyduyu/DataWarehouse