文章來源于Python大數(shù)據(jù)分析，作者費弗里

本文對應(yīng)代碼已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

在上一篇文章中我們對geopandas中的數(shù)據(jù)結(jié)構(gòu)展開了較為全面的學習，其中涉及到面積長度等計算的過程中提到了具體的計算結(jié)果與所選擇的投影坐標系關(guān)系密切，投影坐標系選擇的不恰當會帶來計算結(jié)果的偏差，直接關(guān)乎整個分析過程的有效與否。

作為基于geopandas的空間數(shù)據(jù)分析系列文章的第二篇，通過本文你將會學習到geopandas中的坐標參考系管理。

2 坐標參考系基礎(chǔ)

2.1 CRS

在一個二維的平面中，我們可以使用如圖1所示的坐標系統(tǒng)，通過坐標唯一確定點的位置：

圖1

現(xiàn)實世界中的地球作為一個球體，當我們想要用同樣的方式利用坐標來唯一確定地球球面上的某個位置時，需要一套適應(yīng)球體形狀的坐標系統(tǒng)。

而當我們想要在紙面或電腦屏幕上繪制平面地圖時，就又需要有一套將地球球面展平的方法。

上述的這些用于在不同情況下定義對象位置信息的坐標系統(tǒng)，就稱為坐標參考系統(tǒng)（Coordinate Reference System，下文統(tǒng)稱CRS）：

圖2

CRS可細分為地理坐標系和投影坐標系。

2.1.1 地理坐標系

以弧度制下度數(shù)為單位的地理坐標系（Geographic Coordinate Systems）幫助我們定位物體在地球球面上的具體位置以及繪制球體地圖：

圖3　WGS84地理坐標系示意圖

地理坐標系以地表上確定的某一個點為原點，創(chuàng)建了包裹全球的網(wǎng)格，譬如WGS84，將本初子午線與赤道的交點作為原點（圖4）：

圖4 WGS84地理坐標系及其經(jīng)緯網(wǎng)格

2.1.2 投影坐標系

地理坐標系雖然解決了我們在地球球面上定位的問題，但緯度和經(jīng)度位置沒有使用統(tǒng)一的測量單位。

因為經(jīng)度不變的情況下，緯度每變化1單位因為是對固定弧長的映射，所以真實距離是固定不變的，緯度變化1度的真實距離恒等于：

地球極半徑千米

可是經(jīng)度每變化1單位對應(yīng)的真實距離要隨著緯度的變化而變化，經(jīng)度變化1度的真實距離為：

這就導(dǎo)致我們既不能直接在地理坐標系下精確度量幾何對象的長度、面積，也無法直接用地理坐標系在平面上繪制出幾何對象真實的形狀。

為了解決上述問題，各種各樣的投影坐標系（Projected Coordinate Systems）被開發(fā)出來（圖5，其中右下角為地理坐標系，其余均為投影坐標系）：

圖5 各種CRS

投影坐標系指的是從將3D球面展平為2D平面的一套數(shù)學計算方法，利用它可以優(yōu)化形狀、比例/距離以及面積的失真情況。

但實際情況中沒有在整個地球表面都能“三全其美”的投影坐標系，有些投影坐標系優(yōu)化形狀上的失真，有些投影坐標系優(yōu)化距離上的失真，有些投影坐標系專門針對面積失真進行優(yōu)化，而有些投影坐標系可以對局部區(qū)域進行三個方面上的優(yōu)化。

圖6　投影坐標系變換過程示意

常用的投影坐標系如橫軸墨卡托（Universal Transverse Mercator，簡稱UTM），基于經(jīng)度將全球等分為編號0-60的區(qū)域，且每個區(qū)域又進一步細分為南半球區(qū)域或北半球區(qū)域，譬如圖7所示為美國本土跨過的區(qū)域：

圖7

劃分出的每個區(qū)域，其原點位于左下角頂點，距離區(qū)域中軸線500千米（圖8）：

圖8

