基于Python的地理空间分析（七）：ArcGIS API for Python

发表于 分类于 阅读次数： Disqus：

引言

本文将介绍ArcGIS的Python API，以及ArcGIS API的安装、使用和分析等。ArcGIS API for Python是一个用于处理地图和地理空间数据的Python库。我们可以用conda在本地安装API，然后与Esri的云GIS进行交互（无论是ArcGIS Online(SaaS)还是Portal for ArcGIS(为组织部署本地私有云)）。ArcGIS API能够配合Jupter使用，为使用Python脚本进行基于网络的地图分析提供了现代化的解决方案。

ArcGIS API for Python和ArcGIS Online介绍

Esri是一家以ArcGIS平台而闻名的地理空间软件公司，它将Python集成到他们的ArcGIS桌面软件及其后续ArcGIS Pro中。 Esri开发的第一个Python包是ArcPy，它是Python模块的集合，提供所有现有的以及扩展的ArcMap和ArcGIS Pro功能。 Python现在可以用作脚本和编程语言来自动执行涉及与图形用户界面（GUI）的大量交互的重复性任务。使用ArcPy，可以通过Python脚本，附加组件或工具箱执行这些任务。

Python成功地嵌入ArcGIS软件当中，而GIS本身正在向云平台发展——不仅是地理空间数据，还有软件本身。 Esri为组织提供了使用公共、私有或混合云服务通过各种云环境产品实现此目的的可能性。我们将使用ArcGIS Online，即软件即服务（Software as a Service, SaaS）产品，允许用户创建，存储和管理地图，应用程序和数据。在过去几年中，ArcGIS Online已成为Esri ArcGIS系统的关键组件和组成部分。其用户可以通过可用于网络，智能手机和平板电脑的现成可用工具在组织或全世界共享地图。

A Pythohnic Web API

为了能够让用户能够用Python与他们的GIS数据、服务等进行交互，Esri开发了一个全新的Python Web API，即ArcGIS API for Python。ArcGIS API for Python 由一整套子模块、工具和协议构成，用于创建软件和应用程序，建立在ArcGIS REporesentational State Transfer (REST) API和ArcGIS API for JavaScript的基础上，这一API同样用于Python API 2D和3D网络地图的显示。

GIS用户可以免费下载ArcGIS API for Python，用于管理他们的GIS云，无论是ARCGIS Online、ArcGIS Portal还是ArcGIS Enterprise（以前的ArcGIS Server系列）。目前该API要求Python 3.5或更高版本，可以与ArcPy一起使用，也可以在没有ArcPy或者任何桌面GIS产品的情况下工作。

自从2006年首次发布以来，API还在不断地更新完善之中。目前已具有数据处理和地图设计的多种功能，可以用作GIS可视化和分析、空间数据和内容管理以及组织管理等。

ArcGIS API for Python 模块介绍

与其他Python库一样，API具有模块、类、函数和类型，可用于管理和使用ArcGIS平台信息模型的元素。由于API适用于需要自己独特工具的不同用户组，因此API已组织为15个不同的模块（2019年2月）。我们没有必要全部介绍，本文主要用到的如下：

  1. GIS模块。这是最重要的模块，是ArcGIS Online或ArcGIS中托管的GIS的门户。 GIS模块允许您管理GIS中的用户、组和内容。在这里，GIS指的是用于创建，可视化和共享地图，场景，应用，图层，分析和数据的协作环境。
  2. 要素（features)模块。该模块代表API的矢量数据部分。矢量数据通过此模块表示为要素数据，要素图层或要素图层集合。单个数据元素由要素对象表示，而诸如FeatureSet，FeatureLayer和FeatureCollection的类表示要素数据的不同分组。
  3. 栅格（raster）模块。此模块包含用于处理栅格数据和图像层的类和栅格分析功能。而要素模块表示API的矢量数据组件，而栅格模块是栅格数据组件。该模块使用Imagerylayer类显示来自图像服务的数据，并提供用于实时图像处理的光栅功能。可以使用地图窗口小部件可视化图像图层。
  4. 地理处理（geoprocessing）模块。此模块是导入具有地理处理功能的工具箱所必需的，这些功能不属于API，但可通过ArcGIS Online获得。这些地理处理工具箱作为本机Python模块导入，以便您可以调用导入模块中可用的函数来调用这些工具。 API本身还包括丰富的地理处理工具集合，这些工具可通过空间数据类型定义的其他模块获得。地理处理工具是一种从输入数据集开始对GIS数据执行操作的功能。然后，对该数据集执行操作，最后将操作结果作为输出数据集返回。
  5. 小部件（widgets）模块：它提供了可视化GIS数据和分析的组件，并包含MapView Jupyter Notebook小部件。接下来我们将使用此小部件来可视化地图和图层。这不是唯一的可视化模块——单独的mapping模块提供不同的映射层和2D / 3D映射和可视化组件。

