03-空间关系和操作
源代码 请看此处
%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")
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 image
143 image, extent = bounds2img(
--> 144 left, bottom, right, top, zoom=zoom, source=source, ll=False
145 )
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.z
247 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.cache
302 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=timeout
450 )
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-blocking
816 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-6
return _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>
917
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
