本文對應程式碼已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

　　在上一篇文章中我們對geopandas中的資料結構展開了較為全面的學習，其中涉及到面積長度等計算的過程中提到了具體的計算結果與所選擇的投影座標系關係密切，投影座標系選擇的不恰當會帶來計算結果的偏差，直接關乎整個分析過程的有效與否。
　　作為基於geopandas的空間資料分析系列文章的第二篇，通過本文你將會學習到geopandas中的座標參考系管理。

2 座標參考系基礎

2.1 CRS

　　在一個二維的平面中，我們可以使用如圖1所示的座標系統，通過座標\((x_{0},y_{0})\)唯一確定點的位置：

圖1
　　現實世界中的地球作為一個球體，當我們想要用同樣的方式利用座標\((\phi_{0},\lambda_{0})\)來唯一確定地球球面上的某個位置時，需要一套適應球體形狀的座標系統。而當我們想要在紙面或電腦螢幕上繪製平面地圖時，就又需要有一套將地球球面展平的方法，上述的這些用於在不同情況下定義物件位置資訊的座標系統，就稱為座標參考系統（Coordinate Reference System，下文統稱CRS）：
圖2
　　CRS可細分為地理座標系和投影座標系。

2.1.1 地理座標系

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

圖3　WGS84地理座標系示意圖
　　地理座標系以地表上確定的某一個點為原點\((0,0)\)，建立了包裹全球的網格，譬如WGS84，將本初子午線與赤道的交點作為原點（圖5）：
圖4 WGS84地理座標系及其經緯網格

2.1.2 投影座標系

　　地理座標系雖然解決了我們在地球球面上定位的問題，但緯度和經度位置沒有使用統一的測量單位，因為經度不變的情況下，緯度每變化1單位因為是對固定弧長的對映，所以真實距離是固定不變的，緯度變化1度的真實距離恆等於：
\[ \begin{split} 2\times\pi\times地球極半徑/360&\approx2\times3.1415926\times6356.9088/360\\&=110.95(千米) \end{split} \]

可是經度每變化1單位對應的真實距離要隨著緯度的變化而變化，經度變化1度的真實距離為：
\[ \begin{split} (2\times\pi\times地球赤道半徑/360)\times\cos(當地緯度)&\approx(2\times3.1415926\times6377.830/360)\times\cos(當地緯度)\\&=111.314\times\cos(當地緯度)(千米) \end{split} \]
這就導致我們既不能直接在地理座標系下精確度量幾何物件的長度、面積，也無法直接用地理座標系在平面上繪製出幾何物件真實的形狀。為了解決上述問題，各種各樣的投影座標系（Projected Coordinate Systems）被開發出來（圖4，其中右下角為地理座標系，其餘均為投影座標系）：
圖5 各種CRS
　　投影座標系指的是從將3D球面展平為2D平面的一套數學計算方法，利用它可以優化形狀、比例/距離以及面積的失真情況，但實際情況中沒有在整個地球表面都能“三全其美”的投影座標系，有些投影座標系優化形狀上的失真，有些投影座標系優化距離上的失真，有些投影座標系專門針對面積失真進行優化，而有些投影座標系可以對區域性區域進行三個方面上的優化。
圖6　投影座標系變換過程示意
　　常用的投影座標系如橫軸墨卡託（Universal Transverse Mercator，簡稱UTM），基於經度將全球等分為編號0-60的區域，區分南北半球，譬如圖7所示為美國本土跨過的區域：
圖7
　　劃分出的每個區域，其原點\((0,0)\)位於左下角頂點，距離區域中軸線500千米（圖8）：
圖8
　　針對這樣劃分出的獨立區域利用墨卡託投影法建立各自獨立的座標網格，這個過程可以通俗地理解為用圓筒包裹地球球體，從球心發散出的光穿過球體上每個位置點投射到外部圓筒內壁從而完成3D向2D的變換：
圖8
　　當然，這樣做的後果是越靠近極點的幾何物件被拉伸形變得越嚴重（圖9），這也就是為什麼俄羅斯疆域看起來如此龐大的原因：
圖9　世界各國真實大小與墨卡託投影后差別

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對應關係（圖10）：

圖10
　　或在QGIS中檢視：
圖11
　　譬如對於重慶，因為地跨東經105°11~110°11，中軸線距離108E更近，常用如下投影：
圖12
　　對應的EPSG編碼為2381。

3 geopandas中的座標參考系管理

　　至此，我們已經對CRS有了較為全面的瞭解，打好了基礎，接下來我們來正式學習geopandas中的座標參考系管理，geopandas呼叫pyproj作為CRS管理的後端，因此所有可以被pyproj.CRS.from_user_input()接受的合法輸入同樣可以被geopandas識別，譬如針對上文所說的應用於重慶區域繪圖的Xian 1980 / 3-degree Gauss-Kruger CM 108E：

• Proj4

import pyproj

pyproj.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')

圖13

• EPSG

pyproj.CRS.from_user_input(2381)

圖14

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

圖15
　　檢視對應的Proj4資訊，關鍵引數與前面Proj4一致，只是以Proj4形式傳入時系統會視作建立未知CRS一樣，因此相對於官方CRS缺少了一些無關緊要的其他資訊：
圖16

3.1 CRS的設定與再投影

　　在上一篇文章（資料科學學習手札74）基於geopandas的空間資料分析——資料結構篇中我們介紹了建立GeoSeries和GeoDataFrame的方法，實際上，現實的空間分析計算任務中，必須要為資料設定合適的CRS，在geopandas.GeoSeries()和geopandas.GeoDataFrame()中就包含引數crs，下面我們舉例說明，還是先用到geopandas自帶的世界國家地區資料，我們從中選擇中國（堅持一箇中國，我們將臺灣地區組合進國土中）：

import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# 利用name欄位選擇中國區域
china = world.loc[world['name'].isin(['China', 'Taiwan'])]
china

圖17

　　檢視其crs屬性即為其對應CRS，為WGS84對應的EPSG:4326，在當前的CRS下將其繪製出來：

圖18
　　利用to_crs()將其再投影到EPSG:2381並進行繪製：
圖19
　　通過比較可以發現，再投影之後的中國形變失真情況得到緩解，且座標系單位範圍也發生了變化（EPSG:2381單位：米），接下來我們參考谷歌地圖上點擊出的重慶渝中區某地座標：
圖20
　　基於此建立只包含一個點的GeoSeries，嘗試將其與EPSG:2381下的中國地圖一同繪製：

from shapely import geometry
import matplotlib.pyplot as plt

cq = 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)

圖21

　　可以看出我們建立在重慶境內的點並沒有繪製在正確的位置，接下來我們對cq進行再投影，再嘗試將其與EPSG:2381下的中國繪製在一起：

圖22
　　這時我們定義的點被繪製到正確的位置。
　　同樣地，可以在投影后計算更為準確的面積，這裡舉一個粗糙的例子（實際計算國土面積不會這樣粗糙），以中國中軸線東經104.19度最靠近的105度經線對應的EPSG:2380為CRS計算面積：
圖23
　　如果直接用原來的ESPG:4326計算面積結果如下：
圖24
　　可以看出使用ESPG:2380計算出的面積比較接近大家記憶中的960萬平方公里。

　　以上就是本文的全部內容，如有筆誤之處望斧正！
　　下一篇文章將會介紹geopandas中的檔案IO與基礎地圖製作，敬請期待