由此所见，API为不同的任务和用户提供了广泛的模块，包括发布映射数据，执行地理空间分析和数据操作。所有模块都使用Python作为脚本语言来管理GIS数据和功能。现在让我们开始使用API并探讨一些基本功能，然后再继续进行更高级的任务。

guide_api_overview_modules

安装与使用测试

ArcGIS API for Python 的安装与其他包类似，可以通过conda完成，但是其名字为arcgis，channel为esri，在命令行中运行。

1
conda install -c esri arcgis

从1.5.0版本开始，ArcGIS API开始支持Jupyterlab，为了在Jupyterlab中正常运行，还需要安装扩展：

1
2
jupyter labextension install @jupyter-widgets/jupyterlab-manager
jupyter labextension install [email protected]

相关命令参见 ：

https://developers.arcgis.com/python/guide/using-the-jupyter-lab-environment/ ，arcgis-map-ipywidget的版本号可以在 https://www.npmjs.com/package/arcgis-map-ipywidget 。

我们可以直接进行测试：

1
2
3
4
5
6
from arcgis.gis import GIS
from ipywidgets import Layout
mygis = GIS()
m = mygis.map()
m.layout.height = '600px'
m

output_20

1
# help(m.layout)
1
2
3
4
5
6
from ipywidgets import IntSlider
from ipywidgets.embed import embed_minimal_html
from ipywidgets.embed import embed_snippet

embed_minimal_html('map_widgets.html', views=[m], title='ArcGIS Map')
# embed_snippet(m, requirejs=True)

注意：在Python 3.7环境下可能出现： bad escape \u at position 0错误，原因是正则表达式的raw字符串处理，可以将embed.py文件中的return script_escape_re.sub(r'\u003c\1', s)替换为

1
return script_escape_re.sub('\u003c\\1', s)

具体参见：

https://stackoverflow.com/questions/54330673/how-to-fix-error-bad-escape-u-at-position-0

Troubleshooting

如果无法在本地安装和使用API，还可以尝试在云中运行的API的沙盒版本 https://notebooks.esri.com/。 单击此URL，将使用Jupyter Notebook打开浏览器窗口，就可以在其中创建自己的笔记本，运行代码示例以及使用API的所有功能。

有关显示所有模块的在线API参考，带有描述的类和示例，请参阅 http://esri.github.io/arcgis-python-api/apidoc/html/index.html 。

有关API更新，发行说明等信息，请参阅 https://developers.arcgis.com/python/guide/release-notes/ 。

有关API的所有信息的主页面在 https://developers.arcgis.com/python/ ， 它是一个很好的资源，包含大量文档，用户指南和API参考。

连接Esri帐户

正如我们之前所说，API可以管理可以位于云环境中的Web GIS并与之交互。既然我们已经在我们的机器上安装了API，现在是时候讨论如何将它与Esri用户帐户结合使用。 为了能够使用API并与此Web或云GIS进行交互，我们需要注册的sri用户帐户，通过用户名和密码来与此Web GIS建立连接，确保服务器和客户端之间的安全连接以及对正确内容的访问。

不同类型的Esri帐户

Esri为不同类型的Esri帐户提供了不同级别的功能权限：

  1. 匿名用户帐户（Anonymous User Account），无需传递任何用户信息即可访问ArcGIS Online。这是测试某些基本功能的快速解决方案，但不提供个性化帐户附带的任何高级功能；
  2. ArcGIS Online组织帐户（或Portal for ArcGIS帐户），需要（付费）订阅ArcGIS Online或Portal for ArcGIS。此选项为您提供尽可能多的功能，本文不予介绍；
  3. ArcGIS Enterprise试用帐户,免费试用21天，为您提供创建地图和发布内容所需的服务功能，到期之后必须转移到付费帐户才能继续；
  4. 免费的Esri开发者帐户，此帐户是ArcGIS Developer计划的一部分，该计划为开发和测试个人应用程序以及使用ArcGIS Online等提供50个服务；
  5. 创建公共ArcGIS帐户并使用Web浏览器登录ArcGIS Online。使用这些登录详细信息，您现在可以使用API连接到ArcGIS Online，但功能有限，此选项已在API的1.3版本中添加，我们也基本不用。

