首页 > 技术文章 > kaggle地理空间分析-proximity analysis(邻近度分析)<v>

hy1231 2020-09-01 20:04 原文

介绍

在本教程中,您将探索几种用于邻近度分析的技术。 特别是,您将学习如何执行以下操作:

  • 测量地图上点之间的距离,以及
  • 选择某要素半径内的所有点。
import folium
from folium import Marker, GeoJson
from folium.plugins import HeatMap

import pandas as pd
import geopandas as gpd

您将使用来自美国环境保护署(EPA)的数据集,该数据集跟踪美国宾夕法尼亚州费城的有毒化学物质的释放。

releases = gpd.read_file("../input/geospatial-learn-course-data/toxic_release_pennsylvania/toxic_release_pennsylvania/toxic_release_pennsylvania.shp") 
releases.head()


您还将使用包含来自同一城市的空气质量监测站的读数的数据集。

stations = gpd.read_file("../input/geospatial-learn-course-data/PhillyHealth_Air_Monitoring_Stations/PhillyHealth_Air_Monitoring_Stations/PhillyHealth_Air_Monitoring_Stations.shp")
stations.head()


Measuring distance(测量距离)

要测量来自两个不同GeoDataFrame的点之间的距离,我们首先必须确保它们使用相同的坐标参考系统(CRS)。 幸运的是,这里就是这种情况,它们都使用EPSG 2272。

print(stations.crs)
print(releases.crs)
{'init': 'epsg:2272'}
{'init': 'epsg:2272'}

我们还检查CRS,以查看其使用的单位(米,英尺或其他)。 在这种情况下,EPSG 2272具有英尺单位。 (如果您愿意,可以在此处进行检查。)

在GeoPandas中计算距离相对简单。 下面的代码单元计算出last_release中相对较新的版本事件与stationGeoDataFrame中的每个站点之间的距离(以英尺为单位)。

# Select one release incident in particular
recent_release = releases.iloc[360]

# Measure distance from release to each station
distances = stations.geometry.distance(recent_release.geometry)
distances
out:
0     44778.509761
1     51006.456589
2     77744.509207
3     14672.170878
4     43753.554393
5      4711.658655
6     23197.430858
7     12072.823097
8     79081.825506
9      3780.623591
10    27577.474903
11    19818.381002
dtype: float64

使用计算出的距离,我们可以获得诸如到每个站点的平均距离之类的统计信息。

print('Mean distance to monitoring stations: {} feet'.format(distances.mean()))
Mean distance to monitoring stations: 33516.28487007786 feet

或者,我们可以找到最近的监测站。

print('Closest monitoring station ({} feet):'.format(distances.min()))
print(stations.iloc[distances.idxmin()][["ADDRESS", "LATITUDE", "LONGITUDE"]])
Closest monitoring station (3780.623590556444 feet):
ADDRESS      3100 Penrose Ferry Road
LATITUDE                     39.9128
LONGITUDE                   -75.1854
Name: 9, dtype: object

Creating a buffer(创建一个缓冲区)

如果我们想了解地图上所有距某个点一定半径的点,最简单的方法是创建一个缓冲区。

下面的代码单元创建一个GeoSeries two_mile_buffer,其中包含12个不同的Polygon对象。 每个多边形是在不同的空气监测站周围2英里(或2 * 5280英尺)的缓冲区。

two_mile_buffer = stations.geometry.buffer(2*5280)
two_mile_buffer.head()
0    POLYGON ((2721944.640797138 257149.3104284704,...
1    POLYGON ((2682494.289907977 271248.9000113755,...
2    POLYGON ((2744886.638220146 280980.2466903776,...
3    POLYGON ((2703638.579968393 233247.1013432145,...
4    POLYGON ((2726959.772827223 251134.9763285518,...
dtype: object
我们使用`folium.GeoJson()`在地图上绘制每个多边形。 请注意,由于filium需要纬度和经度的坐标,因此在绘制之前必须将CRS转换为EPSG 4326。

Create map with release incidents and monitoring stations

m = folium.Map(location=[39.9526,-75.1652], zoom_start=11)
HeatMap(data=releases[['LATITUDE', 'LONGITUDE']], radius=15).add_to(m)
for idx, row in stations.iterrows():
Marker([row['LATITUDE'], row['LONGITUDE']]).add_to(m)

Plot each polygon on the map

GeoJson(two_mile_buffer.to_crs(epsg=4326)).add_to(m)

Show the map

m

现在,要测试在**任何**监测站的2英里范围内是否发生了有毒物质释放,我们可以对每个多边形运行12种不同的测试(以单独检查其是否包含该点)。

但是更有效的方法是先将所有多边形折叠成一个**MultiPolygon**对象。 我们使用`unary_union`属性来执行此操作。

Turn group of polygons into single multipolygon

my_union = two_mile_buffer.geometry.unary_union
print('Type:', type(my_union))

Show the MultiPolygon object

my_union

Type: <class 'shapely.geometry.multipolygon.MultiPolygon'>

![](https://img2020.cnblogs.com/blog/2043753/202009/2043753-20200902143402462-321060605.png)

我们使用`contains()`方法来检查多边形是否包含一个点。 我们将使用本教程前面部分中的发布事件,我们知道它离最近的监视站大约3781英尺。

The closest station is less than two miles away

my_union.contains(releases.iloc[360].geometry)

`out:True`
但是,并非所有的释放都发生在空气监测站的两英里范围内!

The closest station is more than two miles away

my_union.contains(releases.iloc[358].geometry)

`out:False`

# Your turn
在[最后的练习](https://www.kaggle.com/kernels/fork/5832147)中,您将调查纽约市的医院覆盖率。

推荐阅读