用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>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>'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.671569722.342 4268814.027

基于geopandas的矢量地图投影

import shapely, geopandas, fionaimport seaborn as snsfrom fiona.crs import from_epsg,from_string# Datashp = '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')