本文主要介紹GeoPandas的基本使用方法，以繪製簡單的地圖。GeoPandas是一個Python開源專案，旨在提供豐富而簡單的地理空間資料處理介面。GeoPandas擴充套件了Pandas的資料型別，並使用matplotlib進行繪圖。GeoPandas官方倉庫地址為：GeoPandas。GeoPandas的官方檔案地址為：GeoPandas-doc。

本文所有程式碼見：Python-Study-Notes

GeoPandas推薦使用Python3.7版本及以上，執行環境最好是linux系統。GeoPandas安裝命令如下：

pip install geopandas

如果上述命令安裝出問題，則推薦使用conda安裝GeoPandas，命令如下：

conda install geopandas

或：

conda install --channel conda-forge geopandas

# jupyter notebook環境去除warning
import warnings
warnings.filterwarnings("ignore")

# 檢視geopandas版本
import geopandas as gpd

gpd.__version__

'0.10.2'

1 基礎結構

GeoPandas在Pandas的基礎上增加對地理空間資料的支援。其主要結構為GeoSeries和GeoDataFrame。如下圖所示，整張圖片表示一個GeoDataFrame，行表示地圖上的某一區域，列表示該區域的屬性資料，其中geometry表示該區域的位置資訊。每一個區域的位置形狀用一個或多個幾何物件表示，幾何物件在GeoPandas中以GeoSeries格式資料構成。

1.1 GeoSeries

1.1.1 GeoSeries基礎元素

在Geogeopandas中，GeoSeries用於構成基礎的幾何物件。GeoSeries由Python庫Shapely實現的幾何物件構成，這些幾何元素可以是點（集）、線（集）、面（集）：

• Points / Multi-Points

• Lines / Multi-Lines

• Polygons / Multi-Polygons

在jupyter notebook中， 可以直接以svg格式展示GeoSeries中的單個元素，使用如下：

Points 點

# Point(x,y)用於建立單個點
from shapely.geometry import Point
p = gpd.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)],index=['p1','p2','p3'])
p

p1    POINT (1.00000 1.00000)
p2    POINT (2.00000 2.00000)
p3    POINT (3.00000 3.00000)
dtype: geometry

# 展示單個元素
p[0]

Multi-Points 點集

# MultiPoint([(x1,y1), ... ,(xn,yn)])用於建立多個點構成的集合
from shapely.geometry import MultiPoint
mp = gpd.GeoSeries([MultiPoint([(0, 1), (1, 0)]),
                    MultiPoint([(0, 0), (1, 1)])],
                    index=['mp1', 'mp2'])
mp

mp1    MULTIPOINT (0.00000 1.00000, 1.00000 0.00000)
mp2    MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)
dtype: geometry

# 展示單個元素
mp[0]

Lines 線

# LineString([(x1,y1), ... , (xn,yn)])用於建立單條線段
from shapely.geometry import LineString
l = gpd.GeoSeries([LineString([(0.5, 0), (1, 2), (0, 1)]),
                   LineString([(0, 0), (3.3, 1)])],
                   index=['l1', 'l2'])
l

l1    LINESTRING (0.50000 0.00000, 1.00000 2.00000, ...
l2        LINESTRING (0.00000 0.00000, 3.30000 1.00000)
dtype: geometry

# 展示單個元素
l[0]

MultiLineString 線集

# MultiLineString([LineString1,LineString2])用於建立多條線段的集合
from shapely.geometry import MultiLineString
ml = gpd.GeoSeries([MultiLineString([[(1, 0), (2, 1)], [(-1, 0), (0, 5), (1, 0)]]),
               MultiLineString([[(0.5, 0), (1, 2), (0, 1)]])],
               index=['m1','m2'])
ml

m1    MULTILINESTRING ((1.00000 0.00000, 2.00000 1.0...
m2    MULTILINESTRING ((0.50000 0.00000, 1.00000 2.0...
dtype: geometry

# 展示單個元素
ml[0]

Polygon 面

在Polygon元素，其建立函數為：

# shell表示單個多邊形座標，hole表示是否在內部建立空洞
Polygon(shell=None, holes=None)

# 建立單個多邊形
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                     Polygon([(0, 0), (1, 1), (1, 2)])],
                    index=['pol1','pol2'])
pol

pol1    POLYGON ((0.00000 0.00000, 4.00000 1.00000, 5....
pol2    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 1....
dtype: geometry

# 展示單個元素
pol[0]

# 建立帶空洞的單個多邊形
pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)],
                     [((1, 1),(1,2),(2,2),(2,1)),
                     ((4,2),(4,3),(1,5))])])