综上所述，我们就用免费的开发者帐户好了。

ArcGIS API for Python使用实例

使用API和地图widget

在前面的测试中我们已经导入了API，并且用widget展示出一副地图（可知无需使用帐户）。我们可以传入地名作为参数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
bj_map = mygis.map('Beijing')

# adjust the zoom level
bj_map.zoom = 6

# query the type of current basemap
print(bj_map.basemap)

# query the available basemaps list
print(bj_map.basemaps)

# quer the coordinate reference system
print(bj_map.extent)

# use satellite as basemap
bj_map.basemap = 'satellite'
bj_map.layout.height = '600px'
bj_map
# Requir JS at the first time for future use in html file
# embed_snippet(bj_map, requirejs=True)
default
['dark-gray', 'dark-gray-vector', 'gray', 'gray-vector', 'hybrid', 'national-geographic', 'oceans', 'osm', 'satellite', 'streets', 'streets-navigation-vector', 'streets-night-vector', 'streets-relief-vector', 'streets-vector', 'terrain', 'topo', 'topo-vector']
{'xmin': 116.12223000000003, 'ymin': 39.63250000000003, 'xmax': 116.67223000000004, 'ymax': 40.182500000000026}

output_24

1
bj_map.extent
{'xmin': 116.12223000000003,
 'ymin': 39.63250000000003,
 'xmax': 116.67223000000004,
 'ymax': 40.182500000000026}

实际上，与许多地图一样，我们呈现出的图像也是网络地图，而且是以JSON格式存储、传输和表示的，因此也可以通过在后台查看分析其结构。

看一下bj_map.extent的内容，可以发现地图的坐标系信息存储在spatialReference对象中，还可以看到其空间参考坐标系的代码为'latestWkid': 3857, 'wkid': 10210，查询 http://developers.arcgis.com/ 可知 两者都代表的是网络墨卡托投影（Web Mercator projection），这是目前网络地图应用的事实标准，已经被许多在线地图服务商所使用。

上面的一些查询操作是在匿名情况下完成的，接下来我们将使用开发者帐户进行更多操作。

地理空间内容的搜索、展示和描述

创建Esri帐户的方法非常简单，无需赘述。接下来我们将创建一个地图widget，并且保存到ArcGIS Online的个人文件夹，并且对一个要素图层文件进行一系列的操作。首先要用自己的用户名和密码连接：

1
devgis = GIS(" https://username.maps.arcgis.com", "username", "password")
1
2
3
4
5
6
7
8
# open a map of Beijing
map = devgis.map("Solano County, USA")
map.zoom = 13

# Search trails with content.search() function
# outside_org specifies whether to search for data outside of our own organization
search_result = devgis.content.search(query='san francisco trail', item_type='Feature Layer', outside_org=True)
search_result
[<Item title:"National Park Service - Park Unit Boundaries" type:Feature Layer Collection owner:NationalParkService>,
 <Item title:"National Park Service - Park Unit Centroids" type:Feature Layer Collection owner:NationalParkService>,
 <Item title:"JUBA_RecreationRoute" type:Feature Layer Collection owner:[email protected]_nps>,
 <Item title:"JUBA_HistoricTrail_ln" type:Feature Layer Collection owner:[email protected]_nps>,
 <Item title:"National Park Service - Park Unit Boundaries - copy" type:Feature Layer Collection owner:rfischer_ClimateWest>,
 <Item title:"BikeTrails" type:Feature Layer Collection owner:jlgoicochea>,
 <Item title:"UGA_IFP_Editable Route_Athens to Portland" type:Feature Layer Collection owner:juliaec_USG6>,
 <Item title:"JUBA_ExpeditionCampsites" type:Feature Layer Collection owner:[email protected]_nps>,
 <Item title:"JUBA_HistoricCorridor_py" type:Feature Layer Collection owner:[email protected]_nps>,
 <Item title:"Creek to Peaks Trail" type:Feature Layer Collection owner:sstasio>]

