地理空间分析

1、你的第一张地图

介绍

在本微型课程中，您将学习不同的方法来处理和可视化地理空间数据或具有地理位置的数据。

在此过程中，您将为以下几个实际问题提供解决方案：

一个全球性的非营利组织应该在哪里扩大其在菲律宾偏远地区的影响力？

紫色马丁鸟，一种濒危鸟类，是如何在北美和南美之间旅行的？这些鸟是去保护区的吗？

日本哪些地区可能从额外的地震加固中受益？

加利福尼亚州的哪些星巴克商店是星巴克下一个保留烘焙厂位置的有力候选店？

纽约市是否有足够的医院应对机动车碰撞？该市哪些地区的覆盖率存在差距？

你还可以想象波士顿市的犯罪情况，检查加纳的卫生设施，探索欧洲顶尖大学，跟踪美国有毒化学品的排放。

在第一个教程中，我们将快速介绍完成本微型课程所需的先决条件。如果你想更深入地复习，我们推荐熊猫微课程。

我们还将开始可视化我们的第一个地理空间数据集！

读取数据

第一步是读入一些地理空间数据！为此，我们将使用GeoPandas库。

import geopandas as gpd

有很多不同的地理空间文件格式，如shapefile、GeoJSON、KML和GPKG。我们不会在本微课程中讨论它们的差异，但重要的是要提到：

shapefile是您将遇到的最常见的文件类型，并且

所有这些文件类型都可以通过gpd.read_file（）函数快速加载。

下一个代码单元加载一个形状文件，其中包含有关纽约州环境保护部管理下的森林、荒野地区和其他土地的信息。

# 读入数据
full_data = gpd.read_file("../input/geospatial-learn-course-data/DEC_lands/DEC_lands/DEC_lands.shp")

# 查看数据的前五行
full_data.head()

正如您在“CLASS”列中所看到的，前五行中的每一行对应于不同的林。

在本教程的其余部分，考虑一个场景，在这里你想使用这些数据来计划周末露营旅行。您决定创建自己的地图，而不是依靠在线众包评论。这样，你就可以根据自己的具体兴趣来调整行程。

先决条件

为了查看数据的前五行，我们使用head（）方法。您可能还记得，这也是我们用来预览熊猫数据帧的内容。事实上，您可以在DataFrame中使用的每个命令都可以处理数据！

这是因为数据已加载到（GeoPandas）GeoDataFrame对象中，该对象具有（Pandas）数据帧的所有功能。

type(full_data)