pol

0    POLYGON ((0.00000 0.00000, 4.00000 1.00000, 5....
dtype: geometry

# 展示單個元素
pol[0]

MultiPolygon 面集

# MultiPolygon([Polygon,Polygon])用於建立多個多邊形的集合
from shapely.geometry import Polygon, MultiPolygon

mpol = gpd.GeoSeries([MultiPolygon([Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]),
                                    Polygon([(2, 3), (3, 2), (3, 3)])])])
mpol

0    MULTIPOLYGON (((0.00000 0.00000, 0.00000 1.000...
dtype: geometry

# 展示單個元素
mpol[0]

1.1.2 GeoSeries屬性與方法

GeoSeries常用屬性如下：

• area: 計算GeoSeries每一個元素的面積，面積單位取決於座標系
• bounds: 獲得GeoSeries每一個元素的最小外接矩形的左上角座標和右下角座標
• total_bounds：計算整個GeoSeries物件的最小外接矩形的左上角座標和右下角座標
• length：計算GeoSeries每一個元素的周長，周長單位取決於座標系
• geom_type: 獲得GeoSeries每一個元素的型別
• is_valid: 判斷GeoSeries每一個元素是否構成合理的多邊形
• centroid：獲得GeoSeries每一個元素的幾何中心

# 建立包含兩個多邊形的GeoSeries
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                    Polygon([(0, 0), (1, 1), (1, 8)])],
                    index=['pol1','pol2'])

# 計算各個元素的面積
pol.area

pol1    21.0
pol2     3.5
dtype: float64

# 計算各元素的最小外接矩形的頂點資料
pol.bounds

# 計算整個物件的最小外接矩形的頂點資料
pol.total_bounds

array([0., 0., 5., 8.])

# 計算GeoSeries每一個元素的周長
pol.length

pol1    19.762298
pol2    16.476471
dtype: float64

# 獲得GeoSeries每一個元素的型別
pol.type

pol1    Polygon
pol2    Polygon
dtype: object

# 判斷GeoSeries每一個元素是否構成合理的多邊形
pol.is_valid

pol1    True
pol2    True
dtype: bool

# 獲得GeoSeries每一個元素的幾何中心
pol.centroid

pol1    POINT (1.88889 3.00000)
pol2    POINT (0.66667 3.00000)
dtype: geometry

GeoSeries常用方法如下：

• distance(): 計算GeoSeries每一個元素到某一幾何體的距離，具體使用檢視官方檔案
• representative_point()：返回位於GeoSeries每一個元素中的座標點，該點不會是幾何中心。
• to_crs()：改變GeoSeries的座標參考系
• plot()：繪製整個GeoSeries

# 建立包含兩個多邊形的GeoSeries
from shapely.geometry import Polygon

pol = gpd.GeoSeries([Polygon([(0, 0), (4, 1), (5, 3), (0, 7)]),
                    Polygon([(0, 0), (1, 1), (1, 8)])],
                    index=['pol1','pol2'])

# 計算距離
point = Point(1, 0)
pol.distance(point)

pol1    0.242536
pol2    0.707107
dtype: float64

# 返回位於GeoSeries每一個元素中的座標點
pol.representative_point()

pol1    POINT (1.25000 5.00000)
pol2    POINT (0.78125 4.50000)
dtype: geometry

# 繪製整個GeoSeries
pol.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00a24090>

1.2 GeoDataFrame

GeoDataFrame是包含GeoSeries的表格資料結構，GeoDataFrame最重要的屬性是它始終具有一個表示地理狀態資訊的GeoSeries列，即geometry列。呼叫GeoDataFrame的方法或屬性時，將始終作用於geometry列。

下面用一個世界地圖的範例介紹GeoDataFrame常用函數的使用方法，具體函數的介面請查閱官方檔案。

⚠⚠⚠注意該世界地圖資料集來自geopandas自帶資料集，部分資料資料可能缺失或錯誤，且資料集僅可用於學習不做其他用途

import geopandas as gpd
# 載入自帶世界地圖資料集
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# 檢視資料集
# 人口，大洲，國名，國家縮寫，ISO國家程式碼，gdp，地理位置資料
world.head()

# 繪製GeoDataFrame
# column選擇為地圖著色的資料列
# kind設定圖片型別，某些型別會導致column引數不可用
# cmap設定matplotlib調色盤
# legend設定是否顯示圖例
# figsize設定圖片大小
# legend_kwds設定matplotlib的legend引數，該引數需要根據實際圖例形式設定
world.plot(column="continent",kind="geo", 
            cmap="coolwarm",legend=True,figsize=(9,6), 
            legend_kwds={'loc': (0,1.02),'ncol':4})

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00551790>

# 獲得設定地理位置資訊的列名
world.geometry.name
# 修改地理位置資訊的列名，一般不要修改
# world = world.rename(columns={'geometry': 'borders'}).set_geometry('borders')
# world.geometry.name

'geometry'

# 計算各個區域的幾何中心，並新增列centroid_column
world['centroid_column'] = world.centroid
# 以centroid_column表示地理位置資訊
world = world.set_geometry('centroid_column')
world.head()

# 根據gdp繪圖
world.plot(column="gdp_md_est",legend=True)

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe00489190>

# 疊加幾何中心繪圖
ax = world["geometry"].plot()
world["centroid_column"].plot(ax=ax, color="red")

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe002cfe90>

# 設定資料小數精度
gpd.options.display_precision = 1
world.head()

此外，pandas的資料操作也可以用於geopandas。

# 檢視資料
world.iloc[:4,[3,4]]

# 檢視人口小於10萬的國家地區
world.loc[world['pop_est']<1e5,['pop_est','name','gdp_md_est','geometry']]

2 資料集基本操作

2.1 資料集讀寫

對於包含地理資訊資料的檔案（例如GeoPackage、GeoJSON、Shapefile），geopandas提供了read_file()函數介面來自動讀取對應的型別檔案並返回GeoDataFrame物件。此外，geopandas也提供to_file()函數介面以實現不同型別的地理資訊檔案寫入。下面以一個實際例子來介紹資料集的讀寫。

對於地理資訊資料，通常的獲取方式是從專業網站下載，或從谷歌地圖和高德地圖獲取。計算機領域常用的地理資訊資料檔案格式為GeoJSON，關於GeoJSON格式的中國地理資訊資料，推薦從阿里雲資料視覺化平臺獲取：DataV.GeoAtlas。DataV.GeoAtlas提供了我國各省各市區縣地圖資訊（不包含鄉鎮）。

如下圖所示，開啟網頁後，點選或查詢所需的區域->選擇需要的資料格式->選擇是否包含子區域->複製JSON API->read_file()函數讀取檔案即可。例如本文選擇江蘇省的地理資訊，包含子區域的意思為是否包含江蘇省內地區的地理資料。

import geopandas as gpd

# 讀取資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")
# 展示資料
# abcode：行政地區編碼；
# childrenNum，子區域個數，在這裡表示城市所轄區縣數，
# parent：父區域編碼，320000表示江蘇省的意思
data.head()

# 展示江蘇省地圖，橫軸為經度，縱軸為緯度
data.plot(figsize=(6,6),color='green')

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe002057d0>

# 寫入Shapefile
data.to_file("countries.shp")

# 寫入GeoJSON
# 預設格式是Shapefile,通過driver引數可以設定其他型別
data.to_file("countries.geojson", driver='GeoJSON')

2.2 資料集篩選

GeoPandas提供了各種資料集篩選函數，以提取感興趣區域。其在0.1.0版本中新增了Bounding Box Filter，在0.7.0版本中新增了Geometry Filter和Row Filter，以及其他需要Fiona庫的過濾器。

2.2.1 Bounding Box Filter

GeoPandas在0.1.0版本引入了Bounding Box Filter。Bounding Box Filter用於提取與邊界框相交的子區域，注意邊界框格式為(左下角x, 左下角y, 右上角x, 右上角y)。

# 框選資料
import geopandas as gpd

# 設定框選範圍
bbox = (119, 32, 121, 33)
# 讀取資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json", bbox=bbox)

# 結果是篩選與bbox相交的子區域
data.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbe000e1c90>

# 檢視邊界框與哪些子區域相交
# 在地圖上繪製邊界框
from shapely import geometry

# 設定框選範圍
bbox = (119, 32, 121, 33)
# 讀取資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")
ax = data.plot()
# 在地圖上繪製bbox框
ax = gpd.GeoSeries([geometry.box(minx=bbox[0],miny=bbox[1], 
                                 maxx=bbox[2],maxy=bbox[3]).boundary]).plot(ax=ax, color='red')

2.2.2 Geometry Filter

GeoPandas在0.7.0版本引入了Geometry Filter。與Bounding Box Filter類似，Geometry Filter也是篩選與指定區域相交的子區域，但Geometry Filter輸入的是幾何圖形。

# 檢視幾何圖形與哪些子區域相交
from shapely import geometry

# 設定框選範圍
p = geometry.Polygon([(119,32),(120,33), (121, 33)])
# 讀取資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",mask=p)
# 篩選區域
data.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6fbf590>

# 檢視幾何圖形與哪些子區域相交
from shapely import geometry

# 設定框選範圍
p = geometry.Polygon([(119,32),(120,33), (121, 33)])
# 讀取資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")

ax = data.plot()
ax = gpd.GeoSeries([p.boundary]).plot(ax=ax, color='red')

2.2.3 Row Filter

GeoPandas在0.7.0版本引入了Row Filter。Row Filter用於提取指定行範圍的資料。

# 讀取前3行資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",rows=3)
data

# 讀取第1行和第2行資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json",rows=slice(1, 3))
data

2.2.4 cx索引器

除了標準的Pandas方法外，GeoPandas還提供使用cx索引器進行基於座標的索引。通過cx索引器可以挑選特定區域。

# 讀取世界地圖
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# 篩選與經度範圍為80到200，緯度小於0相交的區域，這裡座標代表經緯度。
southern_world = world.cx[80:200, :0]

southern_world.plot(figsize=(9, 6))

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6f03150>

# 檢視與哪些子區域相交
# 讀取世界地圖
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

ax = world.plot(color='grey',alpha=0.5)

# 篩選與經度範圍為80到200，緯度小於0相交的區域，這裡座標代表經緯度。
southern_world = world.cx[80:200, :0]

southern_world.plot(ax=ax,figsize=(9, 6))

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6ee2450>

3 座標參考系

3.1 座標參考系簡介

地信、計算機和遙感等領域的從業人員或多或少都會接觸地理資訊系統（GIS，Geographic Information System）的相關知識。所謂GIS簡單來說就是一個以計算機為核心，對地理空間位置相關資料進行建立、管理、分析、繪製和展示的多功能整合資訊系統。繪製地圖，則需要了解GIS中的座標參考系。本文只是簡單介紹座標參考系的相關內容。相關內容總結於以下文章，想要具體瞭解更多內容也可以看看這些文章。

• Coordinate Reference System and Spatial Projection
• 聊聊GIS中的座標系
• 基於geopandas的空間資料分析——座標參考系篇
• 前後端都要懂一點的 GIS
• GIS基礎知識

如下圖所示，對於二維平面上的任何一點，我們可以通過如下的笛卡爾座標系確定其唯一座標，並用\((x,y)\)表示。

然而我們生活在一個恰好是「圓形」的三維地球上，一般來說要定義地球上物體的位置，需要一個適應地球形狀的三維座標系。但當我們在紙上或計算機螢幕上展示地球時，受限於儲存媒介（紙張、沙盤、書本、螢幕）原因，不得不使用二維平面來表達三維的地球。這種將三維球面轉換為二維平面所用到的座標系參考就是座標參考系統（CRS，Coordinate Reference System)，如下圖所示。

