用 Python 做地图投影 - 多面孔的世界

(如需转载，请联系作者)

为什么要做地图投影

简而言之，地球表面是一个三维的曲面，在曲面上进行测量是非常困难的。不信你拿个地球仪量一下两点的距离或者计算个夹角试试。将三维的曲面投影到二维平面，这样我们学的平面几何才有用武之地。

youtube - Adventures with GIS: an introduction to IPython and the joy of play

什么是投影

如果要了解地图投影，首先需要知道地理坐标系和投影坐标系。简单的说，地理坐标系是参考平面为椭球面的球面坐标，最常见的是以经纬度来量算的球面坐标系统。投影坐标系是定义在一个二维平面的坐标系统，其地图单位通常为米。将球面坐标转化为平面坐标的过程便称为投影。因此，投影坐标系总是基于地理坐标系的（连坐标系都要讲出身的）。

因为地球是一个赤道略宽两极略扁的不规则的梨形球体，其表面是一个不可展平的曲面，所以运用任何数学方法进行这种投影转换都会产生误差和变形，为按照不同的需求缩小误差，就产生了各种投影方法。比如这篇很有意思的文章你所不知的有趣投影方法。也可以看文末的几幅地球不同投影的图。

EPSG Code

如果经常跟地图打交道，你一定会被各种花式投影变换弄的眼花缭乱。EPSG：European Petroleum Survey Group (EPSG)成立于1986年，并在2005年重组为OGP(Internation Association of Oil & Gas Producers)，它负责维护并发布坐标参照系统的数据集参数，以及坐标转换描述，该数据集被广泛接受并使用，通过一个Web发布平台进行分发，同时提供了微软Acess数据库的存储文件，通过SQL 脚本文件，MySQL, Oracle 和PostgreSQL等数据库也可使用。目前已有的椭球体，投影坐标系等不同组合都对应着不同的ID号，这个号在EPSG中被称为EPSG code，它代表特定的椭球体、单位、地理坐标系或投影坐标系等信息。

开源的QGIS软件中就直接采用了EPSG。

The CRSs available in QGIS are based on those defined by the European Petroleum Search Group (EPSG) and the Institut Geographique National de France (IGNF) and are largely abstracted from the spatial reference tables used in GDAL. EPSG identifiers are present in the database and can be used to specify a CRS in QGIS.

一个查询EPSG code的网站 epsg.io

pyproj小试牛刀

pyproj是PROJ4的python接口封装，直接看一个官网的例子吧。直接利用epsg code来定义投影参数。

from pyproj import Proj,transform

# The Proj class can convert from geographic (longitude,latitude)
# to native map projection (x,y) coordinates and vice versa, 
# or from one map projection coordinate system directly to another.

p1 = Proj(init='epsg:26915')
p2 = Proj(init='epsg:26715')
x1, y1 = p1(-92.199881,38.56694)
x2, y2 = transform(p1,p2,x1,y1)
print '%9.3f %11.3f' % (x1,y1)
print '%9.3f %11.3f' % (x2,y2)

输出为

569704.566 4269024.671
569722.342 4268814.027

基于geopandas的矢量地图投影

import shapely, geopandas, fiona
import seaborn as sns
from fiona.crs import from_epsg,from_string

# Data
shp = 'E:\NationalGISdata\Province.shp'

shp_df = geopandas.GeoDataFrame.from_file(shp)
# #IndexError报错的话，用arcgis将shapefile文件重新导出一遍试试
shp_df.head()

# 根据当前的兰伯特投影绘制
shp_df.plot(column="GDP_1994",colormap='Set1')

# 转换到经纬度坐标
shp_df_wgs84 = shp_df.to_crs(from_epsg(4326))
shp_df_wgs84.plot(column="GDP_1994",colormap='Set1')

# 国家2000坐标系
# EPSG:4508  CGCS2000 / Gauss-Kruger CM 111E
shp_df_4508 = shp_df.to_crs(from_epsg(4508))
shp_df_4508.plot(column="GDP_1994",colormap='Set1')

# 除了直接用ESPG code，也可以自己定义投影参数

ESRI_54024 = """
+proj=bonne +lon_0=0 +lat_1=60 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs 

shp_df_3408 = shp_df.to_crs(from_string(ESRI_54024))
shp_df_3408.plot(column="GDP_1994",colormap='Set1')

代码和矢量地图数据下载

百度网盘下载地址： pan.baidu.com/s/1c1BExH2
密码请在关注微信公众号stdrei后，输入projection直接获取。

拓展：同一个世界，不同的面孔

链接 www.viewsoftheworld.net/?p=752
在不同投影下的这个世界。。。