用Python做地图投影

什么是地图投影

地图投影是利用一定数学法则把地球表面的经、纬线转换到平面上的理论和方法。

由于地球是一个赤道略宽两极略扁的不规则的梨形球体，故其表面是一个不可展平的曲面，所以运用任何数学方法进行这种转换都会产生误差和变形，为按照不同的需求缩小误差，就产生了各种投影方式。

比如这篇很有意思的文章你所不知的有趣投影方法。

为什么要进行地图投影

1. 地理坐标为球面坐标, 不方便进行距离, 方位, 面积等参数的量算.
2. 地球椭球体为不可展曲面.
3. 地图为平面, 符合视觉心理,并进行距离, 方位, 面积等量算和各种空间分析.

常用的投影方法

​ 1、几何透视法：几何透视法是利用透视的关系，将地球体面上的点投影到投影面上的一种投影方法。如假设地球按比例缩小成一个透明的地球仪般的球体，在其球心或球面、球外安置一个光源，将球面上的经纬线投影到球外的一个投影平面上，即将球面经纬线转换成了平面上的经纬线。

​ 2、数学解析法：数学解析法是在球面与投影面之间建立点与点的函数关系，通过数学的方法确定经纬线交点位置的一种投影方法。

几何透视法是一种比较原始的投影方法，有很大的局限性，难于纠正投影变形，精度较低，当前绝大多数地图投影都采用数学解析法。大多数的数学解析法往往是在透视投影的基础上，发展建立球面与投影面之间点与点的函数关系的，因此两种投影方法有一定联系。

地图投影主要分为3个过程

​ 1. 确定地球椭球体（Spheroid/Ellipsoid），需要长半轴、短半轴、曲率三个参数。（模拟地球的形状）

​ 2. 若要逼近某特定地区，则需要大地基准（Geodetic Datum）。（椭球体的原点the position of the origin、方向 the orientation、缩放比例 the scale等）我国现在采取西安1980坐标系基准点，同时也有国家1985高程基准。

​ 3. 关于如何投影，这就涉及高斯-克吕格投影等诸多投影方式的存在，等角等积等距离，方位圆柱圆锥等分类。

由于世界各地区投影类型的不同，因此在叠加、复合不同来源空间数据时，必需首先进行投影转换、配准等设置。

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

更多图示:

(文章参考资料:微信公众号stdrei)