常見情況下，根據不同的轉換方式，CRS又可分為地理座標系和投影座標系。

地理座標系統

參考平面為球面，用球面座標來表示地球上的位置，叫做地理座標系統。這些座標系的座標單位通常為角度，如經緯度。所謂經度就是指離被稱為本初子午線的南北方向走線以東或以西的度數。緯度則是指地球某點與球心的連線和地球赤道面所成的線面角。如下圖所示：

在資料表示中，東經正數，西經為負數，經度範圍為-180度到180度。赤道以北為北緯，以南為南緯，緯度範圍為-90度到90度。按照常用地理座標系統表示的二維地圖如下所示：

投影座標系統

由於我們所在的地球實際不是規則球體，一經度和一緯度代表的實際距離是不相等的（赤道除外）。這導致我們無法通過經緯度座標精確計算某塊地區的面積，長寬。也無法精準表達某一塊地區的實際形狀。因此，有必要採取某些方法將地球表面投影到一個二維平面，這種地理投影方法也就是投影座標系，實際作用為把地球表面經緯度按照一定的數學規則轉換為平面座標。

投影座標系在二維平面中定義，其始終基於地理座標系，位置由二維座標標識，度量單位用米、千米等來表示。不同的投影方法會導致不同的投影結果，因此各種各樣的投影座標系統被提出。但是不管是哪種地圖投影方法，都會導致地圖投影變形。三維投影二維出現投影變形是不可避免的，但不會在任何位置都會產生變形。那些沒有發生變形的點或線稱為標準點線，距離標準點線越近，變形程度則越小。因此在實際投影某塊區域時，最好選擇標標準點線離該區域近的投影座標系。