針對這樣劃分出的獨立區(qū)域利用墨卡托投影法創(chuàng)建各自獨立的坐標網(wǎng)格，這個過程可以通俗地理解為用圓筒包裹地球球體，從球心發(fā)散出的光穿過球體上每個位置點投射到外部圓筒內(nèi)壁從而完成3D向2D的變換：

圖9

當然，這樣做的后果是越靠近極點的幾何對象被拉伸形變得越嚴重（圖10），這也就是為什么俄羅斯疆域看起來如此龐大的原因：

圖10　世界各國真實大小與墨卡托投影后差別

2.2 常用CRS格式

通過前文我們了解到什么是CRS，而在計算機系統(tǒng)中要使用CRS，需要將其文檔化，下面我們來了解CRS兩種常見的文檔存儲格式。

2.2.1 Proj4

Proj4字符串是一種識別空間或坐標參考系統(tǒng)的簡潔方法，通過其定義的語法規(guī)則，將想要定義的CRS全部參數(shù)信息保存到一條字符串中。

• Proj4語法

Proj4字符串包含了一種CRS全部元素信息，用+連接每個元素定義部分，如下面的例子記錄了橫軸墨卡托北11區(qū)CRS對應(yīng)的Proj4字符串：

+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

它記錄了如下信息：

proj=utm：聲明投影方法為墨卡托
zone=11：聲明對應(yīng)北11區(qū)（因為這里是橫軸墨卡托所以擁有獨立分區(qū)，但并不是所有CRS都有分區(qū)，且在Proj4中區(qū)號加S才為南半球分區(qū)如11S，否則默認為北半球分區(qū)）
datum=WGS84：聲明基準面為WGS84(基準面是橢球體用來逼近某地區(qū)用的，因此各個國家都有各自的基準面。國內(nèi)常用的基準面有：BEIJING1954，XIAN1980，WGS84等)
units=m：聲明坐標系單位設(shè)置為米
ellps=WGS84：聲明橢球面（如何計算地球的圓度）使用WGS84

上述例子記錄了投影坐標系的Proj4，下面我們再來看看地理坐標系對應(yīng)的Proj4，如下例：

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

它記錄了如下信息：

proj=longlat：聲明這是一個地理坐標系
datum=WGS84：聲明基準面為WGS84
ellps=WGS84：聲明橢球面使用WGS84

與投影坐標系相比，沒有單位units的信息，因為地理坐標系通常單位為十進制度數(shù)；而上述兩個示例中都帶有towgs84=0,0,0，這是一個轉(zhuǎn)換因子，在需要進行數(shù)據(jù)轉(zhuǎn)換時使用。

2.2.2 EPSG編碼

EPSG（European Petroleum Survey Group）編碼，使用4或5位數(shù)字編碼來唯一確定已存在的一種CRS，可以在http://spatialreference.org/ref/epsg/中查看和搜索所有已知的EPSG與CRS對應(yīng)關(guān)系（圖11）：

圖11

或在QGIS中查看：

圖12

譬如對于重慶，因為地跨東經(jīng)105°11~110°11，中軸線距離108E更近，常用如下投影：

圖13

對應(yīng)的EPSG編碼為2381。

3 geopandas中的坐標參考系管理

至此，我們已經(jīng)對CRS有了較為全面的了解，打好了基礎(chǔ)，接下來我們來正式學習geopandas中的坐標參考系管理。

geopandas調(diào)用pyproj作為CRS管理的后端，因此所有可以被pyproj.CRS.from_user_input()接受的合法輸入同樣可以被geopandas識別，譬如針對上文所說的應(yīng)用于重慶區(qū)域繪圖的Xian 1980 / 3-degree Gauss-Kruger CM 108E：

• Proj4

import pyprojpyproj.CRS.from_user_input('+proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 +y_0=0 +ellps=IAU76 +units=m +no_defs') 圖14

• EPSG

pyproj.CRS.from_user_input(2381) 圖15

直接傳入字符串格式的EPSG亦可：

圖16

