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

01-29

(如需轉載，請在顯著位置註明個人微信公眾號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,transformnn# The Proj class can convert from geographic (longitude,latitude)n# to native map projection (x,y) coordinates and vice versa, n# or from one map projection coordinate system directly to another.nnp1 = Proj(init=epsg:26915)np2 = Proj(init=epsg:26715)nx1, y1 = p1(-92.199881,38.56694)nx2, y2 = transform(p1,p2,x1,y1)nprint %9.3f %11.3f % (x1,y1)nprint %9.3f %11.3f % (x2,y2)n

輸出為

569704.566 4269024.671n569722.342 4268814.027n

基於geopandas的矢量地圖投影

import shapely, geopandas, fionanimport seaborn as snsnfrom fiona.crs import from_epsg,from_stringnn# Datanshp = E:NationalGISdataProvince.shpnnshp_df = geopandas.GeoDataFrame.from_file(shp)n# #IndexError報錯的話，用arcgis將shapefile文件重新導出一遍試試nshp_df.head()n

# 根據當前的蘭伯特投影繪製nshp_df.plot(column="GDP_1994",colormap=Set1)n

# 轉換到經緯度坐標nshp_df_wgs84 = shp_df.to_crs(from_epsg(4326))nshp_df_wgs84.plot(column="GDP_1994",colormap=Set1)n

# 國家2000坐標系n# EPSG:4508 CGCS2000 / Gauss-Kruger CM 111Enshp_df_4508 = shp_df.to_crs(from_epsg(4508))nshp_df_4508.plot(column="GDP_1994",colormap=Set1)n

# 除了直接用ESPG code，也可以自己定義投影參數nnESRI_54024 = """n+proj=bonne +lon_0=0 +lat_1=60 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs nnshp_df_3408 = shp_df.to_crs(from_string(ESRI_54024))nshp_df_3408.plot(column="GDP_1994",colormap=Set1)n

代碼和矢量地圖數據下載

百度網盤下載地址： http://pan.baidu.com/s/1c1BExH2

密碼請在關注微信公眾號stdrei後，輸入projection直接獲取。

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

鏈接 One world, many faces: A brief look at map projections - Views of the World
在不同投影下的這個世界。。。