1
2
3
4
from IPython.display import display

for item in search_result:
    display(item)

我们对其中的自行车道（bike trails）最感兴趣，那么需要把它添加到我们的地图当中：

1
2
3
4
# reference the data source
bike_trails_item = search_result[5]
bike_trails_item.layers

[<FeatureLayer url:" https://services2.arcgis.com/SCn6czzcqKAFwdGU/arcgis/rest/services/BikeTrails/FeatureServer/0">,
 <FeatureLayer url:" https://services2.arcgis.com/SCn6czzcqKAFwdGU/arcgis/rest/services/BikeTrails/FeatureServer/1">]

1
2
3
# get the names of layers
for lyr in bike_trails_item.layers:
    print(lyr.properties.name)
BikeTrails
Parks

1
2
3
4
5
6
7
bike_trails_layer = bike_trails_item.layers[0]
for field in bike_trails_layer.properties['fields']:
    print(field['name'])
# store the attribute table in a pandas dataframe
bike_df = bike_trails_layer.query().sdf
print(bike_df.shape)
bike_df.head()
OBJECTID
SHAPESTLen
STRNAME
City
NAME_PCASE
FIPS
Shape_Leng
LASTUPDATE
LASTEDITOR
ParkName
Shape_Le_2
LastDateEd
GlobalID
Shape__Length
(406, 15)

从AcrGIS API for Python 1.5开始，启用 the Spatially Enabled DataFrame (SEDF)，关于SEDF的信息参见： https://developers.arcgis.com/python/guide/introduction-to-the-spatially-enabled-dataframe/#Example:-Feature-Layer-Query-Results-to-a-Spatially-Enabled-DataFrame

由于AcrGIS API for Python还在活跃开发中，因此其数据结构、接口等可能经常变化，使用时请注意版本并以官方技术文档为准。

接下来，将自行车道的要素添加到我们的地图当中。

1
2
3
4
map.add_layer(bike_trails_layer)
map.layout.height = '600px'
map
# JS not required from the second time for future use in html file

output_12

作为ArcGIS Python API新增加的一项功能，这对于在ArcGIS Online上定制自己的网络地图非常有用。除此以外，我们还可以用mapping模块中的WebMap创建并存储地图，具体如下：

1
2
3
4
5
6
7
8
9
10
11
from arcgis.mapping import WebMap
# log into your GIS
devgis = GIS(" https://username.maps.arcgis.com", "username", "password")
wm = WebMap()
web_map_properties = {'title':'Bike Trails ',
                      'snippet':'This map service shows bike trails in Solano County',
                      'tags':'ArcGIS Python API'}
wm.add_layer(bike_trails_layer, options={'title':'Bike trails layer'})
web_map_item = wm.save(item_properties=web_map_properties)
web_map_item

Bike Trails
This map service shows bike trails in Solano CountyWeb Map by [email protected]
Last Modified: 二月 25, 2019
0 comments, 0 views

从ArcGIS Online分享地图：

## API处理栅格数据

前面提到ArcGIS API for Python具有raster子模块，这次我们只需要使用匿名帐户处理Landsat 8卫星图像的栅格数据。

1
2
3
4
5
6
7
8
9
10
11
12
# import package and modules
import arcgis
from arcgis.gis import GIS
from IPython.display import display
gis = GIS()

# search for content used for the raster module
items = gis.content.search('Landsat 8 Views', item_type='Imagery Layer', max_items=10)

# dispalyt the items
for item in items:
    display(item)
