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

(如需轉載,請在顯著位置註明我的微信公衆號stdrei)python

爲何要作地圖投影

簡而言之,地球表面是一個三維的曲面,在曲面上進行測量是很是困難的。不信你拿個地球儀量一下兩點的距離或者計算個夾角試試。將三維的曲面投影到二維平面,這樣咱們學的平面幾何纔有用武之地。git

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

什麼是投影

若是要了解地圖投影,首先須要知道地理座標系投影座標系。簡單的說,地理座標系是參考平面爲橢球面的球面座標,最多見的是以經緯度來量算的球面座標系統。投影座標系是定義在一個二維平面的座標系統,其地圖單位一般爲米。將球面座標轉化爲平面座標的過程便稱爲投影。 所以,投影座標系老是基於地理座標系的(連座標系都要講出身的)。數據庫

由於地球是一個赤道略寬兩極略扁的不規則的梨形球體,其表面是一個不可展平的曲面,因此運用任何數學方法進行這種投影轉換都會產生偏差和變形,爲按照不一樣的需求縮小偏差,就產生了各類投影方法。好比這篇頗有意思的文章你所不知的有趣投影方法。也能夠看文末的幾幅地球不一樣投影的圖。微信

EPSG Code

若是常常跟地圖打交道,你必定會被各類花式投影變換弄的眼花繚亂。EPSGEuropean 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。ide

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.iospa

pyproj小試牛刀

pyproj是PROJ4的python接口封裝,直接看一個官網的例子吧。直接利用epsg code來定義投影參數。.net

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)

輸出爲code

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...

在不一樣投影下的這個世界。。。

相關文章
相關標籤/搜索