GeoPandas如何實現坐標參考系統

這篇文章主要介紹了 GeoPandas如何實現坐標參考系統，具有一定借鑒價值，感興趣的朋友可以參考下，希望大家閱讀完這篇文章之后大有收獲，下面讓小編帶著大家一起了解一下。

02-坐標參考系統

源代碼請看此處

%matplotlib inline

import pandas as pd
import geopandas

countries = geopandas.read_file("zip://data/ne_110m_admin_0_countries.zip")
cities = geopandas.read_file("zip://data/ne_110m_populated_places.zip")
rivers = geopandas.read_file("zip://data/ne_50m_rivers_lake_centerlines.zip")

2.1 坐標參考系統

到目前為止，我們已經使用了具有特定坐標的幾何數據，而沒有進一步想知道這些坐標是什么意思或者它們是如何表達的

**坐標參考系統(CRS)**將坐標與地球上的特定位置相關聯

有關坐標參考系統的詳細解釋：https://docs.qgis.org/2.8/en/docs/gentle_gis_introduction/coordinate_reference_systems.html

2.1.1 地理坐標系

即經緯度

例如 48°51′N, 2°17′E

最常見的坐標類型是地理坐標，緯度和經度來定義地球上相對于赤道和本初子午線的位置

有了這個坐標系統，我們可以很容易地指定地球上的任何位置，例如在全球定位系統中，如果在谷歌地圖上查看一個位置的坐標，會看到緯度和經度

注意 Python中的使用(緯度，經度)(lon,lat)表示地理坐標

• Longitude: [-180, 180]

• Latitude: [-90, 90]

2.1.2 投影坐標系

(x,y)坐標通常以米或英尺為單位

雖然地球是一個球體，但實際上我們通常在一個平面上表示它，如生活中的地圖，或者我們用Python在電腦屏幕上制作的地圖

從立體地球到平面地圖就是我們所說的投影

<table><tr> <td> <img src="img/projection.png"/> </td> </tr></table>

我們把地球的表面投影到2D平面上，這樣我們就可以在平面上用笛卡爾的x和y坐標表示位置。在這個平面上，我們通常使用長度單位，如米，而不是度，這使得分析更加方便和有效

然而，有一點很重要：三維地球永遠不可能完美地呈現在二維地圖上，所以投影不可避免地會引入失真。為了最大限度地減少這種錯誤，有不同的投影方法，每種方法都有特定的優點和缺點

一些投影系統將試圖保持幾何圖形的面積大小，如艾伯斯等面積投影（Albers Equal Area）。其他投影系統試圖保留角度，如墨卡托投影，但會看到該地區的大變形。每個投影系統總會有一些面積、角度或距離的失真

<table><tr> <td> <img src="img/projections-AlbersEqualArea.png"/>AlbersEqualArea </td> <td> <img src="img/projections-Mercator.png"/>Mercator </td> </tr> <tr> <td> <img src="img/projections-Robinson.png"/>Robinson </td> </tr></table>

投影尺寸與實際尺寸的比較 <table><tr> <td> <img src="https://github.com/jorisvandenbossche/geopandas-tutorial/blob/master/img/mercator_projection_area.gif?raw=true"/> </td> </tr></table>

2.2 Python和GeoPandas中的坐標參考系統

GeoDataFrame或GeoSeries具有一個.crs屬性，該屬性保存(可選)幾何坐標參考系統的描述:

countries.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
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

可以看出，對于countries數據框，使用的是EPSG 4326 / WGS84 lon/lat參考系統，這是地理坐標系中最常用的系統之一

它使用坐標作為緯度和經度的度數，從圖上的x/y標簽可以看出:

countries.plot()

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

.crs屬性的返回值是字典類型，在這種情況下，它只指示EPSG代碼，但它也可以包含完整的proj4字符串(以字典的形式)

可能的CRS表達：

• proj4字符串<br> 例如：+proj=longlat +datum=WGS84 +no_defs或者{'proj': 'longlat', 'datum': 'WGS84', 'no_defs': True}

• EPSG 代碼<br> 例如：EPSG:4326 = WGS84地理坐標系(longitude, latitude)

• 富文本（Well-Know-Text，WKT）表示（在下一個GeoPandas版本中，PROJ6提供了更好的支持）<br> 例如：https://epsg.io/4326

在后端，GeoPandas使用pyproj/PROJ庫來處理重投影

更多信息：http://geopandas.readthedocs.io/en/latest/projections.html

2.3 坐標系統的轉換

在GeoDataFrame中，使用to_crs函數進行坐標系統的轉換