Pansharpened Landsat
Landsat 8 OLI, 15m multitemporal PanSharpened Natural Color images. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
0 comments, 38,321 views
Panchromatic Landsat
Landsat 8 OLI, 15m multitemporal Panchromatic images. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
0 comments, 26,288 views
Landsat 8 Pansharpened
Landsat 8 OLI, 15m multitemporal PanSharpened Natural Color images. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
1 comments, 118,721 views
Multispectral Landsat
Landsat 8 OLI, 30m multispectral and multitemporal 8-band imagery, with on-the-fly renderings and indices. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
0 comments, 116,508 views
Landsat 8 Panchromatic
Landsat 8 OLI, 15m multitemporal Panchromatic images. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
0 comments, 45,193 views
MDA NaturalVue Satellite Imagery
This image service presents NaturalVue 15m satellite imagery of the world created by MDA Information Systems Inc.Imagery Layer by esri
Last Modified: 八月 25, 2017
0 comments, 35,011 views
Landsat 8 Views
Landsat 8 OLI, 30m multispectral and multitemporal 8-band imagery, with on-the-fly renderings and indices. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: 五月 04, 2018
0 comments, 155,547 views
Landsat Arctic Views - Geology with DRA
Landsat 8 OLI, 30m Multispectral 8 imagery in the Arctic, with visual renderings and indices. Based on the Landsat on AWS collections, this layer is in beta.Imagery Layer by esride_content
Last Modified: 十月 06, 2016
0 comments, 84 views
Color Infrared Multi-spectral Image Service
Multi-spectral Landsat using the bands near-infrared, red, green (5, 4, 3) with fixed stretch applied on apparent reflectance. Healthy vegetation is bright red while stressed vegetation is dull red.Imagery Layer by boulw2_bpresponse
Last Modified: 二月 01, 2018
0 comments, 90 views

我们选择第7个试试：

1
2
3
4
l8_view = items[6]
print(l8_view.layers)
l8_layer = l8_view.layers[0]
l8_layer

output_47_1

查看栅格数据的相关属性:

1
2
3
4
5
from IPython.display import HTML
import pandas as pd
l8_layer.properties
display(HTML(l8_layer.properties.description))
pd.DataFrame(l8_layer.key_properties()['BandProperties'])

Multispectral Landsat 8 OLI Image Service covering the landmass of the World. This service includes 8-band multispectral scenes, at 30 meter resolution. It can be used for mapping and change detection of agriculture, soils, vegetation health, water-land features and boundary studies. Using on-the-fly processing, the raw DN values are transformed to scaled (0 - 10000) apparent reflectance values and then different service based renderings for band combinations and indices are applied. The service is updated on a daily basis to include the latest best scenes from the USGS.

BandName WavelengthMax WavelengthMin
0 CoastalAerosol 450.0 430.0
1 Blue 510.0 450.0
2 Green 590.0 530.0
3 Red 670.0 640.0
4 NearInfrared 880.0 850.0
5 ShortWaveInfrared_1 1650.0 1570.0
6 ShortWaveInfrared_2 2290.0 2110.0
7 Cirrus 1380.0 1360.0
8 NaN NaN NaN

raster模块的方法信息

1
2
for fn in l8_layer.properties.rasterFunctionInfos:
    print(fn['name'])
Agriculture with DRA
Bathymetric with DRA
Color Infrared with DRA
Natural Color with DRA
Short-wave Infrared with DRA
Geology with DRA
Agriculture
Bathymetric
Color Infrared
Geology
Natural Color
Short-wave Infrared
NDVI Colorized
Normalized Difference Moisture Index Colorized
NDVI Raw
NBR Raw
None

1
2
3
# add our satellite image to the map widget,
bj_map.add_layer(l8_layer)
bj_map

output_18

1
2
3
4
5
6
7
8
9
10
11
12
# to apply the raster functions
from arcgis.raster.functions import apply

# create a natural color image with dynamic range adjustment（DRA）, using bands 4, 3, and 2:
natural_color = apply(l8_layer, 'Natural Color with DRA')
bj_map.add_layer(natural_color)

# visualize the agricultural map, using the bands 6, 5 and 2, referring to
# shortwave IR-1, near-IR and blue, respectively.
agric = apply(l8_layer, 'Agriculture')
bj_map.add_layer(agric)
bj_map

output_26

小结

本文介绍了全新的ArcGIS API for Python， 包括API的不同模块、地图窗口小部件的使用、ArcGIS Online以及矢量和栅格数据的分析。我们可以利用存储在基于云的ArcGIS Online系统中的数据，用ArcGIS API for Python for Python来执行基本的地理空间分析和创建ArcGIS Online Web地图等任务。

接下来将介绍用于基于AWS云基础架构的Elasticsearch和MapD GPU数据库。我们将通过Python与二者进行交互，完成地理空间搜索、地理定位数据处理等任务。