例如，如果我们不打算使用所有列，我们可以选择其中的一个子集(要查看选择数据的其他方法，请查看Pandas micro课程中的本教程。）

data = full_data.loc[:, ["CLASS", "COUNTY", "geometry"]].copy()

我们使用value_counts（）方法查看不同土地类型的列表，以及它们在数据集中出现的次数(要查看此教程（以及相关方法），请查看Pandas micro课程中的本教程。）

# 每种类型有多少块地？
data.CLASS.value_counts()

您还可以使用loc（和iloc）和isin来选择数据的子集(要回顾这一点，请查看Pandas micro课程中的本教程。）

# 选择属于“野生森林”或“荒野”类别的土地
wild_lands = data.loc[data.CLASS.isin(['WILD FOREST', 'WILDERNESS'])].copy()
wild_lands.head()

如果您不熟悉上面的命令，建议您将此页面标记为书签以供参考，以便根据需要查找这些命令(或者，您也可以学习Pandas micro课程。）在创建地图之前，我们将在整个micro课程中使用这些命令来理解和过滤数据。

创建您的第一张地图！

我们可以使用plot（）方法快速可视化数据。

wild_lands.plot()

每个GeoDataFrame都包含一个特殊的“几何体”列。它包含调用plot（）方法时显示的所有几何对象。

# 查看“几何图形”列中的前五个条目
wild_lands.geometry.head()

虽然此列可以包含各种不同的数据类型，但每个条目通常是点、线字符串或多边形。

数据集中的“几何体”列包含2983个不同的多边形对象，每个对象对应于上图中的不同形状。

在下面的代码单元中，我们再创建三个地理数据框，其中包含营地位置（点）、步道（线串）和县边界（多边形）。

# 纽约州的营地（点）
POI_data = gpd.read_file("../input/geospatial-learn-course-data/DEC_pointsinterest/DEC_pointsinterest/Decptsofinterest.shp")
campsites = POI_data.loc[POI_data.ASSET=='PRIMITIVE CAMPSITE'].copy()

# 纽约州的步道（线条）
roads_trails = gpd.read_file("../input/geospatial-learn-course-data/DEC_roadstrails/DEC_roadstrails/Decroadstrails.shp")
trails = roads_trails.loc[roads_trails.ASSET=='FOOT TRAIL'].copy()

# 纽约州的县边界（多边形）
counties = gpd.read_file("../input/geospatial-learn-course-data/NY_county_boundaries/NY_county_boundaries/NY_county_boundaries.shp")

接下来，我们从所有四个GeodataFrame创建一个地图。

plot（）方法将几个可用于自定义外观的参数作为（可选）输入。最重要的是，为ax设置值可确保所有信息都绘制在同一地图上。

# 定义具有县边界的底图
ax = counties.plot(figsize=(10,10), color='none', edgecolor='gainsboro', zorder=3)

# 将荒地、营地和步道添加到基础地图
wild_lands.plot(color='lightgreen', ax=ax)
campsites.plot(color='maroon', markersize=2, ax=ax)
trails.plot(color='black', markersize=1, ax=ax)

看起来该州的东北部将是露营旅行的绝佳选择！

2、坐标参考系统

介绍

您在本课程中创建的地图以二维形式描绘地球表面。但是，正如你所知，世界实际上是一个三维的地球。因此，我们必须使用一种称为贴图投影的方法将其渲染为平面。

地图投影不可能100%准确。每个投影都会以某种方式扭曲地球表面，同时保留一些有用的特性。例如，

等面积投影（如“兰伯特圆柱等面积”或“非洲阿尔伯斯等面积圆锥”）保留面积。例如，如果你想计算一个国家或城市的面积，这是一个很好的选择。

等距投影（如“方位等距投影”）保留距离。这将是计算飞行距离的一个好选择。

我们使用坐标参考系（CRS）来显示投影点如何对应于地球上的实际位置。在本教程中，您将了解有关坐标参照系的更多信息，以及如何在GeoPandas中使用坐标参照系。

import geopandas as gpd
import pandas as pd

设置CRS

当我们从shapefile创建GeoDataFrame时，CRS已经为我们导入。

# 加载包含加纳区域的GeoDataFrame
regions = gpd.read_file("../input/geospatial-learn-course-data/ghana/ghana/Regions/Map_of_Regions_in_Ghana.shp")
print(regions.crs)

你如何解释？

欧洲石油勘探集团（EPSG）规范引用了坐标参考系。

此GeoDataFrame使用EPSG 32630，通常称为“墨卡托”投影。此投影保留角度（使其对海上导航有用），并略微扭曲区域。

但是，当从CSV文件创建GeoDataFrame时，我们必须设置CRS。EPSG 4326对应于纬度和经度的坐标。

# 与加纳的医疗机构建立一个数据框架
facilities_df = pd.read_csv("../input/geospatial-learn-course-data/ghana/ghana/health_facilities.csv")

# 将数据帧转换为地理数据帧
facilities = gpd.GeoDataFrame(facilities_df, geometry=gpd.points_from_xy(facilities_df.Longitude, facilities_df.Latitude))

# 将坐标参考系（CRS）设置为EPSG 4326
facilities.crs = {'init': 'epsg:4326'}

# 查看GeoDataFrame的前五行
facilities.head()

在上面的代码单元中，要从CSV文件创建GeoDataFrame，我们需要同时使用Pandas和GeoPandas：

我们首先创建一个数据框，其中包含具有纬度和经度坐标的列。

要将其转换为GeoDataFrame，我们使用gpd.GeoDataFrame（）。

函数的作用是：通过纬度和经度列创建点对象。

重新投影

重新投影是指更改CRS的过程。这是在GeoPandas中使用to_crs（）方法完成的。

绘制多个地理数据框时，重要的是它们都使用相同的CRS。在下面的代码单元中，我们在绘制之前更改设施地理数据框的CRS，以匹配区域的CRS。

# 创建地图
ax = regions.plot(figsize=(8,8), color='whitesmoke', linestyle=':', edgecolor='black')
facilities.to_crs(epsg=32630).plot(markersize=1, ax=ax)

to_crs（）方法仅修改“geometry”列：所有其他列保持原样。

# “纬度”和“经度”列保持不变
facilities.to_crs(epsg=32630).head()

如果EPSG代码在GeoPandas中不可用，我们可以使用CRS的“proj4字符串”更改CRS。例如，要转换为纬度/经度坐标的proj4字符串如下所示：

# 将CRS更改为EPSG 4326
regions.to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs").head()

几何对象的属性

正如您在第一个教程中所了解的，对于任意的GeoDataFrame，“geometry”列中的类型取决于我们试图显示的内容：例如，我们可以使用：

地震震中的一个点，

街道的线条，或

显示国家边界的多边形。

所有三种类型的几何对象都具有内置属性，可用于快速分析数据集。例如，可以分别从x和y属性获取点的x坐标和y坐标。

# 获取每个点的x坐标
facilities.geometry.x.head()

并且，您可以从“长度”属性中获取线条字符串的长度。

或者，可以从“面积”属性获取多边形的面积。

# 计算GeoDataFrame中每个多边形的面积（以平方米为单位）
regions.loc[:, "AREA"] = regions.geometry.area / 10**6

print("Area of Ghana: {} square kilometers".format(regions.AREA.sum()))
print("CRS:", regions.crs)
regions.head()

在上面的代码单元中，由于regions GeoDataFrame的CRS设置为EPSG 32630（“墨卡托”投影），因此面积计算的精度略低于使用“非洲阿尔伯斯等面积圆锥”等面积投影时的精度。

但加纳的面积约为239585平方公里，与正确答案相差不远。

3、交互式地图

介绍

在本教程中，您将学习如何使用folium软件包创建交互式地图。在此过程中，您将应用新技能可视化波士顿犯罪数据。

import pandas as pd
import geopandas as gpd
import math

import folium
from folium import Choropleth, Circle, Marker
from folium.plugins import HeatMap, MarkerCluster

您的第一张交互式地图

我们首先用folium.map（）创建一个相对简单的映射。

# 创建地图
m_1 = folium.Map(location=[42.32,-71.0589], tiles='openstreetmap', zoom_start=10)

# 显示地图
m_1

几个参数可自定义映射的外观：

位置设置地图的初始中心。我们使用波士顿市的纬度（42.32°N）和经度（-71.0589°E）。

平铺更改地图的样式；在本例中，我们选择OpenStreetMap样式。如果你好奇，你可以找到这里列出的其他选项。

zoom_start设置贴图的初始缩放级别，其中值越高，缩放距离越近。

现在就花点时间通过放大和缩小，或者通过向不同方向拖动地图来探索。

数据

现在，我们将向地图添加一些犯罪数据！

我们将不关注数据加载步骤。相反，您可以想象您已经在数据帧中拥有了数据。数据的前五行如下所示。

# 加载数据
crimes = pd.read_csv("../input/geospatial-learn-course-data/crimes-in-boston/crimes-in-boston/crime.csv", encoding='latin-1')

# 删除缺少位置的行
crimes.dropna(subset=['Lat', 'Long', 'DISTRICT'], inplace=True)

# 2018年重点关注重大犯罪
crimes = crimes[crimes.OFFENSE_CODE_GROUP.isin([
    'Larceny', 'Auto Theft', 'Robbery', 'Larceny From Motor Vehicle', 'Residential Burglary',
    'Simple Assault', 'Harassment', 'Ballistics', 'Aggravated Assault', 'Other Burglary', 
    'Arson', 'Commercial Burglary', 'HOME INVASION', 'Homicide', 'Criminal Harassment', 
    'Manslaughter'])]
crimes = crimes[crimes.YEAR>=2018]

# 打印表格的前五行
crimes.head()

标绘点

为了减少我们需要在地图上显示的数据量，我们（暂时）将注意力集中在日间抢劫上。

daytime_robberies = crimes[((crimes.OFFENSE_CODE_GROUP == 'Robbery') & \
                            (crimes.HOUR.isin(range(9,18))))]

folium.Marker
We add markers to the map with folium.Marker(). Each marker below corresponds to a different robbery.

# 创建地图
m_2 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)

# 在地图上添加点
for idx, row in daytime_robberies.iterrows():
    Marker([row['Lat'], row['Long']]).add_to(m_2)

# 显示地图
m_2

folium.plugins.MarkerCluster

如果我们有很多标记要添加，folium.plugins.MarkerCluster（）可以帮助分离地图。每个标记都将添加到MarkerCluster对象。

# 创建地图
m_3 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)

# 在地图上添加点
mc = MarkerCluster()
for idx, row in daytime_robberies.iterrows():
    if not math.isnan(row['Long']) and not math.isnan(row['Lat']):
        mc.add_child(Marker([row['Lat'], row['Long']]))
m_3.add_child(mc)

# 显示地图
m_3

气泡图

气泡贴图使用圆而不是标记。通过改变每个圆的大小和颜色，我们还可以显示位置和其他两个变量之间的关系。

我们通过使用folium.Circle（）迭代添加圆来创建气泡图。在下面的代码单元格中，9-12小时发生的抢劫用绿色表示，而13-17小时发生的抢劫用红色表示。

# 创建底图
m_4 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)

def color_producer(val):
    if val <= 12:
        return 'forestgreen'
    else:
        return 'darkred'

# 将气泡贴图添加到基础贴图
for i in range(0,len(daytime_robberies)):
    Circle(
        location=[daytime_robberies.iloc[i]['Lat'], daytime_robberies.iloc[i]['Long']],
        radius=20,
        color=color_producer(daytime_robberies.iloc[i]['HOUR'])).add_to(m_4)

# 显示地图
m_4

请注意，folium.Circle（）接受几个参数：

位置是包含圆心的列表，以纬度和经度表示。

半径设置圆的半径。

请注意，在传统的气泡贴图中，允许每个圆的半径变化。我们可以通过定义一个类似于color_producer（）函数的函数来实现这一点，该函数用于改变每个圆的颜色。

颜色设置每个圆的颜色。

color_producer（）函数用于可视化时间对抢劫地点的影响。

热图

为了创建热图，我们使用folium.plugins.heatmap（）。这显示了城市不同区域的犯罪密度，红色区域的犯罪事件相对较多。

正如我们对大城市的预期，大多数犯罪发生在市中心附近。

# 创建底图
m_5 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=12)