查看對應(yīng)的Proj4信息，關(guān)鍵參數(shù)與前面Proj4一致，只是以Proj4形式傳入時系統(tǒng)會視作創(chuàng)建未知CRS一樣，因此相對于官方CRS缺少了一些無關(guān)緊要的其他信息：

圖17

3.1 CRS的設(shè)置與再投影

在上一篇文章（數(shù)據(jù)科學學習手札74）基于geopandas的空間數(shù)據(jù)分析——數(shù)據(jù)結(jié)構(gòu)篇中我們介紹了創(chuàng)建GeoSeries和GeoDataFrame的方法。

實際上，現(xiàn)實的空間分析計算任務(wù)中，必須要為數(shù)據(jù)設(shè)置合適的CRS，在geopandas.GeoSeries()和geopandas.GeoDataFrame()中就包含參數(shù)crs。

下面我們舉例說明，還是先用到geopandas自帶的世界國家地區(qū)數(shù)據(jù)，我們從中選擇中國（堅持一個中國，我們將臺灣地區(qū)組合進國土中）：

import geopandas as gpdworld = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) # 利用name字段選擇中國區(qū)域 china = world.loc[world['name'].isin(['China', 'Taiwan'])] china 圖18

查看其crs屬性即為其對應(yīng)CRS，為WGS84對應(yīng)的EPSG:4326，在當前的CRS下將其繪制出來：

圖19

利用to_crs()將其再投影到EPSG:2381并進行繪制：

圖20

通過比較可以發(fā)現(xiàn)，再投影之后的中國形變失真情況得到緩解，且坐標系單位范圍也發(fā)生了變化（EPSG:2381單位：米），接下來我們參考谷歌地圖上點擊出的重慶渝中區(qū)某地坐標：

圖21

基于此創(chuàng)建只包含一個點的GeoSeries，嘗試將其與EPSG:2381下的中國地圖一同繪制：

from shapely import geometry import matplotlib.pyplot as pltcq = gpd.GeoSeries([geometry.Point([106.561203, 29.558078])],crs='EPSG:4326')fig, ax = plt.subplots() china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8) cq.plot(ax=ax, color='orange', markersize=100, marker='x') plt.xticks(rotation=20) 圖22

可以看出我們創(chuàng)建在重慶境內(nèi)的點并沒有繪制在正確的位置，接下來我們對cq進行再投影，再嘗試將其與EPSG:2381下的中國繪制在一起：

fig, ax = plt.subplots() china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8) # 先再投影到EPSG:2381 cq.to_crs(crs='EPSG:2381').plot(ax=ax, color='orange', markersize=100, marker='x') plt.xticks(rotation=20) 圖23

這時我們定義的點被繪制到正確的位置。

同樣地，可以在投影后計算更為準確的面積，這里舉一個粗糙的例子（實際計算國土面積不會這樣粗糙），以中國中軸線東經(jīng)104.19度最靠近的105度經(jīng)線對應(yīng)的EPSG:2380為CRS計算面積：

圖24

如果直接用原來的ESPG:4326計算面積結(jié)果如下：

圖25

可以看出使用ESPG:2380計算出的面積比較接近大家記憶中的960萬平方公里。

以上就是本文的全部內(nèi)容，如有筆誤之處望斧正！

下一篇文章將會介紹geopandas中的文件IO與基礎(chǔ)地圖制作，敬請期待。

-END-

往期精彩回顧適合初學者入門人工智能的路線及資料下載機器學習在線手冊深度學習在線手冊AI基礎(chǔ)下載（pdf更新到25集）備注：加入本站微信群或者qq群，請回復(fù)“加群”獲取一折本站知識星球優(yōu)惠券，請回復(fù)“知識星球”喜歡文章，點個在看

與50位技術(shù)專家面對面20年技術(shù)見證，附贈技術(shù)全景圖

總結(jié)

以上是生活随笔為你收集整理的Python地信专题 | 基于geopandas的空间数据分析-坐标参考系篇的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。