關於常見投影方法的介紹可以閱讀本節開頭給出的參考文章。想要感受不同投影方法帶來的效果，可以使用一個非常酷的座標系投影介紹網站：worldmapcreator。在worldmapcreator中，我們可以選擇不同的投影方法來實時顯示投影變換效果，以及自定義投影結果，具體如下：

3.2 EPSG編碼

EPSG：European Petroleum Survey Group（歐洲石油調查小組）是一個涉及測地學、測量、製圖學與石油勘探相關的科學組織，它維護和釋出的EPSG編碼系統提供了各種座標參照系統crs的資料集引數。我們通過如下網站查詢EPSG投影座標的編碼：espg.io，該網站介面介紹如下：

網際網路常見的座標參考系如下圖所示。WGS84和CGCS2000都為地理座標系統。其中WGS84（World Geodetic System 1984）是目前應用最為廣泛的座標系統。在沒有特別宣告的時候，預設使用WGS84座標系統儲存地圖資料或提供地圖資料。WGS84座標單位為度，WGS84的EPSG程式碼為EPSG:4326。CGCS2000是我國當前最新的國家大地座標系，它的EPSG編碼為EPSG:4490，在大多數精度要求不高的情況下，我們可以預設WGS84座標系等同於CGCS2000座標系。

除此之外，WGS84 Web Mercator（EPSG:3857）座標系統也較為常用。EPSG:3857為投影座標系統，單位為米。EPSG:3857在WGS84座標系基礎上進行偽墨卡託投影（Pseudo-Mercator）得到。關於Pseudo-Mercator介紹見墨卡託投影。墨卡託投影能夠保證方向的正確，所以EPSG:3857廣泛用於導航海航。但越到高緯度，墨卡託投影會導致大小扭曲越嚴重，到兩極會被放到無限大。在長度表示上EPSG:3857有很大偏差，無法用於實際距離和麵積的測算。墨卡託投影下每個國家的大小和實際大小的差異見下圖：

