摘要：如需轉(zhuǎn)載，請(qǐng)?jiān)陲@著位置注明個(gè)人微信公眾號(hào)為什么要做地圖投影簡(jiǎn)而言之，地球表面是一個(gè)三維的曲面，在曲面上進(jìn)行測(cè)量是非常困難的。拓展同一個(gè)世界，不同的面孔鏈接在不同投影下的這個(gè)世界。。。

(如需轉(zhuǎn)載，請(qǐng)?jiān)陲@著位置注明個(gè)人微信公眾號(hào)stdrei)

簡(jiǎn)而言之，地球表面是一個(gè)三維的曲面，在曲面上進(jìn)行測(cè)量是非常困難的。不信你拿個(gè)地球儀量一下兩點(diǎn)的距離或者計(jì)算個(gè)夾角試試。將三維的曲面投影到二維平面，這樣我們學(xué)的平面幾何才有用武之地。

如果要了解地圖投影，首先需要知道地理坐標(biāo)系和投影坐標(biāo)系。簡(jiǎn)單的說，地理坐標(biāo)系是參考平面為橢球面的球面坐標(biāo)，最常見的是以經(jīng)緯度來量算的球面坐標(biāo)系統(tǒng)。投影坐標(biāo)系是定義在一個(gè)二維平面的坐標(biāo)系統(tǒng)，其地圖單位通常為米。將球面坐標(biāo)轉(zhuǎn)化為平面坐標(biāo)的過程便稱為投影。?因此，投影坐標(biāo)系總是基于地理坐標(biāo)系的（連坐標(biāo)系都要講出身的）。

因?yàn)榈厍蚴且粋€(gè)赤道略寬兩極略扁的不規(guī)則的梨形球體，其表面是一個(gè)不可展平的曲面，所以運(yùn)用任何數(shù)學(xué)方法進(jìn)行這種投影轉(zhuǎn)換都會(huì)產(chǎn)生誤差和變形，為按照不同的需求縮小誤差，就產(chǎn)生了各種投影方法。比如這篇很有意思的文章你所不知的有趣投影方法。也可以看文末的幾幅地球不同投影的圖。

如果經(jīng)常跟地圖打交道，你一定會(huì)被各種花式投影變換弄的眼花繚亂。EPSG：European Petroleum Survey Group (EPSG)成立于1986年，并在2005年重組為OGP(Internation Association of Oil & Gas Producers)，它負(fù)責(zé)維護(hù)并發(fā)布坐標(biāo)參照系統(tǒng)的數(shù)據(jù)集參數(shù)，以及坐標(biāo)轉(zhuǎn)換描述，該數(shù)據(jù)集被廣泛接受并使用，通過一個(gè)Web發(fā)布平臺(tái)進(jìn)行分發(fā)，同時(shí)提供了微軟Acess數(shù)據(jù)庫的存儲(chǔ)文件，通過SQL 腳本文件，MySQL, Oracle 和PostgreSQL等數(shù)據(jù)庫也可使用。目前已有的橢球體，投影坐標(biāo)系等不同組合都對(duì)應(yīng)著不同的ID號(hào)，這個(gè)號(hào)在EPSG中被稱為EPSG code，它代表特定的橢球體、單位、地理坐標(biāo)系或投影坐標(biāo)系等信息。

開源的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.

一個(gè)查詢EPSG code的網(wǎng)站 epsg.io

pyproj是PROJ4的python接口封裝，直接看一個(gè)官網(wǎng)的例子吧。直接利用epsg code來定義投影參數(shù)。

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

import shapely, geopandas, fiona
import seaborn as sns
from fiona.crs import from_epsg,from_string

# Data
shp = "E:NationalGISdataProvince.shp"

shp_df = geopandas.GeoDataFrame.from_file(shp)
# #IndexError報(bào)錯(cuò)的話，用arcgis將shapefile文件重新導(dǎo)出一遍試試
shp_df.head()

# 根據(jù)當(dāng)前的蘭伯特投影繪制
shp_df.plot(column="GDP_1994",colormap="Set1")

# 轉(zhuǎn)換到經(jīng)緯度坐標(biāo)
shp_df_wgs84 = shp_df.to_crs(from_epsg(4326))
shp_df_wgs84.plot(column="GDP_1994",colormap="Set1")

# 國(guó)家2000坐標(biāo)系
# 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，也可以自己定義投影參數(shù)

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

百度網(wǎng)盤下載地址： http://pan.baidu.com/s/1c1BExH2
密碼請(qǐng)?jiān)陉P(guān)注微信公眾號(hào)stdrei后，輸入projection直接獲取。

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

在不同投影下的這個(gè)世界。。。