Python地信专题 | 基于geopandas的空间数据分析-坐标参考系篇
文章來源于Python大數據分析，作者費弗里
本文對應代碼已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes
1 簡介
在上一篇文章中我們對geopandas中的數據結構展開了較為全面的學習，其中涉及到面積長度等計算的過程中提到了具體的計算結果與所選擇的投影坐標系關系密切，投影坐標系選擇的不恰當會帶來計算結果的偏差，直接關乎整個分析過程的有效與否。
作為基于geopandas的空間數據分析系列文章的第二篇，通過本文你將會學習到geopandas中的坐標參考系管理。
2 坐標參考系基礎
2.1 CRS
在一個二維的平面中，我們可以使用如圖1所示的坐標系統，通過坐標唯一確定點的位置：
圖1現實世界中的地球作為一個球體，當我們想要用同樣的方式利用坐標來唯一確定地球球面上的某個位置時，需要一套適應球體形狀的坐標系統。
而當我們想要在紙面或電腦屏幕上繪制平面地圖時，就又需要有一套將地球球面展平的方法。
上述的這些用于在不同情況下定義對象位置信息的坐標系統，就稱為坐標參考系統（Coordinate Reference System，下文統稱CRS）：
圖2CRS可細分為地理坐標系和投影坐標系。
2.1.1 地理坐標系
以弧度制下度數為單位的地理坐標系（Geographic Coordinate Systems）幫助我們定位物體在地球球面上的具體位置以及繪制球體地圖：
圖3 WGS84地理坐標系示意圖地理坐標系以地表上確定的某一個點為原點，創建了包裹全球的網格，譬如WGS84，將本初子午線與赤道的交點作為原點（圖4）：
圖4 WGS84地理坐標系及其經緯網格2.1.2 投影坐標系
地理坐標系雖然解決了我們在地球球面上定位的問題，但緯度和經度位置沒有使用統一的測量單位。
因為經度不變的情況下，緯度每變化1單位因為是對固定弧長的映射，所以真實距離是固定不變的，緯度變化1度的真實距離恒等于：
地球極半徑千米
可是經度每變化1單位對應的真實距離要隨著緯度的變化而變化，經度變化1度的真實距離為：
這就導致我們既不能直接在地理坐標系下精確度量幾何對象的長度、面積，也無法直接用地理坐標系在平面上繪制出幾何對象真實的形狀。
為了解決上述問題，各種各樣的投影坐標系（Projected Coordinate Systems）被開發出來（圖5，其中右下角為地理坐標系，其余均為投影坐標系）：
圖5 各種CRS投影坐標系指的是從將3D球面展平為2D平面的一套數學計算方法，利用它可以優化形狀、比例/距離以及面積的失真情況。
但實際情況中沒有在整個地球表面都能“三全其美”的投影坐標系，有些投影坐標系優化形狀上的失真，有些投影坐標系優化距離上的失真，有些投影坐標系專門針對面積失真進行優化，而有些投影坐標系可以對局部區域進行三個方面上的優化。
圖6 投影坐標系變換過程示意常用的投影坐標系如橫軸墨卡托（Universal Transverse Mercator，簡稱UTM），基于經度將全球等分為編號0-60的區域，且每個區域又進一步細分為南半球區域或北半球區域，譬如圖7所示為美國本土跨過的區域：
圖7劃分出的每個區域，其原點位于左下角頂點，距離區域中軸線500千米（圖8）：
圖8針對這樣劃分出的獨立區域利用墨卡托投影法創建各自獨立的坐標網格，這個過程可以通俗地理解為用圓筒包裹地球球體，從球心發散出的光穿過球體上每個位置點投射到外部圓筒內壁從而完成3D向2D的變換：
圖9當然，這樣做的后果是越靠近極點的幾何對象被拉伸形變得越嚴重（圖10），這也就是為什么俄羅斯疆域看起來如此龐大的原因：
圖10 世界各國真實大小與墨卡托投影后差別2.2 常用CRS格式
通過前文我們了解到什么是CRS，而在計算機系統中要使用CRS，需要將其文檔化，下面我們來了解CRS兩種常見的文檔存儲格式。
2.2.1 Proj4
Proj4字符串是一種識別空間或坐標參考系統的簡潔方法，通過其定義的語法規則，將想要定義的CRS全部參數信息保存到一條字符串中。
Proj4語法
Proj4字符串包含了一種CRS全部元素信息，用+連接每個元素定義部分，如下面的例子記錄了橫軸墨卡托北11區CRS對應的Proj4字符串：
+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0它記錄了如下信息：
proj=utm：聲明投影方法為墨卡托
zone=11：聲明對應北11區（因為這里是橫軸墨卡托所以擁有獨立分區，但并不是所有CRS都有分區，且在Proj4中區號加S才為南半球分區如11S，否則默認為北半球分區）
datum=WGS84：聲明基準面為WGS84(基準面是橢球體用來逼近某地區用的，因此各個國家都有各自的基準面。國內常用的基準面有：BEIJING1954，XIAN1980，WGS84等)
units=m：聲明坐標系單位設置為米
ellps=WGS84：聲明橢球面（如何計算地球的圓度）使用WGS84
上述例子記錄了投影坐標系的Proj4，下面我們再來看看地理坐標系對應的Proj4，如下例：
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0它記錄了如下信息：
proj=longlat：聲明這是一個地理坐標系
datum=WGS84：聲明基準面為WGS84
ellps=WGS84：聲明橢球面使用WGS84
與投影坐標系相比，沒有單位units的信息，因為地理坐標系通常單位為十進制度數；而上述兩個示例中都帶有towgs84=0,0,0，這是一個轉換因子，在需要進行數據轉換時使用。
2.2.2 EPSG編碼
EPSG（European Petroleum Survey Group）編碼，使用4或5位數字編碼來唯一確定已存在的一種CRS，可以在http://spatialreference.org/ref/epsg/中查看和搜索所有已知的EPSG與CRS對應關系（圖11）：
圖11或在QGIS中查看：
圖12譬如對于重慶，因為地跨東經105°11~110°11，中軸線距離108E更近，常用如下投影：
圖13對應的EPSG編碼為2381。
3 geopandas中的坐標參考系管理
至此，我們已經對CRS有了較為全面的了解，打好了基礎，接下來我們來正式學習geopandas中的坐標參考系管理。
geopandas調用pyproj作為CRS管理的后端，因此所有可以被pyproj.CRS.from_user_input()接受的合法輸入同樣可以被geopandas識別，譬如針對上文所說的應用于重慶區域繪圖的Xian 1980 / 3-degree Gauss-Kruger CM 108E：
Proj4
EPSG
直接傳入字符串格式的EPSG亦可：
圖16查看對應的Proj4信息，關鍵參數與前面Proj4一致，只是以Proj4形式傳入時系統會視作創建未知CRS一樣，因此相對于官方CRS缺少了一些無關緊要的其他信息：
圖173.1 CRS的設置與再投影
在上一篇文章（數據科學學習手札74）基于geopandas的空間數據分析——數據結構篇中我們介紹了創建GeoSeries和GeoDataFrame的方法。
實際上，現實的空間分析計算任務中，必須要為數據設置合適的CRS，在geopandas.GeoSeries()和geopandas.GeoDataFrame()中就包含參數crs。
下面我們舉例說明，還是先用到geopandas自帶的世界國家地區數據，我們從中選擇中國（堅持一個中國，我們將臺灣地區組合進國土中）：
import geopandas as gpdworld = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) # 利用name字段選擇中國區域 china = world.loc[world['name'].isin(['China', 'Taiwan'])] china 圖18查看其crs屬性即為其對應CRS，為WGS84對應的EPSG:4326，在當前的CRS下將其繪制出來：
圖19利用to_crs()將其再投影到EPSG:2381并進行繪制：
圖20通過比較可以發現，再投影之后的中國形變失真情況得到緩解，且坐標系單位范圍也發生了變化（EPSG:2381單位：米），接下來我們參考谷歌地圖上點擊出的重慶渝中區某地坐標：
圖21基于此創建只包含一個點的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可以看出我們創建在重慶境內的點并沒有繪制在正確的位置，接下來我們對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這時我們定義的點被繪制到正確的位置。
同樣地，可以在投影后計算更為準確的面積，這里舉一個粗糙的例子（實際計算國土面積不會這樣粗糙），以中國中軸線東經104.19度最靠近的105度經線對應的EPSG:2380為CRS計算面積：
圖24如果直接用原來的ESPG:4326計算面積結果如下：
圖25可以看出使用ESPG:2380計算出的面積比較接近大家記憶中的960萬平方公里。
以上就是本文的全部內容，如有筆誤之處望斧正！
下一篇文章將會介紹geopandas中的文件IO與基礎地圖制作，敬請期待。
-END-
往期精彩回顧適合初學者入門人工智能的路線及資料下載機器學習在線手冊深度學習在線手冊AI基礎下載（pdf更新到25集）備注：加入本站微信群或者qq群，請回復“加群”獲取一折本站知識星球優惠券，請回復“知識星球”喜歡文章，點個在看 與50位技術專家面對面20年技術見證，附贈技術全景圖
總結
以上是生活随笔為你收集整理的Python地信专题 | 基于geopandas的空间数据分析-坐标参考系篇的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 抗击肺炎：新冠肺炎疫情数据可视化及疫情预
- 下一篇： Python地信专题 | 基于geopa