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

(如需转载，请在显著位置注明我的微信公众号stdrei)

为何要作地图投影

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

什么是投影

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

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

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

代码和矢量地图数据下载

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

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

连接 http://www.viewsoftheworld.ne...

在不一样投影下的这个世界。。。