02-坐标参考系统

源代码 请看此处

%matplotlib inlineimport 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" style="font-size:120%">

注意:

所有发生在地表上的计算都假设数据是在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 再次计算所有区的面积(结果现在用m²表示)

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 style="text-align: right;"> <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 ctxax = stations_webmercator.plot(figsize=(10, 10), markersize=5, alpha=0.5, edgecolor='k')
ctx.add_basemap(ax, source=ctx.providers.Stamen.TonerLite)