# 将热图添加到基础地图
HeatMap(data=crimes[['Lat', 'Long']], radius=10).add_to(m_5)

# 显示地图
m_5

正如您在上面的代码单元中所看到的，folium.plugins.HeatMap（）接受几个参数：

数据是一个包含我们想要绘制的位置的数据框。

半径控制热图的平滑度。值越高，热图看起来越平滑（即间隙越小）。

乔洛佩斯地图

为了了解不同警区的犯罪情况，我们将创建一个choropleth地图。

作为第一步，我们创建一个GeoDataFrame，其中每个地区都分配了不同的行，“geometry”列包含地理边界。

# 波士顿警区地理边界的地理数据框架
districts_full = gpd.read_file('../input/geospatial-learn-course-data/Police_Districts/Police_Districts/Police_Districts.shp')
districts = districts_full[["DISTRICT", "geometry"]].set_index("DISTRICT")
districts.head()

我们还制作了一个名为plot_dict的熊猫系列，显示每个地区的犯罪数量。

# 每个警区的罪案数目
plot_dict = crimes.DISTRICT.value_counts()
plot_dict.head()

plot_dict与District具有相同的索引，这一点非常重要-这是代码知道如何用适当的颜色匹配地理边界的方式。

使用folium.Choropleth（）类，我们可以创建一个Choropleth映射。如果下面的地图无法为您呈现，请尝试在其他web浏览器中查看该页面。

