(如需轉(zhuǎn)載圈膏,請在顯著位置注明個人微信公眾號stdrei)
為什么要做地圖投影
簡而言之狈癞,地球表面是一個三維的曲面洁灵,在曲面上進行測量是非常困難的暮屡。不信你拿個地球儀量一下兩點的距離或者計算個夾角試試。將三維的曲面投影到二維平面温眉,這樣我們學的平面幾何才有用武之地缸匪。
什么是投影
如果要了解地圖投影,首先需要知道地理坐標系和投影坐標系类溢。簡單的說凌蔬,地理坐標系是參考平面為橢球面的球面坐標,最常見的是以經(jīng)緯度來量算的球面坐標系統(tǒng)闯冷。投影坐標系是定義在一個二維平面的坐標系統(tǒng)砂心,其地圖單位通常為米。將球面坐標轉(zhuǎn)化為平面坐標的過程便稱為投影蛇耀。 因此辩诞,投影坐標系總是基于地理坐標系的(連坐標系都要講出身的)。
因為地球是一個赤道略寬兩極略扁的不規(guī)則的梨形球體纺涤,其表面是一個不可展平的曲面译暂,所以運用任何數(shù)學方法進行這種投影轉(zhuǎn)換都會產(chǎn)生誤差和變形,為按照不同的需求縮小誤差撩炊,就產(chǎn)生了各種投影方法外永。比如這篇很有意思的文章你所不知的有趣投影方法。也可以看文末的幾幅地球不同投影的圖拧咳。
EPSG Code
如果經(jīng)常跟地圖打交道伯顶,你一定會被各種花式投影變換弄的眼花繚亂。EPSG:European Petroleum Survey Group (EPSG)成立于1986年骆膝,并在2005年重組為OGP(Internation Association of Oil & Gas Producers)祭衩,它負責維護并發(fā)布坐標參照系統(tǒng)的數(shù)據(jù)集參數(shù),以及坐標轉(zhuǎn)換描述谭网,該數(shù)據(jù)集被廣泛接受并使用汪厨,通過一個Web發(fā)布平臺進行分發(fā),同時提供了微軟Acess數(shù)據(jù)庫的存儲文件愉择,通過SQL 腳本文件劫乱,MySQL, Oracle 和PostgreSQL等數(shù)據(jù)庫也可使用织中。目前已有的橢球體,投影坐標系等不同組合都對應著不同的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的網(wǎng)站 epsg.io
pyproj小試牛刀
pyproj是PROJ4的python接口封裝谦趣,直接看一個官網(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
基于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()
# 根據(jù)當前的蘭伯特投影繪制
shp_df.plot(column="GDP_1994",colormap='Set1')
# 轉(zhuǎn)換到經(jīng)緯度坐標
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摘悴,也可以自己定義投影參數(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')
代碼和矢量地圖數(shù)據(jù)下載
百度網(wǎng)盤下載地址: http://pan.baidu.com/s/1c1BExH2
密碼請在關注微信公眾號stdrei
后,輸入projection
直接獲取舰绘。
拓展:同一個世界蹂喻,不同的面孔
鏈接 http://www.viewsoftheworld.net/?p=752
在不同投影下的這個世界。捂寿。口四。