當然CRS還有其他常見的資料編碼格式。計算機領域更多地使用EPSG編碼。

3.3 在GeoPandas中管理座標參考系

GeoPandas提供了豐富的函數介面以管理座標參考系，具體使用如下：

3.3.1 座標參考系管理

檢視座標管理系

# 框選資料
import geopandas as gpd

# 讀取江蘇省資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")

# 檢視座標管理系，可以看到預設座標參考系為EPSG:4326，單位為度。
data.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# 計算面積，可以看到面積計算結果不對，因為單位是度
# 直接計算面積會報錯，提示UserWarning: Geometry is in a geographic CRS
sum(data.area)

9.997823062190006

# 繪圖
data.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6e5aed0>

設定座標參考系

對於有些幾何形狀或者座標點沒有設定座標參考系，為了告訴GeoPandas如何解釋這些幾何形狀，需要設定預設座標系。GeoPandas提供set_crs函數來執行這一操作。

p = gpd.GeoSeries([geometry.Point([134.1233,23.456])])
# 沒有預設座標系
p.crs

# 重新設定座標系
p = p.set_crs("EPSG:4326")
p.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

重設座標參考系

GeoPandas提供了crs函數以重新設定座標參考系。

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

world.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6d79d90>

# 將座標參考系轉換為EPSG:3857
# 這裡捨去北極圈和南極圈資料，在這些地區WGS84/ Pseudo-Mercator變形嚴重
world = world.to_crs('EPSG:3857')
world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
world.crs

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# 橫縱座標尺度變為米
world.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6d65c10>