# 创建底图
m_6 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=12)

# 将choropleth贴图添加到基础贴图
Choropleth(geo_data=districts.__geo_interface__, 
           data=plot_dict, 
           key_on="feature.id", 
           fill_color='YlGnBu', 
           legend_name='Major criminal incidents (Jan-Aug 2018)'
          ).add_to(m_6)

# 显示地图
m_6
Make this No

请注意folium.Choropleth（）接受几个参数：

geo_数据是一个GeoJSON特性集合，包含每个地理区域的边界。

在上面的代码中，我们将districts GeoDataFrame转换为带有_geo_interface__属性的GeoJSON FeatureCollection。

数据是一个系列，包含用于对每个地理区域进行颜色编码的值。

“打开”键将始终设置为feature.id。

这是指用于geo_数据的GeoDataFrame和数据中提供的Pandas系列具有相同的索引。为了理解细节，我们必须更仔细地查看GeoJSON功能集合的结构（其中与“features”键对应的值是一个列表，其中每个条目都是一个包含“id”键的字典）。

填充颜色设置颜色比例。

图例\名称在地图右上角标记图例。

4、操作地理空间数据

介绍

在本教程中，您将了解地理空间数据的两种常见操作：地理编码和表联接。

!pip install geopy==1.22.0
import pandas as pd
import geopandas as gpd
import numpy as np
import folium
from folium import Marker

