03-空间关系和操作

源代码 请看此处

%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")

3.1 空间关系

地理空间数据的一个重要方面是我们可以观察 spatial relationships ：两个空间对象如何相互关联(它们是否重叠、相交、包含……)

地理信息系统中的拓扑关系、集合论关系通常基于DE-9IM模型，更多信息请见 https://en.wikipedia.org/wiki/Spatial_relation

(Image by Krauss, CC BY-SA 3.0 )

3.2 单个对象之间的空间关系

首先创建一些空间对象:

多边形：使用 .squeeze() 从 GeoSeries 中提取几何对象向量

belgium=countries.loc[countries['name']=='Belgium','geometry']#.squeeze()
# .squeeze()从数组的形状中删除单维度条目，即把shape中为1的维度去掉

type(belgium)

geopandas.geoseries.GeoSeries

belgium=countries.loc[countries['name']=='Belgium','geometry'].squeeze()

type(belgium)

shapely.geometry.polygon.Polygon

belgium

两个点对象

paris=cities.loc[cities['name']=='Paris','geometry'].squeeze()
brussels = cities.loc[cities['name'] == 'Brussels', 'geometry'].squeeze()

paris

brussels

线段

from shapely.geometry import LineString
line = LineString([paris, brussels])

line

把这4个几何对象一起可视化（把它们放在一个GeoSeries中，以便于用 geopandas.plot() 方法把它们一起显示出来）

geopandas.GeoSeries([belgium,paris,brussels,line]).plot(cmap='tab10')

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

可以认出比利时的抽象形状

比利时首都布鲁塞尔位于比利时境内。这是一种空间关系，我们可以使用以下单个的几何对象进行测试

brussels.within(belgium)

True

belgium.contains(brussels)

True

belgium.contains(paris)

False

paris.within(belgium)

False

从巴黎到布鲁塞尔所画的直线并不完全位于比利时境内，但它确实与比利时相交

belgium.contains(line)

False

line.intersects(belgium)

True

3.3 GeoDataFrame中的空间关系

如上所述，在单个 shapely 矢量对象上可用的方法，也可作为 GeoSeries / GeoDataFrame 对象的方法

例如，如果我们在世界数据集上调用带有 paris 点的 contains 方法，它将对 world 数据框架中的每个国家进行这种空间检查

countries.contains(paris)

0      False
1      False
2      False
3      False
4      False...
172    False
173    False
174    False
175    False
176    False
Length: 177, dtype: bool

因为上面给了我们一个布尔结果，我们可以用它来筛选GeoDataFrame

countries[countries.contains(paris)]