3.3.2 選擇合適的投影座標參考系

座標參考系選擇方法

如果我們獲得是地理座標系資料，如WGS84，在沒有更多資訊的情況下，如何轉換到投影座標參考系。一個好的辦法就是根據所在區域的幾何中心點，從前面提到的espg.io尋找合適EPSG編碼。這種方法很不專業，只是低精度使用！

以投影江蘇省地圖為例，一種辦法是找到整個國家的幾何中心，一種是找到當地的幾何中心。

# 框選資料
import geopandas as gpd

# 讀取整個國家的資料，無子區域資料
data_china = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/100000.json")

# 檢視經緯度幾何中心，POINT (103.9 36.4)
# 東經103.9
data_china.centroid

0    POINT (103.9 36.4)
dtype: geometry

# 框選資料
import geopandas as gpd

# 讀取江蘇省資料，無子區域資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000.json")

# 檢視經緯度幾何中心，POINT (119.5 33.0)
# 東經119.5
data.centroid

0    POINT (119.5 33.0)
dtype: geometry

接下來我們可以在 espg.io中，搜尋New Beijing kind:PROJCRS。這個搜尋詞表示搜尋以New Beijing為關鍵詞，投影座標(kind:PROJCRS)的結果。當然以beijing kind:PROJCRS和xian kind:PROJCRS也是可以的。這些都是我國常用的座標系，不新增kind:PROJCRS也沒關係，搜尋結果Coordinate reference systems中有相應的選項。然後我們檢視地區的幾何中心在是否座標系簡介的Area of use中。如下圖所示：

比如以整個中國的幾何中心(103.9 36.4)查詢的結果為EPSG:4573，如下所示：

以江蘇省的幾何中心(119.5 33.0)查詢的結果為EPSG:4586，如下所示：

座標參考系轉換範例