地理编码

地理编码是将地名或地址转换为地图上某个位置的过程。例如，如果您曾经使用谷歌地图、必应地图或百度地图根据地标描述查找地理位置，那么您使用的是地理编码器！

我们将使用geopandas来完成所有的地理编码。

from geopandas.tools import geocode

要使用地理编码器，我们只需提供：

作为Python字符串的名称或地址，以及

供应商的名称；为了避免提供API密钥，我们将使用OpenStreetMap命名地理编码器。

如果地理编码成功，它将返回一个包含两列的GeoDataFrame：

“几何图形”列包含（纬度、经度）位置，以及

“地址”列包含完整地址。

result = geocode("The Great Pyramid of Giza", provider="nominatim")
result

“geometry”列中的条目是一个点对象，我们可以分别从y和x属性获得纬度和经度。

point = result.geometry.iloc[0]
print("Latitude:", point.y)
print("Longitude:", point.x)

通常情况下，我们需要对许多不同的地址进行地理编码。例如，假设我们想要获得欧洲100所顶尖大学的位置。

universities = pd.read_csv("../input/geospatial-learn-course-data/top_universities.csv")
universities.head()

然后，我们可以使用lambda函数将地理编码器应用于数据帧中的每一行(我们使用try/except语句来解释地理编码不成功的情况。）

def my_geocoder(row):
    try:
        point = geocode(row, provider='nominatim').geometry.iloc[0]
        return pd.Series({'Latitude': point.y, 'Longitude': point.x, 'geometry': point})
    except:
        return None

universities[['Latitude', 'Longitude', 'geometry']] = universities.apply(lambda x: my_geocoder(x['Name']), axis=1)

print("{}% of addresses were geocoded!".format(
    (1 - sum(np.isnan(universities["Latitude"])) / len(universities)) * 100))

# 删除未成功地理编码的大学
universities = universities.loc[~np.isnan(universities["Latitude"])]
universities = gpd.GeoDataFrame(universities, geometry=universities.geometry)
universities.crs = {'init': 'epsg:4326'}
universities.head()

接下来，我们可视化地理编码器返回的所有位置。请注意，其中一些位置肯定不准确，因为它们不在欧洲！

# 创建地图
m = folium.Map(location=[54, 15], tiles='openstreetmap', zoom_start=2)

# 在地图上添加点
for idx, row in universities.iterrows():
    Marker([row['Latitude'], row['Longitude']], popup=row['Name']).add_to(m)

# 显示地图
m

表联接

现在，我们将切换主题，思考如何组合来自不同来源的数据。

属性连接

您已经知道如何使用pd.DataFrame.join（）将来自多个数据帧的信息与共享索引组合在一起。我们将这种连接数据的方式（通过简化索引中的匹配值）称为属性连接。

当使用GeoDataFrame执行属性联接时，最好使用gpd.GeoDataFrame.merge（）。为了说明这一点，我们将使用包含欧洲每个国家边界的GeoDataFrame europe_Bounders。此GeoDataFrame的前五行打印在下面。

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
europe = world.loc[world.continent == 'Europe'].reset_index(drop=True)

europe_stats = europe[["name", "pop_est", "gdp_md_est"]]
europe_boundaries = europe[["name", "geometry"]]

europe_boundaries.head()

我们将加入一个数据框架europe_stats，其中包含每个国家的估计人口和国内生产总值（GDP）。

europe_stats.head()

我们在下面的代码单元中进行属性连接。on参数设置为列名，用于将europe_边界中的行与europe_统计中的行进行匹配。

# 使用属性联接合并有关欧洲国家的数据
europe = europe_boundaries.merge(europe_stats, on="name")
europe.head()

空间连接

另一种类型的连接是空间连接。通过空间连接，我们根据“几何体”列中对象之间的空间关系来组合GeoDataFrames。例如，我们已经有了一个包含欧洲大学地理编码地址的地理数据框架。

然后，我们可以使用一个空间连接来匹配每个大学与其对应的国家。我们使用gpd.sjoin（）来实现这一点。

# 使用空间连接将大学与欧洲国家匹配
european_universities = gpd.sjoin(universities, europe)

# 调查结果
print("We located {} universities.".format(len(universities)))
print("Only {} of the universities were located in Europe (in {} different countries).".format(
    len(european_universities), len(european_universities.name.unique())))

european_universities.head()

上面的空间连接查看两个GeodataFrame中的“geometry”列。如果大学地理数据框中的点对象与欧洲数据框中的多边形对象相交，则相应的行将合并并添加为欧洲大学数据框的单行。否则，没有匹配大学的国家（和没有配套国家的大学）将被忽略。

通过how和op参数，可以为不同类型的联接自定义gpd.sjoin（）方法。例如，您可以通过设置how='left'（或how='right'）来等效于SQL左（或右）联接。我们不会在本微型课程中详细介绍，但您可以在文档中了解更多信息。

5、邻近分析

介绍

在本教程中，您将探讨几种邻近性分析技术。特别是，您将学习如何执行以下操作：

测量地图上各点之间的距离，以及

选择特征某个半径内的所有点。

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()

测量距离

要测量来自两个不同地理坐标系的点之间的距离，我们首先必须确保它们使用相同的坐标参考系（CRS）。谢天谢地，这里就是这种情况，两者都使用EPSG2272。

print(stations.crs)
print(releases.crs)

我们还检查CRS以查看它使用的单位（米、英尺或其他）。在本例中，EPSG 2272具有英尺单位(如果您愿意，可以在此处查看。）

在GeoPandas中计算距离相对简单。下面的代码单元格计算最近发布中相对较新的发布事件与stations GeoDataFrame中的每个桩号之间的距离（英尺）。

# 特别选择一个发布事件
recent_release = releases.iloc[360]

# 测量从释放装置到每个工位的距离
distances = stations.geometry.distance(recent_release.geometry)
distances

使用计算出的距离，我们可以获得到每个站点的平均距离等统计信息。

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

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

print('Closest monitoring station ({} feet):'.format(distances.min()))
print(stations.iloc[distances.idxmin()][["ADDRESS", "LATITUDE", "LONGITUDE"]])

创建缓冲区

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

下面的代码单元将创建一个GeoSeries两英里缓冲区，其中包含12个不同的多边形对象。每个多边形是围绕不同空气监测站2英里（或2*5280英尺）的缓冲区。

two_mile_buffer = stations.geometry.buffer(2*5280)
two_mile_buffer.head()

我们使用folium.GeoJson（）在地图上绘制每个多边形。请注意，由于folium需要纬度和经度坐标，因此在绘图之前，我们必须将CRS转换为EPSG 4326。

# 创建包含发布事件和监视站的地图
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)
    
# 在地图上绘制每个多边形
GeoJson(two_mile_buffer.to_crs(epsg=4326)).add_to(m)

# 显示地图
m

现在，为了测试有毒物质是否在任何监测站的2英里范围内释放，我们可以对每个多边形运行12个不同的测试（分别检查它是否包含该点）。

但更有效的方法是首先将所有多边形折叠成一个多多边形对象。我们使用一元联合属性来实现这一点。

# 将多边形组转换为单个多多边形
my_union = two_mile_buffer.geometry.unary_union
print('Type:', type(my_union))

# 显示多多边形对象
my_union

我们使用contains（）方法检查multipolygon是否包含一个点。我们将使用本教程前面的发布事件，我们知道它距离最近的监测站大约3781英尺。

#最近的车站离这里不到两英里
my_union.contains(releases.iloc[360].geometry)

但并非所有的释放都发生在距空气监测站两英里的范围内！

# 最近的车站在两英里以外
my_union.contains(releases.iloc[358].geometry)

False