使用 Python 和 Esri 的免费资源进行地理空间可视化

区域增长的一个指标是新建筑发放的许可证数量。探索和分析许可证活动有助于区域规划者确保开发符合该区域的长期目标。蒙哥马利县是马里兰州的一个郊区县，靠近华盛顿特区，近期经历了快速增长。县规划者希望观察空间和时间上的增长趋势，找出某些区域增长速度快于其他区域的原因，并向公众传达有关该县增长的关键信息。

在本例中，我们将使用 Esri（全球领先的位置智能平台）提供的免费资源，探索马里兰州蒙哥马利县的建筑许可数据。

我们将使用 ArcGIS API for Python，这是一个强大且易于使用的 Python 库，用于在空间背景下管理、可视化和分析数据，同时还可以访问 ArcGIS Online 的内容——一个基于 SaaS 的空间分析和制图解决方案。

您可以使用以下方式将其安装到您的 Conda 环境中：

conda install -c esri arcgis

准备好 ArcGIS API for Python 后，打开一个 Jupyter notebook，并导入 pandas 和 ArcGIS API for Python。

from arcgis.gis import GIS
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from arcgis.features import GeoAccessor, GeoSeriesAccessor

我们将使用 GIS 模块连接到 ArcGIS Online 的免费内容。

gis = GIS()

搜索“Commercial Permits since 2010”（2010 年以来的商业许可证）图层。我们将指定所有者姓名以获得更具体的结果。

data = gis.content.search(
    'title: Commercial Permits since 2010 owner: rpeake_LearnGIS',
    'Feature layer')
data[0]

获取结果中的第一项。

permits = data[0]

地理空间数据通常由多个图层组成，使用 ArcGIS API for Python，我们可以使用此要素图层集合的 layers 属性，它为我们提供了一个集合内的 FeatureLayer 对象列表。许可证图层是此项目中的第一个图层。

permit_layer = permits.layers[0]

接下来，我们将在马里兰州蒙哥马利县的地图上可视化此图层。为此，我们首先创建一个空白地图，然后将图层添加到地图中。

permit_map = gis.map('Montgomery County, Maryland’)
permit_map

permit_map.add_layer(permit_layer)

该地图是完全交互式的，可以平移、缩放，并且地图上的点可以被点击以显示有关每个点的属性弹出窗口。但是，与其单独点击每个点，不如探索空间数据框中的数据。

数据探索

空间数据框的外观和感觉与 pandas 数据框相似，这使得它们非常容易上手。

首先，让我们将空间图层转换为空间数据框。

sdf = pd.DataFrame.spatial.from_layer(permit_layer)
sdf.head()

for col in sdf.columns:
    print(col)

我们也可以使用图表来探索数据。

plt.figure(figsize=(9,15))
ax=sns.countplot(y='Use_Code', data=sdf)

plt.axis('equal') 
permits_by_status.plot(kind='pie', legend=False, label='Permits by Status');

上面的饼图显示了四种许可证状态，每种状态的大小由许可证数量决定。绝大多数许可证已“Issued”（已发布）或“Finaled”（已完成）。已完成的许可证是已发布并且已执行必要检查的许可证。

在地图上可视化许可证属性的空间分布很有帮助。您将更改地图，使每个许可证的符号代表其状态。

permits_by_status_map = gis.map('Montgomery County, Maryland')
permits_by_status_map

sdf.spatial.plot(kind='map', map_widget=permits_by_status_map,
        renderer_type='u', # specify the unique value renderer using its notation 'u'
        col='Status')  # column to get unique values from

可视化时间和空间趋势

数据显示了许可证的打开和发布位置，但这些许可证对何时何地正在该县发生增长有何说明？数据还包含时间属性字段，例如 Added_Date（添加日期），表示许可证首次添加到系统的时间。该字段有多个值，按年、月甚至小时细分数据。

我们将拆分 Added_date 以获得 year（年份）、month（月份）、week_of_day（星期几）。

sdf['datetime'] = pd.to_datetime(sdf['Added_Date'], unit='ms')
sdf['year'], sdf['month'], sdf['day_of_week'] = sdf.datetime.dt.year, sdf.datetime.dt.month, sdf.datetime.dt.dayofweek

sns.countplot(x="year", data=sdf);

该图表显示了自 2010 年以来每年发放的许可证数量。（2017 年的许可证数量显著减少，因为该数据集仅涵盖 2017 年的一部分。）您可以通过每个条形的大小来直观地比较许可证数量。虽然每年都会有些波动，但大多数年份的许可证活动相似。

同样，您也可以按月份可视化。

sns.countplot(x="month", data=sdf);

我们可以结合这些信息来可视化许可证发放随时间变化的趋势和模式。

ddf['num'] = 1
ddf['num'].resample('M').sum().plot();

当我们检查此图表时，我们看到 2011 年中期出现了巨大的许可证活动高峰。

是什么导致了这个高峰？是整体许可证活动的增加，还是主要是某种特定类型许可证的增加？为了更好地理解这个异常值，我们将根据 Use_Code（使用代码）绘制许可证数量，以找出是否存在任何特定类型的许可证导致了高峰。

fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(1, 1, 1)

ax.plot(ddf['num'].resample('M').sum(), 'k', label='Total permits')
for use_code in filtered_permits:
    x = ddf[ddf.Use_Code == use_code]['num'].resample('M').sum()
    ax.plot(x, label=use_code)
ax.legend();

根据图例，2011 年的许可证活动高峰是由于多户住宅许可证数量急剧增加。这很可能意味着 2011 年出现了大规模的居民区增长。

我们要回答的最后一个问题是，2011 年的许可证高峰在哪里最密集？为此，我们将使用我们之前提取的时间信息和一个热点地图渲染器。

hotspot_2011_map = gis.map('Germantown, Montgomery County, Maryland')
hotspot_2011_map

sdf[sdf.year==2011].spatial.plot(kind='map', map_widget=hotspot_2011_map,
        renderer_type='h',
        col='Status')  # column to get unique values from

有关 ArcGIS 中可用的空间分析和数据科学工具的更多信息，请访问 Esri 的空间分析和数据科学页面以获取详细信息，或参加我们即将举行的为期 6 周的免费空间数据科学大规模开放在线课程（MOOC）。

此外，如果您有兴趣使用 ArcGIS API for Python 启动您自己的空间分析工作流程，请通过ArcGIS for Developers 基础计划免费开始。