例如，將countries轉換城墨卡托投影(http://epsg.io/3395):

# 因為墨卡托投影無法處理兩極，因此移除南極
countries=countries[(countries['name']!="Antarctica")]

countries_mercator=countries.to_crs(epsg=3395)# 或者使用 .to_crs({'init': 'epsg:3395'})

countries_mercator.plot()

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

注意比較此處的x和y軸刻度與地理坐標系中投影的區別

2.4 使用不同CRS的原因

• 具有不同CRS的不同源數據->需要轉換到相同的CRS下

  df1 = geopandas.read_file(...)
  df2 = geopandas.read_file(...)
  
  df2 = df2.to_crs(df1.crs)

  
  

• 地圖映射(由于形狀和距離的扭曲)

• 基于距離/面積的計算->確保使用適當的投影坐標系，以有意義的單位表示，如米或英尺(而不是度)

<div class="alert alert-info" >

注意:

所有發生在地表上的計算都假設數據是在2D笛卡爾平面上，因此這些計算的結果只有在數據被正確投影的情況下才是正確的

</div>

2.5 練習

再次使用巴黎數據集，到目前為止，我們為練習提供了一個合適的投影CRS中的數據集

但是原始數據實際上是地理坐標

在下面的練習中，我們將從原始數據開始

原始數據集，現在以地理坐標系的GeoJSON("data/paris_districts.geojson")文件提供

為了轉換為投影坐標，我們將使用法國的標準投影坐標系統，即RGF93 / Lambert-93參考系，參考號為EPSG:2154(比利時為 Lambert 72，EPSG為31370)。

2.5.1 練習一：對GeoDataFrame進行投影

• 將地區數據集("data/paris_districts.geojson")讀入名為districts的GeoDataFrame中

• 查看districts的CRS屬性

• 對districts數據集進行簡單繪圖

• 計算所有地區的面積

• 將districts進行投影(法國使用EPSG:2154)，將新數據集稱為districts_RGF93

• 繪制一個類似的districts_RGF93圖

• 對districts_RGF93再次計算所有區的面積(結果現在用m2表示)

districts=geopandas.read_file("data/paris_districts.geojson")

districts.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
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

districts.plot(figsize=(12,6))

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

districts.geometry.area

0     0.000107
1     0.000051
2     0.000034
3     0.000033
4     0.000023
        ...   
75    0.000159
76    0.000099
77    0.000182
78    0.000196
79    0.000256
Length: 80, dtype: float64

districts_RGF93=districts.to_crs(epsg=2154)

districts_RGF93.plot(figsize=(12,6))

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

districts_RGF93.geometry.area

0     8.690007e+05
1     4.124585e+05
2     2.736968e+05
3     2.694568e+05
4     1.880122e+05
          ...     
75    1.294988e+06
76    8.065686e+05
77    1.486971e+06
78    1.599002e+06
79    2.090904e+06
Length: 80, dtype: float64

可以看到在不同坐標系統下同一地區的面積值不同

在01-地理數據介紹中，我們做了一個練習，繪制了巴黎自行車站的位置，并使用contextily包為其添加了背景地圖

contextily”假設數據坐標系統為網絡墨卡托投影(Web Mercator projection)中，這是大多數網絡地圖服務使用的系統。在第一個練習中，我們在適當的CRS中提供了數據，所以不需要關心這個方面

但是，通常情況下，數據的參考坐標系不會是網絡墨卡托投影(Web Mercator projection)，因此我們必須自己將它們與網絡背景圖進行對齊

2.5.2 練習二：對GeoDataFrame進行投影

• 將自行車站數據集(data/paris_bike_stations.geojson)讀入名為stations的GeoDataFrame

• 將stations數據集轉換為網絡墨卡托投影(EPSG:3857)命名為stations_webmercator

• 繪制該投影數據集地圖圖(指定標記大小為5)，并使用contextily添加背景地圖

stations=geopandas.read_file('data/paris_bike_stations.geojson')

stations_webmercator=stations.to_crs(epsg=3857)

stations_webmercator.head()

<div> <style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style> <table border="1" class="dataframe"> <thead> <tr > <th></th> <th>name</th> <th>bike_stands</th> <th>available_bikes</th> <th>geometry</th> </tr> </thead> <tbody> <tr> <th>0</th> <td>14002 - RASPAIL QUINET</td> <td>44</td> <td>4</td> <td>POINT (259324.887 6247620.771)</td> </tr> <tr> <th>1</th> <td>20503 - COURS DE VINCENNES PYRéNéES</td> <td>21</td> <td>3</td> <td>POINT (267824.377 6249062.894)</td> </tr> <tr> <th>2</th> <td>20011 - PYRéNéES-DAGORNO</td> <td>21</td> <td>0</td> <td>POINT (267742.135 6250378.469)</td> </tr> <tr> <th>3</th> <td>31008 - VINCENNES (MONTREUIL)</td> <td>56</td> <td>0</td> <td>POINT (271326.638 6250750.824)</td> </tr> <tr> <th>4</th> <td>43006 - MINIMES (VINCENNES)</td> <td>28</td> <td>27</td> <td>POINT (270594.689 6248007.705)</td> </tr> </tbody> </table> </div>

import contextily as ctx

ax = stations_webmercator.plot(figsize=(10, 10), markersize=5, alpha=0.5, edgecolor='k')
ctx.add_basemap(ax, source=ctx.providers.Stamen.TonerLite)

感謝你能夠認真閱讀完這篇文章，希望小編分享的“ GeoPandas如何實現坐標參考系統”這篇文章對大家有幫助，同時也希望大家多多支持億速云，關注億速云行業資訊頻道，更多相關知識等著你來學習!