<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>iso_a3</th> <th>name</th> <th>continent</th> <th>pop_est</th> <th>gdp_md_est</th> <th>geometry</th> </tr> </thead> <tbody> <tr> <th>55</th> <td>FRA</td> <td>France</td> <td>Europe</td> <td>67106161.0</td> <td>2699000.0</td> <td>MULTIPOLYGON (((2.51357 51.14851, 2.65842 50.7...</td> </tr> </tbody> </table> </div>

事实上，法国是世界上唯一一个巴黎所在的国家

又如，提取南美洲亚马逊河的线路图，我们可以查询河流流经哪些国家。

amazon = rivers[rivers['name'] == 'Amazonas'].geometry.squeeze()

amazon

countries[countries.crosses(amazon)]

<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>iso_a3</th> <th>name</th> <th>continent</th> <th>pop_est</th> <th>gdp_md_est</th> <th>geometry</th> </tr> </thead> <tbody> <tr> <th>22</th> <td>BRA</td> <td>Brazil</td> <td>South America</td> <td>207353391.0</td> <td>3081000.0</td> <td>POLYGON ((-57.62513 -30.21629, -56.29090 -28.8...</td> </tr> <tr> <th>35</th> <td>COL</td> <td>Colombia</td> <td>South America</td> <td>47698524.0</td> <td>688000.0</td> <td>POLYGON ((-66.87633 1.25336, -67.06505 1.13011...</td> </tr> <tr> <th>124</th> <td>PER</td> <td>Peru</td> <td>South America</td> <td>31036656.0</td> <td>410400.0</td> <td>POLYGON ((-69.52968 -10.95173, -68.66508 -12.5...</td> </tr> </tbody> </table> </div>

countries[countries.intersects(amazon)]

<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>iso_a3</th> <th>name</th> <th>continent</th> <th>pop_est</th> <th>gdp_md_est</th> <th>geometry</th> </tr> </thead> <tbody> <tr> <th>22</th> <td>BRA</td> <td>Brazil</td> <td>South America</td> <td>207353391.0</td> <td>3081000.0</td> <td>POLYGON ((-57.62513 -30.21629, -56.29090 -28.8...</td> </tr> <tr> <th>35</th> <td>COL</td> <td>Colombia</td> <td>South America</td> <td>47698524.0</td> <td>688000.0</td> <td>POLYGON ((-66.87633 1.25336, -67.06505 1.13011...</td> </tr> <tr> <th>124</th> <td>PER</td> <td>Peru</td> <td>South America</td> <td>31036656.0</td> <td>410400.0</td> <td>POLYGON ((-69.52968 -10.95173, -68.66508 -12.5...</td> </tr> </tbody> </table> </div>

空间关系函数 :

检查空间关系的不同函数概述（ 空间谓词函数 ）

• equals
• contains
• crosses
• disjoint
• intersects
• overlaps
• touches
• within
• covers

关于这些方法的概述，见 https://shapely.readthedocs.io/en/stable/manual.html#predicates-and-relationships

关于这些操作的语义，详见 https://en.wikipedia.org/wiki/DE-9IM

3.3 练习

我们将再次使用巴黎的数据集来做一些练习，再次导入

districts = geopandas.read_file("data/paris_districts.geojson").to_crs(epsg=2154)
stations = geopandas.read_file("data/paris_bike_stations.geojson").to_crs(epsg=2154)

3.3.1 练习一：埃菲尔铁塔

埃菲尔铁塔是一座建于19世纪的铁格子塔，可能是巴黎最具代表性的景观

埃菲尔铁塔的位置是：x：648237.3，y：6862271.9

• 以埃菲尔铁塔的坐标创建一个 Shapely 点对象，并将其分配给一个名为 eiffel_tower 的变量，输出结果
• 检查埃菲尔铁塔是否位于Montparnasse区
• 检查 Montparnasse 区是否包含自行车站
• 计算埃菲尔铁塔和自行车站之间的距离（注意：在这种情况下，距离以米为单位返回）。

from shapely.geometry import Point

eiffel_tower=Point(648237.3,6862271.9)

print(eiffel_tower)

POINT (648237.3 6862271.9)

district_montparnasse=districts[districts['district_name']=='Montparnasse']['geometry'].squeeze()

district_montparnasse

# 检查埃菲尔铁塔是否位于Montparnasse区
district_montparnasse.contains(eiffel_tower)

False

bike_station = stations.loc[293, 'geometry']

bike_station

# 检查 Montparnasse 区是否包含自行车站
district_montparnasse.contains(bike_station)

True

# 计算埃菲尔铁塔和自行车站之间的距离
eiffel_tower.distance(bike_station)

3540.1534488921966

3.3.2 练习二：埃菲尔铁塔位于哪个区？

在上一个练习中，我们为埃菲尔铁塔的位置构造了一个 Point 对象，我们检查了它不在montparnasse蒙帕纳斯区，现在我们的问题是确定它位于巴黎的哪个区

• 创建一个布尔掩模（或筛选器），表示每个区是否包含埃菲尔铁塔。调用结果 mask 。
• 用布尔掩码筛选 districts 数据框并打印结果。

mask=districts.contains(eiffel_tower)

districts[mask]

<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>id</th> <th>district_name</th> <th>population</th> <th>geometry</th> </tr> </thead> <tbody> <tr> <th>27</th> <td>28</td> <td>Gros-Caillou</td> <td>25156</td> <td>POLYGON ((649336.752 6861767.761, 649110.815 6...</td> </tr> </tbody> </table> </div>

3.3.3 练习三：找出距离埃菲尔铁塔最近的自行车站点

现在，我们可能对埃菲尔铁塔附近的自行车站感兴趣。为了探索它们，让我们将埃菲尔铁塔本身以及1公里内的自行车站可视化。

要做到这一点，我们可以计算出每个站点到埃菲尔铁塔的距离。根据这个结果，我们就可以创建一个掩模，如果站点在1km以内，则取 True ，否则取 False ，然后用它来过滤站点的GeoDataFrame。最后，我们对这个子集进行可视化处理。

• 计算每个站点到埃菲尔铁塔的距离，并调用结果 dist_eiffel
• 打印距离最近的车站的距离（这是 dist_eiffel 的最小值）
• 选择距离埃菲尔铁塔小于1公里的 "站点 "GeoDataFrame行（注意，距离是以米为单位）。调用结果 stations_eiffel

dist_eiffel=stations.distance(eiffel_tower)

print(dist_eiffel.min())

232.34672323449053

stations_eiffel=stations[dist_eiffel<1000]

ax = stations_eiffel.to_crs(epsg=3857).plot()geopandas.GeoSeries([eiffel_tower], crs='EPSG:2154').to_crs(epsg=3857).plot(ax=ax, color='red')

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

import contextily
contextily.add_basemap(ax)

---------------------------------------------------------------------------KeyboardInterrupt                         Traceback (most recent call last)<ipython-input-43-2923a97d25c0> in <module>1 import contextily
----> 2 contextily.add_basemap(ax)/usr/local/lib/python3.6/dist-packages/contextily/plotting.py in add_basemap(ax, zoom, source, interpolation, attribution, attribution_size, reset_extent, crs, resampling, url, **extra_imshow_args)142         # Download image143         image, extent = bounds2img(
--> 144             left, bottom, right, top, zoom=zoom, source=source, ll=False145         )146         # Warping/usr/local/lib/python3.6/dist-packages/contextily/tile.py in bounds2img(w, s, e, n, zoom, source, ll, wait, max_retries, url)246         x, y, z = t.x, t.y, t.z247         tile_url = _construct_tile_url(provider, x, y, z)
--> 248         image = _fetch_tile(tile_url, wait, max_retries)249         tiles.append(t)250         arrays.append(image)/usr/local/lib/python3.6/dist-packages/joblib/memory.py in __call__(self, *args, **kwargs)563 564     def __call__(self, *args, **kwargs):
--> 565         return self._cached_call(args, kwargs)[0]566 567     def __getstate__(self):/usr/local/lib/python3.6/dist-packages/joblib/memory.py in _cached_call(self, args, kwargs, shelving)529 530         if must_call:
--> 531             out, metadata = self.call(*args, **kwargs)532             if self.mmap_mode is not None:533                 # Memmap the output at the first call to be consistent with/usr/local/lib/python3.6/dist-packages/joblib/memory.py in call(self, *args, **kwargs)725         if self._verbose > 0:726             print(format_call(self.func, args, kwargs))
--> 727         output = self.func(*args, **kwargs)728         self.store_backend.dump_item(729             [func_id, args_id], output, verbose=self._verbose)/usr/local/lib/python3.6/dist-packages/contextily/tile.py in _fetch_tile(tile_url, wait, max_retries)301 @memory.cache302 def _fetch_tile(tile_url, wait, max_retries):
--> 303     request = _retryer(tile_url, wait, max_retries)304     with io.BytesIO(request.content) as image_stream:305         image = Image.open(image_stream).convert("RGB")/usr/local/lib/python3.6/dist-packages/contextily/tile.py in _retryer(tile_url, wait, max_retries)441     """442     try:
--> 443         request = requests.get(tile_url, headers={"user-agent": USER_AGENT})444         request.raise_for_status()445     except requests.HTTPError:/usr/local/lib/python3.6/dist-packages/requests/api.py in get(url, params, **kwargs)74 75     kwargs.setdefault('allow_redirects', True)
---> 76     return request('get', url, params=params, **kwargs)77 78 /usr/local/lib/python3.6/dist-packages/requests/api.py in request(method, url, **kwargs)59     # cases, and look like a memory leak in others.60     with sessions.Session() as session:
---> 61         return session.request(method=method, url=url, **kwargs)62 63 /usr/local/lib/python3.6/dist-packages/requests/sessions.py in request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)528         }529         send_kwargs.update(settings)
--> 530         resp = self.send(prep, **send_kwargs)531 532         return resp/usr/local/lib/python3.6/dist-packages/requests/sessions.py in send(self, request, **kwargs)641 642         # Send the request
--> 643         r = adapter.send(request, **kwargs)644 645         # Total elapsed time of the request (approximately)/usr/local/lib/python3.6/dist-packages/requests/adapters.py in send(self, request, stream, timeout, verify, cert, proxies)447                     decode_content=False,448                     retries=self.max_retries,
--> 449                     timeout=timeout450                 )451 /usr/local/lib/python3.6/dist-packages/urllib3/connectionpool.py in urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw)675                 body=body,676                 headers=headers,
--> 677                 chunked=chunked,678             )679 /usr/local/lib/python3.6/dist-packages/urllib3/connectionpool.py in _make_request(self, conn, method, url, timeout, chunked, **httplib_request_kw)379         # Trigger any extra validation we need to do.380         try:
--> 381             self._validate_conn(conn)382         except (SocketTimeout, BaseSSLError) as e:383             # Py2 raises this as a BaseSSLError, Py3 raises it as socket timeout./usr/local/lib/python3.6/dist-packages/urllib3/connectionpool.py in _validate_conn(self, conn)974         # Force connect early to allow us to validate the connection.975         if not getattr(conn, "sock", None):  # AppEngine might not have  `.sock`
--> 976             conn.connect()977 978         if not conn.is_verified:/usr/local/lib/python3.6/dist-packages/urllib3/connection.py in connect(self)368             ca_cert_data=self.ca_cert_data,369             server_hostname=server_hostname,
--> 370             ssl_context=context,371         )372 /usr/local/lib/python3.6/dist-packages/urllib3/util/ssl_.py in ssl_wrap_socket(sock, keyfile, certfile, cert_reqs, ca_certs, server_hostname, ssl_version, ciphers, ssl_context, ca_cert_dir, key_password, ca_cert_data)375     ) or IS_SECURETRANSPORT:376         if HAS_SNI and server_hostname is not None:
--> 377             return context.wrap_socket(sock, server_hostname=server_hostname)378 379         warnings.warn(/usr/lib/python3.6/ssl.py in wrap_socket(self, sock, server_side, do_handshake_on_connect, suppress_ragged_eofs, server_hostname, session)405                          suppress_ragged_eofs=suppress_ragged_eofs,406                          server_hostname=server_hostname,
--> 407                          _context=self, _session=session)408 409     def wrap_bio(self, incoming, outgoing, server_side=False,/usr/lib/python3.6/ssl.py in __init__(self, sock, keyfile, certfile, server_side, cert_reqs, ssl_version, ca_certs, do_handshake_on_connect, family, type, proto, fileno, suppress_ragged_eofs, npn_protocols, ciphers, server_hostname, _context, _session)815                         # non-blocking816                         raise ValueError("do_handshake_on_connect should not be specified for non-blocking sockets")
--> 817                     self.do_handshake()818 819             except (OSError, ValueError):/usr/lib/python3.6/ssl.py in do_handshake(self, block)1075             if timeout == 0.0 and block:1076                 self.settimeout(None)
-> 1077             self._sslobj.do_handshake()1078         finally:1079             self.settimeout(timeout)/usr/lib/python3.6/ssl.py in do_handshake(self)687     def do_handshake(self):688         """Start the SSL/TLS handshake."""
--> 689         self._sslobj.do_handshake()690         if self.context.check_hostname:691             if not self.server_hostname:KeyboardInterrupt:

3.4 空间操作

除了返回布尔值的空间谓词判断外，Shapely和GeoPandas还提供了返回新几何对象的操作

二元操作

<table><tr> <td> <img src="img/spatial-operations-base.png"/> </td> <td> <img src="img/spatial-operations-intersection.png"/> </td> </tr> <tr> <td> <img src="img/spatial-operations-union.png"/> </td> <td> <img src="img/spatial-operations-difference.png"/> </td> </tr></table>

缓冲区:

<table><tr> <td> <img src="img/spatial-operations-buffer-point1.png"/> </td> <td> <img src="img/spatial-operations-buffer-point2.png"/> </td> </tr> <tr> <td> <img src="img/spatial-operations-buffer-line.png"/> </td> <td> <img src="img/spatial-operations-buffer-polygon.png"/> </td> </tr></table>

详见 https://shapely.readthedocs.io/en/stable/manual.html#spatial-analysis-methods。

例如，使用上面的实验数据，在布鲁塞尔周围构造一个缓冲区（返回一个多边形）

geopandas.GeoSeries([belgium,brussels.buffer(1)]).plot(alpha=0.5,cmap='tab10')

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

现在取这两个多边形的交、并、差

brussels.buffer(1).intersection(belgium)

brussels.buffer(1).union(belgium)

brussels.buffer(1).difference(belgium)

另一个有用的属性是 .unary_union ，它通过取所有这些矢量对象的并集，将GeoDataFrame中的一组矢量对象转换为一个单一的矢量对象（即矢量数据的 融合 ）

例如，我们可以为非洲大陆构造一个单一对象

africa_countries = countries[countries['continent'] == 'Africa']

africa=africa_countries.unary_union

africa

print(str(africa)[:1000])

MULTIPOLYGON (((32.83012047702888 -26.7421916643362, 32.58026492689768 -27.47015756603182, 32.46213260267845 -28.30101124442056, 32.20338870619304 -28.75240488049007, 31.52100141777888 -29.25738697684626, 31.325561150851 -29.40197763439891, 30.90176272962535 -29.90995696382804, 30.62281334811382 -30.42377573010613, 30.05571618014278 -31.14026946383296, 28.92555260591954 -32.1720411109725, 28.2197558936771 -32.77195281344886, 27.46460818859597 -33.2269637997788, 26.41945234549283 -33.61495045342619, 25.90966434093349 -33.6670402971764, 25.7806282895007 -33.94464609144834, 25.17286176931597 -33.79685149509358, 24.67785322439212 -33.98717579522455, 23.59404340993464 -33.79447437920815, 22.98818891774474 -33.91643075941698, 22.57415734222224 -33.86408253350531, 21.54279910654103 -34.25883879978294, 20.689052768647 -34.41717538832523, 20.07126102059763 -34.79513681410799, 19.61640506356457 -34.81916635512371, 19.19327843595872 -34.46259897230979, 18.85531456876987 -34.44430551527847, 18.424

<div class="alert alert-info" style="font-size:120%">

注意

GeoPandas（和Shapely的单个对象）提供了大量的基本方法来分析地理空间数据（距离、长度、中心点、边界、凸壳、简化、变换......），比我们在本教程中能接触到的几种方法多得多。

• GeoPandas提供的所有方法的概述可以在这里找到： http://geopandas.readthedocs.io/en/latest/reference.html

</div>

3.5 空间操作练习

3.5.1 练习一：塞纳河流经哪些区

下面，以类似于GeoJSON的特征字典（创建于 http://geojson.io） 的方式提供了巴黎附近的塞纳河的坐标

根据这个 seine 对象，我们想知道哪些地区离塞纳河很近（最多150米）

• 在塞纳河周围建立一个150米的缓冲区
• 检查哪些区域与该缓冲对象相交
• 将各区的情况可视化，说明哪些区位于塞纳河附近

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

# 根据http://geojson.io的geojson格式创建一条曲线
s_seine = geopandas.GeoDataFrame.from_features({"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"LineString","coordinates":[[2.408924102783203,48.805619828930226],[2.4092674255371094,48.81703747481909],[2.3927879333496094,48.82325391133874],[2.360687255859375,48.84912860497674],[2.338714599609375,48.85827758964043],[2.318115234375,48.8641501307046],[2.298717498779297,48.863246707697],[2.2913360595703125,48.859519915404825],[2.2594070434570312,48.8311646245967],[2.2436141967773438,48.82325391133874],[2.236919403076172,48.82347994904826],[2.227306365966797,48.828339513221444],[2.2224998474121094,48.83862215329593],[2.2254180908203125,48.84856379804802],[2.2240447998046875,48.85409863123821],[2.230224609375,48.867989496547864],[2.260265350341797,48.89192242750887],[2.300262451171875,48.910203080780285]]}}]},crs={'init': 'epsg:4326'})
s_seine_utm = s_seine.to_crs(epsg=2154)

/usr/local/lib/python3.6/dist-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6return _prepare_from_string(" ".join(pjargs))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(16, 7))
districts.plot(ax=ax, color='grey', alpha=0.4, edgecolor='k')
s_seine_utm.plot(ax=ax)

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

seine=s_seine_utm.geometry.squeeze()

seine_buffer=seine.buffer(150)

seine_buffer

districts_seine = districts[districts.intersects(seine_buffer)]

districts_seine.plot()

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

fig, ax = plt.subplots(figsize=(16, 7))
districts.plot(ax=ax,color='grey',alpha=0.4,edgecolor='k')
districts_seine.plot(ax=ax,color='red',alpha=0.4,edgecolor='k')
s_seine_utm.plot(ax=ax,color='blue')

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