(如需轉(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)單的說(shuō)，地理坐標(biāo)系是參考平面為橢球面的球面坐標(biāo)，最常見(jiàn)的是以經(jīng)緯度來(lái)量算的球面坐標(biāo)系統(tǒng)。投影坐標(biāo)系是定義在一個(gè)二維平面的坐標(biāo)系統(tǒng)，其地圖單位通常為米。將球面坐標(biāo)轉(zhuǎn)化為平面坐標(biāo)的過(guò)程便稱(chēng)為投影。 因此，投影坐標(biāo)系總是基于地理坐標(biāo)系的（連坐標(biāo)系都要講出身的）。
因?yàn)榈厍蚴且粋€(gè)赤道略寬兩極略扁的不規(guī)則的梨形球體，其表面是一個(gè)不可展平的曲面，所以運(yùn)用任何數(shù)學(xué)方法進(jìn)行這種投影轉(zhuǎn)換都會(huì)產(chǎn)生誤差和變形，為按照不同的需求縮小誤差，就產(chǎn)生了各種投影方法。比如這篇很有意思的文章你所不知的有趣投影方法。也可以看文末的幾幅地球不同投影的圖。
EPSG Code
如果經(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ù)集被廣泛接受并使用，通過(guò)一個(gè)Web發(fā)布平臺(tái)進(jìn)行分發(fā)，同時(shí)提供了微軟Acess數(shù)據(jù)庫(kù)的存儲(chǔ)文件，通過(guò)SQL 腳本文件，MySQL, Oracle 和PostgreSQL等數(shù)據(jù)庫(kù)也可使用。目前已有的橢球體，投影坐標(biāo)系等不同組合都對(duì)應(yīng)著不同的ID號(hào)，這個(gè)號(hào)在EPSG中被稱(chēng)為EPSG code，它代表特定的橢球體、單位、地理坐標(biāo)系或投影坐標(biāo)系等信息。
開(kāi)源的QGIS軟件中就直接采用了EPSG。
The CRSs available in QGIS are based>pyproj小試牛刀
pyproj是PROJ4的python接口封裝，直接看一個(gè)官網(wǎng)的例子吧。直接利用epsg code來(lái)定義投影參數(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>'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報(bào)錯(cuò)的話，用arcgis將shapefile文件重新導(dǎo)出一遍試試shp_df.head()
# 根據(jù)當(dāng)前的蘭伯特投影繪制shp_df.plot(column="GDP_1994",colormap='Set1')