江蘇省的面積為10.72萬平方公里，通過EPSG:4573計算的面積為10.84萬平方公里，EPSG:4586計算的面積為10.38萬平方公里。這兩種都是近似估算，實際不會這麼粗糙。但是可以看到選擇合適的EPSG編碼結果都大差不差，想要更準確的結果，可以自行查詢其他EPSG編碼。

此外，為了展示結果可靠性 我們可以獲得某一座標位置進行繪製。通過地圖查詢蘇州市上的某一經緯度座標(120.601868 31.332525)，如下所示:

# 框選資料
import geopandas as gpd
from shapely import geometry

# 讀取江蘇省資料
data = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/320000_full.json")


ax = data.plot()
# 繪製點，該座標點需要預先設定crs
p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.plot(ax=ax, color='red', markersize=100)

<matplotlib.axes._subplots.AxesSubplot at 0x7fbbc6c62550>

# EPSG:4573下計算面積，EPSG:4573單位尺度為米
# 1平方公里等於10的6次方米，進行計算後結果為10.84萬平方公里
data.to_crs("EPSG:4573").area.sum()/1e6/1e4

10.849724500925806

ax = data.to_crs("EPSG:4573").plot()

p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4573").plot(ax=ax, color='red', markersize=100)
p.to_crs("EPSG:4573")

0    POINT (19993305.0 3575282.1)
dtype: geometry

# EPSG:4586下計算面積
data.to_crs("EPSG:4586").area.sum()/1e6/1e4

10.381446380451157

# 注意繪圖結果是有差別的
ax = data.to_crs("EPSG:4586").plot()

p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4586").plot(ax=ax, color='red', markersize=100)
p.to_crs("EPSG:4586")

0    POINT (842904.5 3473513.1)
dtype: geometry

從上面可以看到EPSG:4573和EPSG:4586表示地圖的橫座標單位不一樣。EPSG:4573代表的是New Beijing / Gauss-Kruger zone 18，這個zone就是指帶號的意思。EPSG:4586代表的是New Beijing / Gauss-Kruger CM 117E，表示沒有加帶號。簡單解釋：例如常見的6度分帶投影，表示從零度子午線開始，自西向東每個經差6度為一投影帶，全球共分60個帶，用1，2，3，4，5，... ，60表示。比如沒帶號的橫座標為500000米。有帶號zone 18，表示在沒帶號的橫座標前加上帶號18，變為18500000。具體可以查閱投影分帶。

在本例中，如果以江蘇省的幾何中心查詢加帶號的座標系，則應該是New Beijing / Gauss-Kruger zone 20（EPSG:4575）。與EPSG:4586相比，以EPSG:4575為基準繪圖，橫座標加了帶號（20），面積計算結果和是相等的。這些只是GIS基礎知識，本文說得非常簡單也有諸多錯誤（因為作者非專業人士），低精度軟體繪製地圖也不需要用到過多專業的知識。如果想高精度繪圖，請從頭開始學習GIS相關知識。

data.to_crs("EPSG:4575").area.sum()/1e6/1e4

10.381446380451155

ax = data.to_crs("EPSG:4575").plot()
p = gpd.GeoSeries([geometry.Point([120.601868,31.332525])]).set_crs("EPSG:4326")
p.to_crs("EPSG:4575").plot(ax=ax, color='red', markersize=100)

# EPSG:4586 為POINT (842904.5 3473513.1)
p.to_crs("EPSG:4575")

0    POINT (20842904.5 3473513.1)
dtype: geometry

4 參考

4.1 GeoPandas使用

• GeoPandas
• GeoPandas-doc
• DataV.GeoAtlas
• espg.io
• EPSG:4573
• EPSG:4586

4.2 相關資料

• Coordinate Reference System and Spatial Projection
• 聊聊GIS中的座標系
• 基於geopandas的空間資料分析——座標參考系篇
• 前後端都要懂一點的 GIS
• GIS基礎知識
• worldmapcreator
• 墨卡託投影
• 投影分帶