ArcPy-和-ArcGIS-指南-全-

ArcPy 和 ArcGIS 指南(全)

原文:zh.annas-archive.org/md5/b868c8d13da5b4bad684622354841e5d

译者:飞龙

协议:CC BY-NC-SA 4.0

前言

ArcGIS,来自行业领导者 ESRI 的 GIS 软件,允许分析和展示地理空间数据。

将 Python 集成到 ArcGIS 中使得 ArcPy 模块成为 GIS 学生和专业人士的重要工具。ArcPy 模块提供了一种强大的方式来提高进行地理空间分析时的生产力。从基本的 Python 脚本编写到高级的 ArcPy 方法和属性,ArcPy 和其他 Python 模块将提高任何 GIS 工作流程的速度和可重复性。

本书将指导您从基本的 Python 脚本编写到高级脚本工具。它专注于地理空间分析脚本编写,并涉及自动化地图输出。到本书结束时,您将能够创建可重用模块,将可重复分析作为脚本工具添加到 ArcToolbox 中,并自动导出地图。通过将 GIS 的耗时性质从几天减少到几小时,一位 GIS 专业人士可以变得像整个团队一样强大。

本书涵盖内容

第一章, ArcGIS 的 Python 简介,快速介绍了 Python 的基础知识,包括编程语言的其他用途。它涵盖了 Python 数据类型和本书中使用的的重要模块。

第二章, 配置 Python 环境,介绍了 Python 的工作原理:其文件夹结构、可执行文件和模块。它还解释了将模块导入脚本、内置模块,并涵盖了集成开发环境(IDE),这些是强大的编程辅助工具。

第三章, 编写第一个 Python 脚本,展示了如何使用 ArcGIS ModelBuilder 来建模第一个分析,然后将其导出为 Python 脚本。还介绍了字符串操作以及如何在 Python 中使用文件路径。

第四章, 复杂的 ArcPy 脚本和泛化函数,探讨了如何执行分析并生成使用 ModelBuilder 无法实现的结果。通过使用函数或可重用代码块,避免了重复代码。

第五章, ArcPy 游标 – 搜索、插入和更新,涵盖了 ArcPy 数据访问游标及其如何用于在要素类和表中搜索、更新或插入记录。它解释了使用游标迭代的怪癖,以及如何仅选择或更新感兴趣的记录。

第六章, 使用 ArcPy 几何对象,探讨了 ArcPy 几何对象及其如何与游标结合以执行空间分析。它展示了如何使用数据游标和 Arcpy 几何类型进行缓冲、裁剪、重投影等操作,而不使用 ArcToolbox。

第七章, 创建脚本工具,解释了如何将脚本转换为在 ArcToolbox 中出现的工具,并且这些工具是动态的。它解释了工具和脚本如何通信,以及如何设置 ArcTool 对话框以正确传递参数给脚本。

第八章, ArcPy 映射简介,探讨了强大的 Arcpy.Mapping 模块以及如何修复损坏的图层链接、开启和关闭图层、动态调整标题和文本。它展示了如何根据地理空间分析创建动态地图输出。

第九章, 更多 ArcPy 映射技术,介绍了图层对象及其方法和属性。它演示了如何控制数据帧的地图比例和范围,并涵盖了自动地图导出。

第十章, 高级几何对象方法,扩展了 ArcPy 几何对象的方法和属性。它还解释了如何创建模块以保存代码,以便在后续脚本中重复使用,并演示了如何创建包含地理空间分析结果的 Excel 电子表格。

第十一章, 使用 ArcPy 进行网络分析和空间分析,介绍了使用 ArcGIS for Desktop 网络分析器和空间分析扩展进行高级地理空间分析的基本方法。

第十二章, 起点结束,涵盖了需要理解的其他重要主题,以便全面掌握 ArcPy。这些主题包括环境设置、XY 值和 Z 和 M 分辨率、空间参考系统(投影)、Describe 函数等。

您需要这本书的内容

您需要 ArcGIS 10.1/10.2/10.3 的专有版或免费版。为了支持您的环境,您需要 2GB RAM、32 位或 64 位机器硬件配置,以及 Windows 7/8。需要 Python 2.7 来进行编程,并且它将与 ArcGIS 一起安装。

这本书面向的对象

这本书旨在为需要了解如何使用 ArcPy 来减少重复性任务和提高分析速度的 GIS 学生和专业人员提供帮助。对于希望了解如何使用行业标准 ArcGIS 软件及其 ArcPy 模块来自动化地理空间分析的 Python 程序员来说,这也是一本有价值的书。

惯例

在这本书中,您将找到多种文本样式,用于区分不同类型的信息。以下是一些这些样式的示例及其含义的解释。

文本中的代码单词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 用户名如下所示:"这两个数据片段,BusStopID 和averatePop变量随后被添加到一个列表中。"

代码块应如下设置:

with arcpy.da.SearchCursor(Intersect71Census, ["STOPID","POP10"]) as cursor:
    for row in cursor:
        busStopID = row[0]
        pop10 = row[1]
        if busStopID not in dataDictionary.keys():
            dataDictionary[busStopID] = [pop10]
        else:
            dataDictionary[busStopID].append(pop10)

任何命令行输入或输出都应如下所示:

>>> aString = "This is a string"
>>> bString = " and this is another string"
>>> aString + bString 

新术语重要词汇将以粗体显示。您在屏幕上看到的单词,例如在菜单或对话框中,将以如下方式显示:“通过单击它,然后单击编辑按钮来选择它。”

注意

警告或重要提示将以如下框中的形式出现。

小贴士

小技巧和窍门将以如下形式出现。

读者反馈

我们始终欢迎读者的反馈。请告诉我们您对这本书的看法——您喜欢什么或可能不喜欢什么。读者反馈对我们开发您真正能从中受益的标题非常重要。

要向我们发送一般性反馈,只需发送电子邮件至 <feedback@packtpub.com>,并在邮件主题中提及书名。

如果你在某个领域有专业知识,并且对撰写或参与一本书籍感兴趣,请参阅我们的作者指南,网址为 www.packtpub.com/authors

客户支持

现在您已经是 Packt 书籍的骄傲拥有者,我们有一些事情可以帮助您从您的购买中获得最大收益。

下载示例代码

您可以从您在 www.packtpub.com 的账户中下载您购买的所有 Packt 书籍的示例代码文件。如果您在其他地方购买了这本书,您可以访问 www.packtpub.com/support 并注册,以便将文件直接通过电子邮件发送给您。

下载本书的颜色图像

我们还为您提供了一个包含本书中使用的截图/图表的颜色图像的 PDF 文件。这些颜色图像将帮助您更好地理解输出的变化。您可以从 www.packtpub.com/sites/default/files/downloads/8662OS_ColorImages.pdf 下载此文件。

错误清单

尽管我们已经尽最大努力确保内容的准确性,但错误仍然可能发生。如果您在我们的书中发现错误——可能是文本或代码中的错误——如果您能向我们报告这一点,我们将不胜感激。这样做可以帮助其他读者避免挫败感,并帮助我们改进本书的后续版本。如果您发现任何错误清单,请通过访问 www.packtpub.com/submit-errata,选择您的书籍,点击错误提交表单链接,并输入您的错误详情。一旦您的错误清单得到验证,您的提交将被接受,错误清单将被上传到我们的网站或添加到该标题的错误部分下的现有错误清单中。

要查看之前提交的错误清单,请访问 www.packtpub.com/books/content/support,并在搜索字段中输入书名。所需信息将出现在错误清单部分下。

侵权

互联网上对版权材料的盗版是一个跨所有媒体的持续问题。在 Packt,我们非常重视我们版权和许可证的保护。如果您在互联网上发现我们作品的任何非法副本,无论形式如何,请立即提供位置地址或网站名称,以便我们可以寻求补救措施。

请通过链接<版权@packtpub.com>与我们联系,提供疑似盗版材料的链接。

我们感谢您在保护我们作者以及为我们提供有价值内容的能力方面的帮助。

问题

如果您在本书的任何方面遇到问题,请通过链接<问题@packtpub.com>与我们联系,我们将尽力解决。

第一章:ArcGIS 的 Python 简介

在本章中,我们将讨论 Python 作为编程语言的发展,从 20 世纪 80 年代末的起点到现在的状态。我们将讨论推动其发展的设计哲学,并简要介绍本书中将使用的重要模块,特别是关注内置在 Python 标准库中的模块。这种对语言及其功能的概述将有助于解释是什么让 Python 成为 ArcGIS 自动化的优秀语言。

本章将涵盖:

  • Python 快速概述:它是什么,做什么,谁创建了它,现在在哪里

  • ArcPy模块和其他重要模块

  • Python 作为通用编程语言

Python 概述

Python,由吉多·范罗瑟姆于 1989 年创建,以他最喜欢的喜剧团体*蒙提·派森命名。当时他的工作组有一个将程序命名为电视节目的传统,他想要一个不拘一格且与众不同的名字——ABC、Pascal、Ada、Eiffel、FORTRAN 和其他。因此,他选择了 Python,感觉这个名字有点儿前卫和吸引人。当然,比起 Python 所基于的语言 C 来说,这个名字听起来更有趣。

现在,Python 是一种主要的编程语言。它被用于网站开发、数据库管理,甚至用于编程机器人。对 GIS 分析师来说最重要的是,Python 可以用来控制 ArcGIS 工具和地图文档,使用优秀的ArcPy模块以有组织和快速的方式产生地理空间数据和地图。

ArcPy 与 ArcGIS 桌面版和 ArcGIS 服务器版一起安装。自 ArcGIS 10.0 以来,ArcPy 一直是官方的 ArcGIS 脚本语言,并在功能和实现方面稳步改进。本书将针对 ArcGIS 桌面版 10.1 及以后的版本,并演示如何利用 Python 及其强大的编程库(或模块)来构建复杂的地理空间分析。

Python 作为编程语言

在过去的 40 年里,编程语言已经从汇编语言和机器代码发展到了更接近英语的高级抽象语言。Python 编程语言的设计是为了克服 1980 年代程序员所抱怨的许多问题:开发速度慢、语法过于复杂、可读性差。范罗瑟姆希望开发一种能够实现快速代码开发和测试、具有简单或至少可读的语法,并且用更少的代码行、更短的时间产生结果的编程语言。Python 的第一个版本(0.9.0)于 1991 年发布,并且从一开始就可以免费获取;在“开源”这个术语被发明之前,Python 就是开源的。

解释型语言

Python 是一种解释型语言。它是用 C 语言编写的,C 是一种编译型语言,代码在执行之前从 Python 解释成 C。实际上,这意味着代码在转换和编译后立即执行。虽然代码解释可能会对基于 Python 的程序的执行速度产生影响,但 Python 允许的快速开发时间使得这个缺点很容易被忽略。在解释型环境中测试代码片段的速度要快得多,非常适合创建自动化基本、可重复的计算任务的脚本。Python 脚本具有 .py 扩展名。一旦代码被解释,就会生成第二个 Python 脚本(具有 .pyc 扩展名),以保存编译后的代码。当原始 .py 脚本中的更改时,.pyc 脚本将自动重新编译。

标准库(内置库)

Python,一旦安装,就有一组基本的功能,被称为标准库。这些工具允许 Python 执行字符串操作、数学计算、HTTP 调用和 URL 解析,以及许多其他功能。一些工具库,Python 程序员称之为模块,是内置的,并且一旦启动 Python 就可用,而其他模块必须使用 import 关键字显式调用,以使它们的函数和类可用。其他模块是由第三方开发的,可以根据需要下载并安装到 Python 安装中。

许多新程序员会 wonder 如果 Python 是一种真正的编程语言,这是一个有争议的问题。答案是肯定的;Python 可以用来创建完整的程序,构建网站,运行计算机网络,以及更多。内置模块和附加模块使 Python 非常强大,它可以(并且已经被)用于计算机的几乎所有部分——操作系统、数据库、网络服务器、桌面应用程序等等。对于这些工具的开发,Python 并非总是最佳选择,但这并没有阻止程序员们尝试,甚至取得成功。

粘合语言

当 Python 被用作粘合语言时,它表现得最好。这个术语描述了使用 Python 来控制其他程序,向它们发送输入并收集输出,然后将这些输出发送到另一个程序或写入磁盘的过程。一个 ArcGIS 的例子是使用 Python 从网站上下载压缩的 shapefile,解压文件,使用 ArcToolbox 处理文件,并将结果编译成 Excel 电子表格。所有这些都可以使用免费提供的模块完成,这些模块要么包含在 Python 的标准库中,要么在安装 ArcGIS 时添加。

包装模块

ArcPy 模块是一个包装模块。在 Python 中,包装模块很常见,之所以这样命名是因为它们将 Python 包装到我们将需要的工具上。它们允许我们使用 Python 与用 C 或其他编程语言编写的其他程序进行接口,使用这些程序的 应用程序编程接口API)。例如,包装器使得从 Excel 电子表格中提取数据并将数据转换或加载到另一个程序(如 ArcGIS)成为可能。并非所有模块都是包装器;一些模块是用纯 Python 编写的,并使用 Python 语法执行分析和计算。无论如何,最终结果是计算机及其程序都可以通过 Python 进行操作和控制。

Python 的禅 是为了简单、易读和简化而创建的,与之前存在的其他语言相比。这种指导哲学由早期 Python 开发者 Tim Peters 组织成一首诗,称为 Python 的禅;它是一个隐藏特性(隐藏功能),包含在每个 Python 安装中,并在输入 import this 时在 Python 解释器中显示:

The Zen of Python, by Tim Peters:
Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

注意

前往 www.python.org/doc/humor/ 获取更多信息。

Python 的基础知识

Python 有许多语言要求和约定,允许控制模块和代码结构。以下是一些重要的基本概念,这些概念将在本书中以及编写用于地理空间分析的脚本时使用。

导入语句

导入语句用于通过调用其他模块来增强 Python 的功能,以便在脚本中使用。这些模块可以是标准 Python 库的一部分,例如 math 模块(用于进行高级数学计算)或,更重要的是,ArcPy,它将使我们能够与 ArcGIS 进行交互。

注意

导入 语句可以在使用模块之前的位置放置,但按照惯例,它们位于脚本的顶部。

创建 导入 语句有三种方式。第一种,也是最标准的方式,是如下导入整个模块:

import arcpy

  • 使用这种方法,我们甚至可以在同一行导入多个模块。以下导入三个模块:arcpyos(操作系统模块)和 sys(Python 系统模块):

    import arcpy, os, sys
    
    
  • 导入脚本的下一个方法是导入模块的特定部分,而不是导入整个模块,使用 from <module> import <submodule> 语法:

    from arcpy import mapping
    
    
  • 当只需要 ArcPy 的一部分代码时,使用这种方法;它实际上限制了模块在调用时使用的内存量。我们也可以以相同的方式导入模块的多个部分:

    from arcpy import mapping, da
    
    
  • 导入模块的第三种方式是使用 from <module> import <submodule> 语法,但可以通过使用星号来导入模块的所有部分:

    from arcpy import *
    
    

这种最后一种方法仍然在使用,但被劝阻,因为它可能产生不可预见的后果。例如,如果模块中的变量没有明确导入,它们可能与另一个模块中的另一个变量冲突。因此,最好避免这种第三种方法。然而,许多现有的脚本包含这种类型的导入语句,因此请注意这些后果。

变量

变量是所有编程语言的一部分。它们用于引用数据并在脚本中存储以供以后使用。关于变量命名的最佳方法有很多争议。Python 脚本编写尚未制定标准。以下是一些在命名变量时应该遵循的最佳实践。

  • 使它们具有描述性:不要仅仅给变量命名为 x;当脚本被审查时,这个变量将毫无用处,因为无法知道它被用于什么,或者为什么。它们应该更长而不是更短,并且应该暗示它们引用的数据或甚至它们引用的对象的数据类型:

    shapefilePath = 'C:/Data/shapefile.shp'
    
    
  • 使用驼峰命名法使变量可读:驼峰命名法是指变量以小写字母开头,中间有大小写字母,类似于骆驼的驼峰:

    camelCase = 'this is a string'
    
    
  • 在变量名中包含数据类型:如果变量包含字符串,可以称其为 variableString。这不是必需的,本书中也不会教条式地使用,但它可以帮助组织脚本,并且对阅读这些脚本的其他人很有帮助。Python 是 动态类型 而不是 静态类型。一种编程语言的区别意味着变量在使用之前不需要声明,这与 Visual Basic 或其他静态类型语言不同。这提高了编写脚本的效率,但在长脚本中可能会出现问题,因为变量的数据类型可能不明显。

注意

ArcGIS 在导出 Python 脚本时不会使用驼峰命名法,许多示例也不会包含它;尽管如此,在编写新脚本时仍然建议使用它。此外,变量不能以数字开头。

循环

编程语言内置了遍历或执行重复过程的能力,以对数据集进行转换或提取满足特定标准的数据。Python 的主要迭代工具被称为 for 循环。for 循环这个术语意味着一个操作将循环或迭代数据集中的项目以对每个项目执行操作。数据集必须是可迭代的才能在 for 循环中使用,这一点将在后面进一步讨论。

我们将在整本书中使用循环。以下是一个简单的例子,它使用 Python 解释器来获取字符串值并以大写格式打印它们,使用循环:

>>> newlist = [ 'a' , 'b' , 'c' , 'd' ]
>>> for item in newlist:
 print item.upper()

输出如下所示:

A
B
C
D

变量 item 是一个通用变量,它在 for 循环中每个对象被输入时分配给每个对象,而不是 Python 所要求的术语。它可以是 x 或 value。在循环中,第一个对象(a)被分配给通用变量 item,并对其应用了字符串函数以产生输出A。一旦执行了这个动作,下一个对象(b)就被分配给通用变量以产生输出。这个循环会重复进行所有数据集 newlist 的成员;一旦完成,变量 item 将仍然携带数据集的最后一个成员的值(在这种情况下是d)。

小贴士

下载示例代码

您可以从您在www.packtpub.com的账户下载您购买的所有 Packt 书籍的示例代码文件。如果您在其他地方购买了这本书,您可以访问www.packtpub.com/support并注册,以便将文件直接通过电子邮件发送给您。

If/Elif/Else 语句

条件语句,在 Python 中称为 if/else 语句,在编程语言中也是标准的。它们用于评估数据;当满足某些条件时,将执行一个动作(初始的 if 语句;如果满足另一个条件,将执行另一个动作;这是一个elif语句),如果数据不满足条件,将分配一个最终的动作来处理这些情况(else 语句)。这些与在 ArcToolbox 中与 Select 工具一起使用的 SQL 语句中的 where 条件类似。以下是如何使用 if/else 语句评估列表中的数据(一个将在后面进一步讨论的数据类型)并使用取模运算符(%)和 Python 的等于运算符(==)找到余数的示例:

>>> data = [1,2,4,5,6,7,10]
>>> for val in data:
 if val % 2 == 0:
 print val,"no remainder"
 elif val % 3 == 2:
 print val, "remainder of two"
 else:

print "final case"

输出如下所示:

final case
2 no remainder
4 no remainder
5 remainder of two
6 no remainder
final case
10 no remainder

While 语句

另一个重要的评估工具是while语句。它用于在条件为真时执行一个动作;当条件为假时,评估将停止。请注意,条件必须变为假,否则动作将始终执行,创建一个无限循环,直到外部关闭 Python 解释器才会停止。以下是一个使用 while 循环执行动作直到条件变为假的示例:

>>> x = 0
>>> while x < 5:
 print x

x+=1

输出如下所示:

0
1
2
3
4

注释

Python 中的注释用于在脚本中添加注释。它们由井号标记,当脚本运行时,Python 解释器会忽略它们。注释有助于解释代码块在执行时的作用,或者添加有用的注释,供脚本作者希望未来的脚本用户阅读:

# This is a comment

虽然在编程中有一个常见的真理,即良好的代码是注释良好的代码,但许多程序员会跳过这个有价值的步骤。此外,过多的注释可能会降低其有用性和脚本的可读性。如果变量描述得足够充分,且代码组织得很好,注释就变得不那么必要;尽可能地将代码写得详细且组织得越好,就越不需要花费太多时间在注释上。

数据类型

GIS 使用点、线、多边形、覆盖和栅格来存储数据。在执行分析时,这些 GIS 数据类型可以以不同的方式使用,并且具有不同的属性和特征。Python 与 GIS 类似,也有组织数据的数据类型。Python 中的主要数据类型包括字符串、整数、浮点数、列表、元组和字典。它们各自具有自己的属性和特征(或属性),并用于代码自动化的特定部分。此外,还有一些内置函数允许将数据类型从一种类型转换为另一种类型;例如,可以使用str()函数将整数 1 转换为字符串 1:

>>> variable = 1
>>> newvar = str(variable)
>>> newvar

输出如下所示:

1

字符串

字符串用于包含任何类型的字符。它们以引号开始和结束,可以使用单引号或双引号,但字符串必须以相同类型的引号开始和结束。在字符串内部,可以出现引号文本;它必须使用相反的引号以避免与字符串冲突。查看以下示例:

>>> quote = 'This string contains a quote: "Here is the quote" '

还有一种第三种字符串类型也被使用,即以三个单引号开始和结束的多行字符串:

>>> multiString = '''This string has
multiple lines and can go for
as long as I want it too'''

整数

整数是没有小数位的整数。在数学运算中使用整数有一些特殊的结果;如果使用整数进行除法,将返回整数结果。查看下面的代码片段以查看此例子的示例:

>>> 5 / 2

输出如下所示:

2

Python 不会返回准确的结果 2.5,而是返回向下取整的值,即任何整数除法计算中的最低整数。这显然可能是一个问题,并可能导致脚本中的小错误,这些错误可能具有重大后果。

小贴士

请在编写脚本时注意此问题,并使用浮点数来避免以下章节中描述的问题。

浮点数

浮点值或浮点数用于 Python 表示小数。在执行除法时推荐使用浮点数:

>>> 5.0 / 2

输出如下所示:

2.5

由于计算机以 2 为基数存储值,因此可能存在在基数为 10 的系统中通常表示的浮点值表示问题。阅读docs.python.org/2/tutorial/floatingpoint.html以进一步讨论此限制的影响。

列表

列表是有序的数据集,包含在方括号([])中。列表可以包含任何其他类型的数据,包括其他列表。单个列表中可以混合不同的数据类型。列表还具有一组方法,允许它们被扩展、反转、排序、求和,或提取最大或最小值,以及许多其他方法。列表中的数据项由逗号分隔。

列表成员是通过它们的索引或列表中的位置来引用的,索引始终从零开始。查看以下示例以更好地理解这一点:

>>> alist = ['a','b','c','d']
>>> alist[0]

输出如下所示:

'a'

此示例展示了如何从名为alist的列表中提取第一个值(索引 0)。一旦列表被填充,列表内的数据通过其索引引用,该索引传递到方括号中的列表。要获取列表中的第二个值(索引 1 的值),使用相同的方法:

>>> alist[1]

输出如下所示:

'b'

要合并两个列表,使用extend方法:

>>> blist = [2,5,6]
>>> alist.extend(blist)
>>> alist

输出如下所示:

['a', 'b', 'c', 'd', 2, 5, 6]

元组

元组与列表相关,并用圆括号(())表示。与列表不同,元组是不可变的——一旦创建,就不能调整或扩展。元组内的数据通过索引引用,与列表相同,索引从零开始:

>>> atuple = ('e','d','k')
>>> atuple[0]

输出如下所示:

'e'

字典

字典用花括号({})表示,并用于创建键:值对。这允许我们将值从键映射到值,因此值可以替换键,并且可以从值中提取数据用于处理。以下是一个简单的示例:

>>> adic = {'key':'value'}
>>> adic['key']

输出如下所示:

'value'

注意,与引用索引位置不同,例如列表或元组,值是通过键来引用的。此外,键可以是任何其他类型的数据,除了列表(因为列表是可变的)。

这在读取 shapefile 或要素类时非常有价值。使用ObjectID作为键,值将是一个与ObjectID关联的行属性列表。查看以下示例以更好地理解此行为:

>>> objectIDdic = { 1 : [ '100' , 'Main' , 'St' ] }
>>> objectIDdic[1]

输出如下所示:

['100', 'Main', 'St']

字典在读取要素类和通过仅调用感兴趣的行等操作轻松解析数据时非常有价值。它们非常适合对数据进行排序和重新排序,以便在脚本中稍后使用,因此请确保继续关注它们。

可迭代数据类型

列表、元组和字符串都是可迭代数据类型,可以在 for 循环中使用。当进入 for 循环时,除非另有说明,否则这些数据类型按顺序操作。对于列表和元组,这很容易理解,因为它们有明显的顺序:

>>> aList = [1,3,5,7]
>>> for value in aList:
 print value * 2

输出如下所示:

2
6
10
14

对于字符串,每个字符都会进行循环:

>>> aString = "esri"
>>> for value in aString:
 print value.upper()

输出如下所示:

E
S
R
I 

字典也是可迭代的,但具有特定的实现,只允许直接访问字典的键(然后可以用来访问值)。此外,键的返回顺序不是特定的:

>>> aDict = {"key1":"value1",
 "key2":"value2"}
>>> for value in aDict:
 print value, aDict[value]

输出如下所示:

key2 value2
key1 value1 

其他重要概念

使用 Python 进行编程需要介绍许多概念,这些概念要么是 Python 特有的但必需的,要么是常见的编程概念,在创建脚本时将反复使用。以下是一些必须涵盖的概念,以便在 Python 中流利地使用。

缩进

与大多数其他编程语言不同,Python 对代码行的缩进执行严格的规则。这一概念再次源于 Guido 希望产生干净、可读的代码的愿望。在创建函数或使用 for 循环、if/else 语句时,需要在后续的代码行上进行缩进。如果 for 循环包含在 if/else 语句中,将有两层缩进。其他语言的资深程序员曾抱怨 Python 缩进的严格性。新程序员通常发现这很有帮助,因为它使代码组织变得容易。请注意,许多新接触 Python 的程序员会在某个时候创建缩进错误,所以请确保注意缩进级别。

函数

函数用于将脚本中反复出现的代码,或跨脚本中的代码,转化为正式的工具。使用关键字def(代表定义函数),函数通过定义的输入和输出创建。在计算机科学中,函数的概念是它将数据从一个状态转换为另一个状态,而不影响脚本的其他部分。这对于自动化 GIS 分析非常有价值。

下面是一个函数的示例,该函数返回任何给定数字的平方:

def square(inVal):
 return inVal ** 2
>>> square(3)

输出如下所示:

9

虽然这当然与数学模块中内置的类似函数重复,但它展示了函数的基本原理。函数(通常)接受数据,按需转换它,然后使用 return 关键字返回数据的新状态。

关键字

Python 内置了许多关键字,在命名变量时应避免使用。这些包括max, min, sum, return, list, tuple, def, del, from, not, in, as, if, else, elif, or, while, and, with等,还有很多。使用这些关键字会导致错误。

命名空间

当函数内部(局部变量)的变量与函数外部(全局变量)的变量同名时,命名空间是一种逻辑上组织变量名的方法。函数内部(无论是脚本还是导入的模块中)包含的局部变量和全局变量可以同名,只要它们不共享命名空间。

当导入的模块中的变量意外地与脚本中的变量同名时,这个问题经常出现。Python 解释器将使用命名空间规则来决定哪个变量被调用,这可能导致不希望的结果。

基于零的索引

如前文所述,在描述列表和元组的章节中提到,Python 的索引和计数从零开始,而不是从一。这意味着一组数据中的第一个成员位于零位置,第二个成员位于第一个位置,依此类推,直到最后一个位置。

当脚本中有 for 循环迭代时,此规则也适用。当迭代开始时,正在迭代的第一个数据成员位于零位置。

此外,还可以从可迭代对象的最后一个成员开始计数时进行索引。在这种情况下,最后一个成员的索引是-1,倒数第二个是-2,依此类推,直到对象的第一个成员。

重要的 GIS 分析 Python 模块

模块,或可以通过脚本调用来增加其编程潜力的代码库,要么是内置在 Python 中的,要么是由第三方创建并在之后添加到 Python 中的。其中大部分是用 Python 编写的,但也有一些是用其他编程语言编写的,然后被封装在 Python 中以在 Python 脚本中使用。模块还用于使其他程序对 Python 可用,例如 Microsoft Word 中内置的工具。

ArcPy 模块

ArcPy 模块是一个包装模块,用于与 ArcGIS 工具交互,这些工具随后由 ArcGIS 以其内部代码格式执行,同时也是一个允许对地理空间分析和地图制作进行额外控制的代码库。ArcPy 用于控制 ArcToolbox 中的工具,但这些工具并没有用 Python 重写;相反,我们能够通过 ArcPy 使用 ArcGIS 工具。ArcPy 还使我们能够控制 ArcGIS 地图文档(MXDs)及其包含的对象:图例、标题、图像、图层以及地图视图本身。ArcPy 还具有 ArcToolbox 中不可用的工具。其中最强大的是数据游标,特别是新的数据分析游标,它为 GIS 数据创建了一个更 Pythonic 的接口。数据游标在第五章 ArcPy 游标:搜索、插入和更新 和第六章 使用 ArcPy 几何对象 中有详细的介绍,对于从数据源中提取数据行进行分析非常有用。

使用 ArcPy 控制地理空间分析的能力允许将 ArcGIS 工具集成到包含其他强大 Python 模块的工作流程中。Python 的粘合语言能力通过减少对地理空间数据进行特殊处理的需求,增加了 ArcGIS 的有用性。

操作系统(OS)模块

操作系统模块(OS 模块)是标准库的一部分,它允许 Python 访问操作系统功能。该模块的常见用途是使用os.path方法通过将路径分为目录路径(即文件夹)和基本路径(即文件)来控制文件路径。还有一个有用的方法,os.walk,它将遍历目录并返回文件夹和子文件夹中的所有文件。在执行 GIS 分析时,会不断访问 OS 模块。

Python 系统(SYS)模块

标准库中的 sys 模块指的是 Python 安装本身。它有许多方法可以获取已安装 Python 版本的信息,以及有关脚本和任何提供给脚本(或参数)的信息,使用sys.argv方法。sys.path方法非常有用,可以追加 Python 文件路径;实际上,这意味着包含脚本的文件夹可以被其他脚本引用,以便它们包含的函数可以被其他脚本导入。

XLRDXLWT模块

XLRDXLWT模块分别用于读取和写入 Excel 电子表格。这些模块在从旧版电子表格中提取数据并将其转换为 GIS 分析可用的数据,或在完成地理空间分析后写入分析结果时非常有用。它们不是 Python 标准库的一部分,但与 ArcGIS 10.2 和 Python 2.7 一起安装。

常用内置函数

在本书中,我们将使用许多内置函数。主要列如下:

  • str:字符串函数用于将任何其他类型的数据转换为字符串

  • int:整数函数用于将字符串或浮点数转换为整数。为了避免错误,传递给整数函数的任何字符串都必须是数字,例如 1。

  • float:浮点函数用于将字符串或整数转换为浮点数,类似于整数函数。

常用标准库模块

以下标准库模块必须导入:

  • datetime:datetime 模块用于获取有关日期和时间的信息,并将字符串日期转换为 Python 日期。

  • math:math 模块用于在必要时进行高级数学函数,例如获取π的值或计算数字的平方。

  • string:字符串模块用于字符串操作。

  • csv:CSV 模块用于创建和编辑逗号分隔值类型文件。

查阅docs.python.org/2/library以获取标准库中所有内置模块的完整列表。

摘要

在本章中,我们讨论了 Python 的禅意,并介绍了使用 Python 进行编程的基础知识。我们开始探索 ArcPy 以及它如何与其他 Python 模块集成以生成完整的工作流程。我们还讨论了 Python 标准库和 Python 的基本数据类型。

接下来,我们将讨论如何配置 Python 以与 ArcGIS 一起使用,并探讨如何使用集成开发环境(IDEs)编写脚本。

第二章:配置 Python 环境

在本章中,我们将配置 Python 和我们的计算机,以便它们协同工作以执行 Python 脚本。我们将配置路径变量和环境变量,以确保导入语句按预期工作,并且当点击脚本时能够运行。我们将讨论 Python 文件夹的结构,以及 ArcPy 模块在 ArcGIS 文件夹结构中的位置。我们还将讨论 集成开发环境IDEs),这些程序旨在协助代码的创建和执行,并比较和对比现有的 IDE,以确定每个 IDE 在编写 Python 脚本时可以提供哪些好处。

本章将涵盖:

  • Python 解释器的位置以及如何调用它来执行脚本

  • 调整计算机的环境变量以确保正确执行代码

  • 集成开发环境

  • Python 的文件夹结构,重点关注模块存储的位置

什么是 Python 脚本?

让我们从编写和执行 Python 脚本的基础开始。什么是 Python 脚本?它是一个包含一系列按顺序编写的、用正式化语言编写的命令的简单文本文件。该文本文件的扩展名为 .py,但除此之外,它与其他文本文件没有区别。它可以使用记事本或 Wordpad 等文本编辑器打开,但 Python 的魔力并不在于 Python 脚本中。没有 Python 解释器,Python 脚本无法运行,其中的命令也无法执行。

Python 如何执行脚本

理解 Python 如何解释脚本并执行其中的命令,与理解 Python 语言本身一样重要。通过花时间正确设置 Python,可以避免数小时的调试和错误检查。Python 的解释性意味着脚本必须首先转换为字节码,然后才能执行。我们将介绍 Python 实现自动化 GIS 分析目标所采取的步骤。

什么是 Python 解释器?

在 Windows 环境下,Python 解释器是一个编译成 Windows 可执行文件的程序,其扩展名为 .exe。Python 解释器 python.exe 是用 C 语言编写的,这是一种较老且广泛使用的编程语言,其语法较为复杂。

用 C 语言编写的程序,最初也是作为文本文件编写的,必须通过编译器转换成可执行文件,编译器是一种将文本命令转换成机器码以创建可执行程序的专用程序。这是一个缓慢的过程,可能会使 C 语言编写简单程序变得费力。好处是产生的程序是独立程序,可以在没有任何依赖的情况下运行。另一方面,Python 快速解释和执行 Python 命令,这使得它成为一款优秀的脚本语言,但脚本必须通过解释器运行,不能单独执行。

如其名称所暗示的,Python 解释器解释 Python 脚本中包含的命令。当运行或执行 Python 脚本时,首先检查其语法以确保它符合 Python 的规则(例如,遵循缩进规则,变量遵循命名约定)。然后,如果脚本有效,其中的命令会被转换成字节码,这是一种由字节码解释器执行的专用代码,字节码解释器是一个用 C 语言编写的虚拟机。字节码解释器进一步将字节码(包含在以.pyc扩展名结尾的文件中)转换成用于计算机的正确机器码,然后 CPU 执行脚本。这是一个复杂的过程,使得 Python 能够保持其简单的外观。

还有一些 Python 解释器的版本是用 Java 编写的(称为 Jython)和用.NET 编写的(称为 IronPython);这些变体用于在其他计算环境中编写 Python 脚本,本书将不会涉及这些内容。ArcGIS 安装程序包括 Python 的标准实现,也称为 CPython,以区分这些变体。

Python 解释器的位置在哪里?

在计算机文件夹结构中 Python 解释器的位置是一个需要掌握的重要细节。Python 通常直接从www.python.org下载并独立于 ArcGIS 安装。然而,每个 ArcGIS 版本都需要特定版本的 Python;考虑到这一要求,将 Python 包含在 ArcGIS 安装包中是有帮助的。对于本书,我们将使用 ArcGIS 10.2,这将需要 Python 2.7。

在 Windows 机器上,Python 文件夹结构直接放置在 C:驱动器上,除非它被明确加载到其他驱动器上。ArcGIS 10.2 的安装过程将在C:\Python27创建一个文件夹,该文件夹将包含一个名为ArcGIS10.2ArcGIS10.2x64的文件夹,具体取决于操作系统和已安装的 ArcGIS 版本。对于本书,我将使用 32 位版本的 ArcGIS,因此最终的文件夹路径将是C:\Python27\ArcGIS10.2

在此文件夹中包含多个子文件夹,以及 python.exe(Python 解释器)。还包括解释器的第二个版本,称为 pythonw.exepythonw.exe 将在没有终端窗口且程序反馈出现的情况下执行脚本。python.exepythonw.exe 都包含所有 Python 命令的完整副本,可以用来执行脚本。

应使用哪个 Python 解释器?

使用 Python 解释器直接执行脚本的通用规则是使用 pythonw.exe,因为没有终端窗口会出现。当需要测试代码片段或需要在终端窗口中查看输出时,通过双击可执行文件来启动 python.exe

python.exe 启动时,将出现一个 Python 解释器控制台:

应使用哪个 Python 解释器?

注意在版本信息说明下方出现的独特的三个箭头(>>>)。这是 Python 提示符,在这里输入代码逐行执行,而不是在一个完整的脚本中执行。这种直接访问解释器的方法对于测试代码片段和理解语法非常有用。自 ArcGIS 10 以来,这种解释器版本,即 Python 窗口,已被集成到 ArcMap 和 ArcCatalog 中。这将在后面的章节中进一步讨论。

计算机是如何知道解释器位置的?

要能够直接执行 Python 脚本(即通过双击脚本使其运行),计算机还需要知道解释器在其文件夹结构中的位置。要完成此操作,需要管理员账户访问权限和对 Windows 如何搜索程序的高级知识。我们必须在高级系统设置对话框中调整一个环境变量,以将解释器注册到系统路径中。

在 Windows 7 机器上,点击开始菜单,然后右键单击 计算机,从菜单中选择 属性。在 Windows 8 机器上,点击 Windows 资源管理器,然后右键单击 此电脑,并从菜单中选择 属性。这些命令是获取 控制面板的系统安全/系统 菜单的快捷方式。从左侧面板中选择 高级系统 设置。点击出现在 系统属性 菜单底部的 环境变量 按钮。在 环境变量 菜单的下半部分,滚动通过 系统变量 窗口,直到出现 Path 变量。通过单击它,然后单击 编辑 按钮来选择它。以下窗口将出现:

计算机是如何知道解释器位置的?

此变量有两个组成部分:变量名(路径)和变量值。变量值是一系列用分号分隔的文件夹路径。这是当 Windows 查找与文件扩展名关联的特定可执行文件时搜索的路径。在我们的情况下,我们将添加包含 Python 解释器的文件夹路径。将C:\Python27\ArcGIS10.2(或您机器上的等效路径)输入到变量值字段中,确保用分号将其与之前的值分开。点击确定退出编辑对话框,然后点击确定退出环境变量菜单,最后点击确定退出系统属性菜单。现在,计算机将知道 Python 解释器的位置,因为它将搜索 Path 变量中包含的所有文件夹以查找名为 Python 的可执行文件。为了测试路径调整是否正确,打开命令窗口(开始菜单/运行cmd),然后输入python。解释器应该直接在命令窗口中运行:

计算机如何知道解释器的位置?

如果出现带有版本信息和三重箭头的 Python 头文件,则路径调整已正确完成。

注意

如果没有管理员访问权限,有一个解决方案。在命令行窗口中,传递 Python 解释器的完整路径(C:\Python27\ArcGIS10.2\python.exe)以启动解释器。

点击时使 Python 脚本可执行

要使脚本在双击时运行(这也意味着它们可以在 ArcGIS 环境之外运行,从而节省大量内存开销),最后一步是将具有.py扩展名的文件与 Python 解释器关联起来。如果脚本尚未与解释器关联,它们将显示为未知类型的文件或文本文件。

要更改此设置,右键单击一个Python脚本。选择打开方式,然后选择选择默认程序。如果python.exepythonw.exe没有作为选项出现,导航到包含它们的文件夹(在这种情况下为C:\Python27\ArcGIS10.2)并选择python.exepythonw.exe。再次强调,两者之间的区别在于使用python.exe运行脚本时会出现终端窗口,该窗口将包含脚本的任何输出(但脚本完成后此窗口将消失)。我建议在执行脚本时使用pythonw.exe,在测试代码时使用python.exe

注意

Python 脚本也可以通过将扩展名调整为.pyw而不是.py来显式调用pythonw.exe

集成开发环境(IDE)

Python 解释器包含执行 Python 脚本或通过交互与 Python 解释器进行交互以测试 Python 代码所需的一切。然而,编写脚本需要文本编辑器。通常,Windows 机器上至少包含两个简单的文本编辑器(记事本和写字板),它们可以在紧急情况下编辑脚本,甚至编写整个脚本。不幸的是,它们非常简单,不允许用户执行使编写多个脚本或非常长的脚本更容易的功能。

为了弥合这一差距,一系列被称为集成开发环境(IDE)的程序被开发出来。所有编程语言都有 IDE,包括变量列表、代码辅助等功能,使它们非常适合编写编程脚本。我们将回顾其中的一些,以评估它们在编写 Python 脚本方面的有用性。以下讨论的三个都是免费的,并在不同的 Python 社区中建立了良好的声誉。

IDLE

当 Python 安装时,它包含一个 IDE。这个 IDE 被称为 IDLE,这个名字是对 IDE 和 Monty Python 的一位杰出成员埃里克·艾德尔(Eric Idle)名字的双关语。在 Windows 7 中,可以通过转到开始菜单并找到 程序 菜单中的 ArcGIS 文件夹来启动它。在 Python 文件夹中,IDLE 将是该文件夹中的一个选项。选择它以启动 IDLE

IDLE

IDLE 包含一个交互式解释器(即三重箭头)和运行完整 Python 脚本的能力。它还使用 Python 内置的 GUI 模块 Tkinter 编写,因此它具有使用执行同一语言的优势。

使用 IDLE 而不是 Python 控制台(python.exe)的另一个优点是,任何打印语句或其他脚本输出都会被定向到 IDLE 交互窗口,该窗口在执行脚本后不会消失。与内存使用相比,IDLE 也非常轻量。脚本通过 文件 菜单中的文件对话框打开,最近运行的脚本列在 文件 菜单的 最近文件 中。

IDLE 的缺点包括有限的代码辅助(或代码自动完成),这是一个有用的 IDE 工具,以及没有方法将脚本组织成逻辑项目。无法找到脚本中包含的所有变量,这是其他 IDE 的另一个有用功能。此外,最近文件菜单对列出的脚本数量有限制,这使得找到几个月未运行的脚本变得更加困难(相信我,这是常见的情况!)如果机器上无法安装其他程序,IDLE 是一个可用的 IDE。它也非常适用于快速测试代码片段。虽然它不是我的主要 IDE,但我发现自己几乎每天都在使用 IDLE。

PythonWin

PythonWin(即 Python for Windows)可在sourceforge.net/projects/pywin32/files/pywin32找到,并包括用于在 Windows 环境中使用 Python 的 IDE 和有用的模块。选择PythonWin的最新构建版本,然后根据已安装的 Python 版本选择正确的 32 位模块(对于我的机器,我选择了pywin32-218.win32-py2.7.exe,这是我的 32 位 Python 2.7 安装的正确版本)。运行可执行文件,如果已下载正确版本,安装 GUI 将识别系统注册表中的 Python 2.7 并自动安装。

PythonWin 包括一个交互式窗口,用户可以直接与 Python 解释器交互。脚本也可以在 PythonWin 中打开,它包括 Windows 菜单中的一组平铺命令,允许用户组织所有打开的脚本和交互式窗口的显示。

PythonWin 相较于 IDLE 的另一个优点是能够在同一脚本窗口中显示脚本的各个部分。如果一个脚本变得过长,在编辑时上下滚动脚本可能会很麻烦。PythonWin 允许用户从脚本顶部下拉以创建第二个脚本窗口,该窗口可以专注于脚本的另一部分。此外,在左侧还可以打开另一个窗口,该窗口将列出 Python 类和变量,这使得导航到脚本特定部分变得更加容易。

PythonWin

PythonWin 的交互式窗口中集成了一个小巧但实用的功能,即能够搜索之前输入的代码语句。在三个箭头提示符下,按住Ctrl键,使用上下箭头键在行之间导航,以找到感兴趣的行。这在测试特定代码片段时可以节省大量时间。

总的来说,PythonWin 是一个有用且易于使用的 IDE,大多数创建 Python 脚本的 ArcGIS 专业人士都使用 PythonWin。我发现 PythonWin 的缺点包括无法将脚本组织到项目中,以及缺少脚本中存在的变量列表,这在导航大型脚本时非常有帮助。

Aptana Studio 3

有时,大型编程社区的工具有时会令新脚本编写者感到畏惧,他们更关注于简单地创建一个可以节省 GIS 分析时间的脚本,而不是使用正确的编程工具。这让我想起了不熟练的计算机用户,他们觉得自己不需要顶级的计算机,因为他们只想浏览互联网和发送电子邮件。

然而,事实正好相反:对计算机不熟悉的用户使用易于使用的顶级计算机会更好,而经验丰富的计算机用户可以用上网本来满足需求。

对于程序员和脚本编写者来说,也是如此。有时,拥有一个功能强大的 IDE 会更有效,它实际上会使脚本编写者更有效率,而经验丰富的程序员可能只需要记事本。Aptana Studio 3 等 IDE 中包含的所有功能都会节省脚本编写者的时间,并且学习起来所需的时间非常少。

Aptana Studio 3 可在 aptana.com 获取。下载并运行提供的安装程序来安装它。选择一个默认的主项目文件夹,该文件夹可以包含所有脚本项目;为此书,我创建了一个名为C:\Projects的文件夹。对于创建的每个项目,Aptana 将创建一个包含每个项目信息的项目文件。在工作时使用 Aptana Studio,使用网络文件夹可能很有用,因为其他人可以使用他们各自的 Aptana 安装访问项目。

安装完成后,下一步是创建一个PyDev项目。转到文件菜单,选择新建,然后选择PyDev项目。在创建第一个项目时,必须将 Python 解释器添加到 Aptana 的 Python 路径中。Aptana 可以支持多个解释器;就我们的目的而言,一个就足够了。转到 PyDev 项目菜单底部,点击点击此处来配置一个解释器。当出现首选项/Python 解释器菜单时,请确保在左侧选择解释器-Python,然后在右上角的菜单中点击新建

Aptana Studio 3

选择新建后,将出现一个小的对话框,要求输入解释器的名称和可执行文件的路径。点击浏览并导航到包含python.exe的文件夹。使用 Aptana Studio 运行 Python 脚本时不会生成终端窗口,因为所有输出都重定向到 Aptana Studio 控制台。选择python.exe并点击确定。接下来,在选择解释器菜单中点击确定,然后在首选项菜单中点击确定。回到PyDev 项目菜单,为项目命名,并使用默认的工作区位置或自定义位置(例如,C:\Projects)。

所有这些配置只需进行一次;一旦完成,创建一个PyDev项目就只需要提供一个名称和位置。现在,与该项目相关的所有脚本都将始终列在左侧菜单中(PyDev 包资源管理器),这是一种非常强大的组织和脚本项目的方法。

确保 Aptana Studio 处于 PyDev 视角(在Windows/打开视角/其他菜单中,选择PyDev)将提供三个主要窗口–左侧的包资源管理器,中间的脚本窗口,以及右侧的大纲窗口,其中列出了脚本中包含的变量。点击右侧的任何一个变量,将脚本窗口移动到代码的该部分,使脚本导航变得快速。此外,我习惯在脚本窗口下方中间添加控制台窗口,以便显示脚本的输出。

打开的脚本在脚本窗口中都有一个标签,这使得在脚本之间切换变得容易。此外,根据需要,窗口可以被关闭,以给脚本窗口腾出更多空间。将鼠标悬停在脚本中的变量上,将弹出一个菜单,描述变量首次创建的位置,这在某些时候可以救命,因为有时很容易忘记哪个变量是哪个(除非,当然,它们根据前一章中描述的规则被清楚地命名;即使如此,有时也会很痛苦)。

IDE 概述

有许多其他 IDE,包括商业和免费的,都可用于 Python 编码。最终,每个 GIS 分析师都必须选择使他们感到高效和舒适的工具。随着编程成为他们日常工作流程中更大的一部分,这可能会改变。务必尝试几个不同的 IDE,以找到易于使用且直观的 IDE。

Python 文件夹结构

Python 的文件夹结构不仅包含 Python 解释器。在子文件夹中驻留了许多重要的脚本、数字链接库,甚至 C 语言模块。并非所有脚本都始终使用,但每个脚本都在使 Python 编程环境成为可能中扮演着角色。最重要的文件夹是site-packages文件夹,其中包含大多数将在 Python 脚本中导入的模块。

Python 文件夹结构

模块所在位置

在每个 Python 文件夹中都有一个名为Lib的文件夹,在该文件夹中有一个名为site-packages的文件夹。在我的机器上,该文件夹位于C:\Python27\ArcGIS10.2\Lib\site-packages.几乎所有的第三方模块都被复制到这个文件夹中,以便按需导入。对于我们来说,这个规则的主要例外是 ArcPy 模块,它存储在Program Files文件夹中的ArcGIS文件夹内(例如,C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy)。为了实现这一点,ArcGIS 安装程序调整了 Python 系统路径(使用 sys 模块),以便 arcPy 模块可导入。

使用 Python 的 sys 模块添加模块

Python 的 sys 模块是一个允许用户利用 Python 解释器内建的系统工具的模块。sys 模块中最有用的函数之一是 sys.path。它是一个文件路径列表,用户可以修改它以调整 Python 查找导入模块的位置,而无需管理员访问权限。

当 ArcGIS 10.2 安装程序安装 Python 2.7 时,安装程序利用 sys.path 函数将 C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy 添加到系统路径。为了测试这一点,启动 Python 解释器或 IDE,并输入以下内容:

>>> import sys
>>> print sys.path

输出如下:

['', 'C:\\WINDOWS\\SYSTEM32\\python27.zip', 'C:\\Python27\\ArcGIS10.2\\Dlls', 'C:\\Python27\\ArcGIS10.2\\lib', 'C:\\Python27\\ArcGIS10.2\\lib\\plat-win', 'C:\\Python27\\ArcGIS10.2\\lib\\lib-tk', 'C:\\Python27\\ArcGIS10.2\\Lib\\site-packages\\pythonwin', 'C:\\Python27\\ArcGIS10.2', 'C:\\Python27\\ArcGIS10.2\\lib\\site-packages', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.2\\bin', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.2\\arcpy', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.2\\ArcToolbox\\Scripts', 'C:\\Python27\\ArcGIS10.2\\lib\\site-packages\\win32', 'C:\\Python27\\ArcGIS10.2\\lib\\site-packages\\win32\\lib']

系统路径(存储在变量 sys.path 中)包括 ArcPy 所需的全部文件夹,以自动化 ArcGIS。系统路径包含 PYTHONPATH 环境变量中列出的所有目录(如果已创建);这不同于之前讨论的 Windows 路径环境变量。这两个独立的路径变量协同工作,以帮助 Python 定位模块。

sys.path.append() 方法

sys.path 函数是一个列表(你注意到前面代码输出中的方括号了吗?)因此可以添加或扩展以包含指向用户想要导入的模块的新文件路径。为了避免调整 sys.path 的需要,将模块复制到 site-packages 文件夹中。当这不可能时,使用 sys.path.append() 方法代替:

>>> sys.path.append("C:\\Projects\\Requests")
>>> sys.path
['', 'C:\\WINDOWS\\SYSTEM32\\python27.zip',
 'C:\\Python27\\ArcGIS10.2\\Dells',
 'C:\\Python27\\ArcGIS10.2\\lib',
..'C:\\Python27\\ArcGIS10.2\\lib\\plat-win',
..'C:\\Python27\\ArcGIS10.2\\lib\\lib-tk',
..'C:\\Python27\\ArcGIS10.2\\Lib\\site-packages\\pythonwin',
..'C:\\Python27\\ArcGIS10.2', ..'C:\\Python27\\ArcGIS10.2\\lib\\site-packages', 'C:\\Program
..Files (x86)\\ArcGIS\\Desktop10.2\\bin', 'C:\\Program Files
..(x86)\\ArcGIS\\Desktop10.2\\arcpy', 'C:\\Program Files
..(x86)\\ArcGIS\\Desktop10.2\\ArcToolbox\\Scripts',
..'C:\\Python27\\ArcGIS10.2\\lib\\site-packages\\win32',
..'C:\\Python27\\ArcGIS10.2\\lib\\site-packages\\win32\\lib',
..'C:\\Projects\\Requests']

当使用 sys.path.append() 方法时,调整是临时的。在 Windows 系统属性 菜单中调整 PYTHONPATH 环境变量(在路径环境变量部分讨论过)以进行永久更改(如果尚未创建,则创建 PYTHONPATH)。

最后一点是,为了在不调整系统路径或复制模块到 site-packages 文件夹的情况下导入模块,请将模块放置在包含导入脚本的文件夹中。只要模块配置正确,它就会正常工作。这在没有对机器的行政访问权限时非常有用。

摘要

在本章中,我们介绍了大量关于 Python 如何执行脚本和命令的知识,以及用于编写脚本的开发环境。特别是,我们讨论了 Python 脚本是如何被 Python 解释器读取和执行的,Python 解释器位于 Python 文件夹结构中的位置,以及不同的 Python 脚本扩展名(.py.pyc.pyw)分别代表什么。我们还介绍了集成开发环境以及它们之间的比较和对比。

在下一章中,我们将介绍如何使用 ModelBuilder 将建模分析转换为 Python 脚本,以及如何使其比导出版本更强大。

第三章。创建第一个 Python 脚本

现在我们已经将 Python 配置得符合我们的需求,我们可以创建 Python 脚本。本章将探讨如何使用 ArcGIS ModelBuilder来建模一个简单的分析作为脚本的基础。ModelBuilder 本身非常有用,并且对于创建 Python 脚本也非常有用,因为它具有操作和可视化组件,并且所有模型都可以输出为 Python 脚本。这将使我们能够比较更熟悉的 ModelBuilder 如何利用 ArcToolbox 中的工具与 Python 如何处理相同的工具。我们还将讨论迭代以及在何时使用 Python 而不是 ModelBuilder 最佳。

在本章中,我们将涵盖以下主题:

  • 使用 ModelBuilder 建模简单分析

  • 将模型导出到 Python 脚本

前提条件

“除了 ArcGIS ModelBuilder 外,还需要数据集和脚本。”

对于本章,相关的数据和脚本应从 Packt Publishing 的网站下载。完成的脚本可用于比较目的,数据将用于本章的分析。

ModelBuilder

ArcGIS 自 1970 年代以来一直在开发中。在这段时间里,它包含各种编程语言和工具,以帮助 GIS 分析师自动化分析和地图制作。这些包括 ArcGIS 3x 系列中的 Avenue 脚本语言以及 ARC/Info 工作站时代的ARC 宏语言(AML),以及直到 ArcGIS 10x 时的 VBScript,当时引入了 Python。ArcGIS 9x 中引入的另一项有用工具是 ModelBuilder,这是一个用于建模分析和创建可重复使用不同输入要素类的工具的可视化编程环境。

ModelBuilder 的另一个有用功能是导出功能,允许模型师直接从模型创建 Python 脚本。这将使比较 ModelBuilder 工具中的输入如何被接受与 Python 脚本如何调用同一工具并向其提供输入,或者创建的要素类如何命名和放置在文件结构中变得更加容易。ModelBuilder 是一个出色的工具,它将使 GIS 分析师轻松地从常规 GIS 工作流程过渡到基于 Python 的自动化工作流程。

创建模型并导出到 Python

本章将依赖于可从 Packt Publishing 网站下载的SanFrancisco.gdb文件地理数据库。旧金山 GDB 包含从data.sfgov.org和位于factfinder2.census.gov的美国人口普查局的美国事实查找网站下载的数据。地理数据库中包含的所有人口普查和地理数据均来自 2010 年人口普查。数据包含在名为SanFrancisco的要素数据集中。该要素数据集中的数据位于 NAD 83 加利福尼亚州平面坐标系 3 区,线性度量单位为美国英尺(这对应于欧洲石油调查组,或 EPSG,格式中的 SRID 2227)。

我们将使用模型创建的分析,并将其最终导出到 Python 中进行进一步优化,将使用旧金山的特定线路上的公交车站。这些公交车站将被缓冲以创建每个公交车站周围的代表性区域。然后,这些缓冲区域将与人口普查区相交,以找出每个代表性区域周围有多少人。

建模选择和缓冲工具

使用模型构建器,我们首先将构建公交车站分析的基础模型。一旦建模完成,它将被导出为自动生成的 Python 脚本。按照以下步骤开始分析:

  1. 打开ArcCatalog并创建一个文件夹连接到包含SanFrancisco.gdb的文件夹。右键单击地理数据库,添加一个名为Chapter3Tools的新工具箱。

  2. 接下来,打开模型构建器并创建一个模型,将其保存在Chapter3Tools工具箱中,命名为Chapter3Model1

  3. ArcToolbox中的分析/提取工具集拖动Bus_Stops要素类和选择工具。

  4. 打开选择工具,并将输出要素类命名为Inbound71。确保要素类被写入到Chapter3Results要素数据集中。

  5. 打开表达式SQL 查询构建器,创建以下 SQL 表达式:NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'建模选择和缓冲工具

  6. 下一步是添加来自分析/邻近工具集的缓冲工具。缓冲工具将用于创建每个公交车站周围的缓冲区。缓冲的公交车站允许我们与以人口普查区形式的人口普查数据进行交集,从而创建每个公交车站周围的代表性区域。

  7. 选择工具的输出(Inbound71)连接到缓冲工具。打开缓冲工具,并将 400 添加到距离字段,单位为英尺。其余选项留空。点击确定并返回模型。

建模选择和缓冲工具

添加交集工具

现在我们已经选择了感兴趣的公交线路,并将车站缓冲以创建代表性区域,我们需要将这些区域与人口普查区相交,以找出每个代表性区域的居民人数:

  1. 首先,将CensusBlocks2010要素类从SanFrancisco要素数据集添加到模型中。

  2. 接下来,添加位于ArcToolbox分析/叠加工具集的交集工具。虽然我们可以使用空间连接来实现类似的结果,但我使用交集工具来捕获交集区域,以便在模型和脚本中稍后使用。

到目前为止,我们的模型应该看起来像这样:

添加交集工具

统计分析结果

在创建了这个简单的分析之后,下一步是确定每个公交车站的结果。找出受每个公交车站 400 英尺缓冲区影响的普查区中居住的人数,需要检查最终要素类中的每一行数据,并选择与公交车站对应的行。一旦这些行被选中,就可以使用 字段计算器汇总 工具计算选中行的总和。所有这些方法都可行,但都不是完美的。它们花费的时间太长,更糟糕的是,如果模型中的假设被调整(例如,缓冲区从 400 英尺调整为 500 英尺),则这些方法无法自动重复。

这就是传统使用 ModelBuilder 的方式开始让分析师感到失败的地方。应该很容易指示模型选择与每个公交车站相关的所有行,并为每个公交车站的代表区域生成人口总和。如果模型能够创建一个包含分析最终结果的电子表格就更好了。现在是时候使用 Python 将这一分析提升到下一个层次了。

导出模型并调整脚本

虽然 ModelBuilder 中的建模分析有其缺点,但 ModelBuilder 中内置了一个非常棒的功能;即创建模型并将其导出到 Python 的能力。结合 ArcGIS 帮助文档,这是发现编写 ArcPy 脚本时应该使用正确 Python 语法的最有效方法。

SanFrancisco 地数据库旁边创建一个可以存放导出脚本的文件夹(例如,C:\Projects\Scripts)。这将包含 ArcGIS 自动生成的导出脚本,以及我们将从这些生成的脚本中构建的版本。

打开名为 Chapter3Model1 的模型,并在左上角的 模型 菜单中点击。从菜单中选择 导出,然后选择 导出到 Python 脚本。将脚本保存在脚本文件夹中,命名为 Chapter3Model1.py

注意

注意,还有将模型导出为图形的选项。创建模型的图形是一个很好的方法,可以在不共享模型和数据的情况下与其他分析师分享模型正在做什么,当共享 Python 脚本时也很有用。

自动生成的脚本

在 IDE 中打开自动生成的脚本。它应该看起来像这样:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8662_Chapter3Model1.py
# Created on: 2014-04-22 21:59:31.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
# Local variables:
Bus_Stops = "C:\\Projects\\PacktDB.gdb\\SanFrancisco\\Bus_Stops"
CensusBlocks2010 = "C:\\Projects\\PacktDB.gdb\\SanFrancisco\\CensusBlocks2010"
Inbound71 = "C:\\Projects\\PacktDB.gdb\\Chapter3Results\\Inbound71"
Inbound71_400ft_buffer = "C:\\Projects\\PacktDB.gdb\\Chapter3Results\\Inbound71_400ft_buffer"
Intersect71Census = "C:\\Projects\\PacktDB.gdb\\Chapter3Results\\Intersect71Census"
# Process: Select
arcpy.Select_analysis(Bus_Stops,
 Inbound71,
 "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'")
# Process: Buffer
arcpy.Buffer_analysis(Inbound71,
 Inbound71_400ft_buffer,
 "400 Feet", "FULL", "ROUND", "NONE", "")
# Process: Intersect
arcpy.Intersect_analysis("C:\\Projects\\PacktDB.gdb\\Chapter3Results\\Inbound71_400ft_buffer #;C:\\Projects\\PacktDB.gdb\\SanFrancisco\\CensusBlocks2010 #", Intersect71Census, "ALL", "", "INPUT")

让我们逐行检查这个脚本。第一行前面有一个井号(#),这意味着这一行是一个注释;然而,当脚本按常规执行时,Python 解释器不会忽略这一行,而是用它来帮助 Python 解释脚本编码,如这里所述:legacy.python.org/dev/peps/pep-0263

第二行注释和第三行注释是为了装饰目的。接下来的四行注释(全部注释)用于向读者提供有关脚本的信息,包括脚本的名称和创建时间,以及从模型属性中提取的描述。还包含了一行装饰性注释,用于在视觉上将信息性标题与脚本的主体分开。虽然将注释信息包含在脚本中供其他脚本用户使用是友好的,但并非必需。

脚本的主体或脚本的执行部分从import arcpy行开始。按照惯例,导入语句包含在脚本主体的顶部。在这种情况下,只导入了 ArcPy 模块。

ModelBuilder 的导出功能不仅创建了一个可执行脚本,还为每个部分添加了注释,以帮助标记脚本的各个部分。这些注释让用户知道变量所在的位置以及 ArcToolbox 工具正在执行的位置。随着读者对代码的理解加深,注释可能会变得多余,但 ESRI 包含这些注释是很友好的。

在导入语句下方是变量。在这种情况下,变量代表输入和输出要素类的文件路径。变量名称是从要素类的名称(文件路径的基本名称)派生出来的。文件路径通过赋值运算符(=)分配给变量,文件路径的部分通过两个反斜杠分隔。

Python 中的文件路径

很好地回顾一下在 Python 中与在 Windows 中表示方式相比,文件路径的使用方式。在 Python 中,文件路径是字符串,Python 中的字符串有特殊字符用于表示制表符(\t)、换行符(\n)或回车符(\r),以及其他许多字符。这些特殊字符都包含单个反斜杠,这使得创建使用单个反斜杠的文件路径变得非常困难。这本来不是什么大问题,但是 Windows 资源管理器中的文件路径全部使用单个反斜杠。

有多种方法可以避免这个问题。Python 是在 Linux 环境中开发的,那里的文件路径使用正斜杠。这种更符合 Python 风格的表示方式在 Windows 环境中使用 Python 时也是可用的,如下所示:

Windows Explorer: "C:\Projects\PacktDB.gdb\Chapter3Results\Intersect71Census"
Pythonic version:  "C:/Projects/PacktDB.gdb/Chapter3Results/Intersect71Census"

在 Python 脚本中,使用正斜杠的文件路径将正常工作,而 Windows 资源管理器版本可能会使脚本抛出异常。

避免特殊字符问题的另一种方法是 ModelBuilder 在从模型自动创建 Python 脚本时使用的方法。在这种情况下,反斜杠是通过使用第二个反斜杠来转义的。前面的脚本使用第二种方法产生了以下结果:

Python escaped version:  "C:\\Projects\\PacktDB.gdb\\Chapter3Results\\Intersect71Census"

第三个方法,我比较喜欢,是创建一个所谓的原始字符串。这与常规字符串相同,但在脚本开始前包含一个r。这个r会通知 Python 解释器,接下来的脚本不包含任何特殊字符或转义字符。以下是一个使用示例:

Python raw string:  r"C:\Projects\PacktDB.gdb\Chapter3Results\Intersect71Census"

使用原始字符串将使从 Windows 资源管理器获取文件路径并将其添加到脚本中的字符串内变得更加容易。它还将使避免在文件路径中意外忘记包含一组双反斜杠变得更加容易,这种情况经常发生,是许多脚本错误的原因。

继续分析脚本:ArcPy 工具

下一个,也是最重要的脚本部分是执行分析的地方。在这个部分中包含了我们在模型中创建的工具,即SelectBufferIntersect工具。模型中提供的相同参数也包含在这里:输入和输出,以及 Select 工具中的 SQL 语句和 Buffer 工具中的缓冲距离。

工具参数按照在模型工具界面中出现的顺序提供给脚本中的工具。以下是脚本中的Select工具:

arcpy.Select_analysis(Bus_Stops, Inbound71, "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'")

它的工作原理是这样的。arcPy 模块有一个名为Select_analysis的方法或特殊属性。当调用此方法时,需要三个参数:输入要素类(或 shapefile)、输出要素类和 SQL 语句。在这个例子中,输入由变量Bus_Stops表示,输出要素类由变量Inbound71表示,这两个变量都在变量部分定义。SQL 语句作为第三个参数包含在内。请注意,如果变量在上文定义,它也可以表示为变量;作为字符串的 SQL 语句可以分配给变量,然后变量可以替换 SQL 语句作为第三个参数。以下是一个使用变量替换参数的示例:

sqlStatement = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
arcpy.Select_analysis(Bus_Stops, Inbound71, sqlStatement)

虽然 ModelBuilder 在将输入和输出要素类分配给变量方面做得很好,但它不会将变量分配给参数的每个部分。当我们调整和构建自己的脚本时,这将是一个重要的纠正点。

Buffer 工具接受与 Select 工具相似的参数集。它有一个由变量表示的输入要素类,一个输出要素类变量,以及我们提供的距离(在这个例子中是 400 英尺),以及一系列默认提供的参数。请注意,这些参数依赖于关键字,并且这些关键字可以在脚本的文本中调整,以调整最终的缓冲输出。例如,英尺可以调整为米,缓冲区会更大。请查看工具的帮助部分,以更好地了解其他参数将如何影响缓冲区,以及 ArcPy 中Buffer工具将接受的键词参数。此外,如前所述,所有参数都可以分配给变量,如果在脚本中重复使用相同的参数,这可以节省时间。

有时提供的参数只是一个空字符串,就像这里最后一个参数的情况:

arcpy.Buffer_analysis(Inbound71,Inbound71_400ft_buffer,
 "400 Feet", "FULL", "ROUND", "NONE", "")

空字符串,在这种情况下表示这个缓冲区没有溶解字段,在 ArcPy 中非常常见。它也可以用两个单引号表示,但 ModelBuilder 已经被构建为使用双引号来包含字符串。

Intersect 工具和字符串操作

最后一个工具,Intersect 工具,在工具执行时使用不同的方法来表示需要相交的文件。因为这个工具接受输入部分中的多个文件(这意味着在一个操作中可以相交的文件数量没有限制),它将所有文件路径存储在一个字符串中。这个字符串使用井号或磅号(#)来分隔输入字符串中的文件路径。如果我们想在脚本工具中使用Intersect工具,我们必须处理这种轻微的偏差。如果我们从这个脚本构建一个工具,我们将在它们运行之前不知道将要相交的文件,因此我们需要了解将变量插入字符串的方法。

有三种方法可以将变量插入到字符串中。每种方法都有其技术上的优缺点。了解这三种方法都很好,因为它们的应用范围超出了我们在这里的需求,所以让我们来回顾一下。

字符串操作方法 1-字符串相加

字符串相加最初可能是一个奇特的概念,因为它看起来不可能将字符串相加,这与整数或浮点数不同,它们是数字。然而,在 Python 和其他编程语言中,这却是一个正常的步骤。使用加号(+),字符串可以相加形成更长的字符串,或者允许变量被添加到现有字符串的中间。以下是一些这个过程的示例:

>>> aString = "This is a string"
>>> bString = " and this is another string"
>>> aString + bString

输出如下:

'This is a string and this is another string'

>>> cString = aString + bString
>>> cString

输出如下:

'This is a string and this is another string'

两个或多个字符串可以相加,甚至可以分配给第三个变量。这个过程在诸如 Intersect 工具的输入字符串等情况下可能很有用。字符串可以被拆分,代表文件路径的变量可以插入到字符串的中间:

filePath1 = r"C:\Projects\Inbound71_400ft_buffer"
filePath2 =  r"C:\Projects\CensusBlocks2010"
arcpy.Intersect_analysis(filePath1 + " #;" + filePath2 + " #", Intersect71Census, "ALL", "", "INPUT")

这是一种强大且有用的方法,可以将文件路径插入到输入字符串中。只要分隔符仍然包含在字符串中,字符串仍然有效,Intersect 工具将按预期运行。以下是字符串相加完成后的字符串看起来像什么:

>>> filePath1 = r"C:\Projects\Inbound71_400ft_buffer"
>>> filePath2 =  r"C:\Projects\CensusBlocks2010"
>>> inputString = filePath1 + " #;" + filePath2 + " #"
>>> print inputString

输出如下:

C:\Projects\Inbound71_400ft_buffer #;C:\Projects\CensusBlocks2010 #

另一种与字符串相加类似的分支是字符串乘法,其中字符串通过一个整数相乘来产生字符串的重复版本:

>>>"string" * 3

输出如下:

'stringstringstring'

字符串操作方法 2–字符串格式化 #1

字符串操作的第二种方法,称为字符串格式化,涉及在字符串中添加将接受特定类型数据的占位符。这意味着这些特殊字符串可以接受其他字符串以及整数和浮点值。这些占位符使用取模 (%) 和一个键字母来指示预期的数据类型。字符串使用 %s 表示,浮点数使用 %f 表示,整数使用 %d 表示。浮点数还可以通过在取模后添加一个修改数字来调整,以限制包含的数字位数。如果字符串中有多个占位符,则值以元组的形式传递到字符串中。

自从 Python 2.6 中引入了下一节讨论的第三种方法以来,这种方法已经变得不那么流行了,但了解它仍然很有价值,因为许多较老的脚本仍在使用它。以下是这个方法的示例:

>>> origString = "This string has as a placeholder %s"
>>> newString = origString % "and this text was added"
>>> print newString

输出如下:

This string has as a placeholder and this text was added

这里是一个使用浮点占位符的示例:

>>> floatString1 = "This string has a float here: %f"
>>> newString = floatString1 % 1.0
>>> print newString

输出如下:

This string has a float here: 1.000000
>>> floatString2 = "This string has a float here: %.1f"
>>> newString2 = floatString2 % 1.0
>>> print newString2

输出如下:

This string has a float here: 1.0

这里是一个使用整数占位符的示例:

>>> intString = "Here is an integer: %d"
>>> newString = intString % 1
>>> print newString

输出如下:

Here is an integer: 1

对于 Intersect 工具,可以使用 %s 符号来接受文件路径字符串变量:

filePath1 = r"C:\Projects\Inbound71_400ft_buffer"
filePath2 =  r"C:\Projects\CensusBlocks2010"
arcpy.Intersect_analysis("%s #;%s #" % (filePath1,filePath2), Intersect71Census, "ALL", "", "INPUT")

字符串操作方法 3–字符串格式化 #2

最后一种方法,也是最近引入的,也称为字符串格式化。它与前面讨论的字符串格式化类似,但增加了不需要特定类型占位符的好处。占位符,或称为标记,只需按顺序排列即可被接受。格式化函数是字符串内建的;通过在字符串中添加 .format 并传递参数,字符串接受这些值:

>>> formatString = "This string has 3 tokens: {0}, {1}, {2}"
>>> newString = formatString.format("String", 2.5, 4)
>>> print newString
The output is as follows:
This string has 3 tokens: String, 2.5, 4

标记不需要在字符串中按顺序排列,甚至可以重复。顺序是从传递给 .format 函数的参数中派生出来的,该函数将值传递到字符串中。

对于 Intersect 工具,字符串格式化看起来像这样:

filePath1 = r"C:\Projects\Inbound71_400ft_buffer"
filePath2 =  r"C:\Projects\CensusBlocks2010"
arcpy.Intersect_analysis("{0} #;{1} #".format(filePath1,filePath2), Intersect71Census, "ALL", "", "INPUT")

第三种方法因为能够重复添加值并使其能够避免向特定占位符提供错误类型的数据,而不同于第二种方法,已经成为我字符串操作的首选方法。

调整脚本

现在是时候将自动生成的脚本进行调整,以适应我们的需求。我们希望脚本既能生成输出数据,又能分析数据并将结果汇总到电子表格中。这个电子表格将包含每个公交车站的平均人口值。平均值将从与公交车站周围的缓冲代表区域相交的每个普查区块中得出。将原始脚本保存为Chapter3Model1Modified.py

将 CSV 模块添加到脚本中

对于这个脚本,我们将使用CSV模块,这是一个用于创建逗号分隔值电子表格的有用模块。其简单的语法将使其成为创建脚本输出的有用工具。需要注意的是,当 ArcGIS for Desktop 安装时,它也会安装xlrdxlwt模块,分别用于读取或生成 Excel 电子表格。

在 import arcPy 行下方添加import csv。这将允许我们使用 csv 模块创建电子表格:

# Import arcpy module
import arcpy
import csv

下一个调整是对Intersect工具进行的。注意,输入字符串中包含的两个路径也在变量部分定义为变量。从输入字符串中删除文件路径,并用编号占位符替换它们,然后添加格式函数并供应变量作为占位符:

# Process: Intersect
arcpy.Intersect_analysis("{0} #;{1}#".format( ..............  Inbound71_400ft_buffer,CensusBlocks2010),
                         Intersect71Census, "ALL", "", "INPUT")

访问数据:使用游标

现在脚本已经到位,可以生成我们需要的原始数据,我们需要一种方法来访问Intersect工具输出特征类中存储的数据。这种访问将允许我们聚合代表每个公交车站的数据行。我们还需要某种东西来在内存中存储聚合数据,以便写入电子表格。

为了完成第二部分,我们将使用 Python 字典。为了完成第一部分,我们将使用 ArcPy 模块中内置的方法:数据访问搜索游标。

Python 字典将在Intersect工具下方添加。在 Python 中,使用花括号创建字典。将以下行添加到脚本中:

dataDictionary = {}

此脚本将使用公交车站 ID 作为字典的键。值将是列表,将包含与每个公交车站 ID 相关的所有人口值。添加以下行以生成数据游标:

with arcpy.da.SearchCursor(Intersect71Census, ["STOPID","POP10"]) as cursor:
    for row in cursor:
        busStopID = row[0]
        pop10 = row[1]
        if busStopID not in dataDictionary.keys():
            dataDictionary[busStopID] = [pop10]
        else:
            dataDictionary[busStopID].append(pop10)

这个迭代结合了 Python 和 ArcPy 中的几个想法。使用 with … as 语句创建一个变量(游标),该变量代表arcpy.da.SearchCursor对象。它也可以写成这样:

cursor = arcpy.da.SearchCursor(Intersect71Census, ["STOPID","POP10"])

注意

with … as 结构的优势在于,当迭代完成后,游标对象将从内存中删除,从而消除了对正在评估的特征类的锁定。

arcpy.da.SearchCursor() 函数需要一个输入要素类和一个要返回的字段列表。可选地,一个 SQL 语句可以限制返回的行数。

下一个行,for row in cursor,是遍历数据。它不是一个正常的 Python 遍历,这种区别在某些情况下会产生影响。例如,然而,它确实允许逐行访问包含在提供的要素类中的数据。请注意,当使用搜索游标时,每一行数据都作为元组返回,无法修改。可以使用索引访问数据,如前述代码所示,其中元组的两个成员被分配给变量。

if/else 条件语句允许对数据进行排序。如前所述,公交站 ID,即元组中包含的数据的第一个成员,将被用作键。条件语句评估公交站 ID 是否包含在字典的现有键中(这些键包含在一个列表中,并使用 dictionary.keys() 方法访问)。如果没有,它将被添加到键中,并分配一个包含(最初)一个数据片段的列表作为值,即该行中包含的人口值。如果它已经存在于键中,则将下一个与该公交站 ID 相关的人口值附加到列表中。使用此代码,我们现在已根据与之关联的公交站对每个普查区人口进行排序。

接下来,我们需要添加代码来创建电子表格。此代码将使用相同的 with ... as 结构,并使用两个内置 Python 函数 sumlen 来生成平均人口值,其中 sum 从数字列表创建总和,而 len 将获取列表、元组或字符串的长度:

with open(r'C:\Projects\Output\Averages.csv', 'wb') as csvfile:
    csvwriter = csv.writer(csvfile, delimiter=',')
    for busStopID in dataDictionary.keys():
        popList = dataDictionary[busStopID]
        averagePop = sum(popList)/len(popList)
        data = [busStopID, averagePop]
        csvwriter.writerow(data)

使用公交站 ID 键从字典中检索平均人口值,并将其分配给变量 averagePop。然后将两个数据片段,BusStopIDaveratePop 变量添加到一个列表中,该列表被提供给一个 CSVwriter 对象,该对象知道如何接受数据并将其写入文件,该文件位于提供给内置 Python open() 函数的文件路径中,用于创建简单文件。

脚本已完成,尽管在末尾添加一行以提供脚本已运行的视觉确认是件好事:

print "Data Analysis Complete"

这将创建一个输出,表明脚本已运行。一旦完成,请转到输出 csv 文件的所在位置,并使用 Excel 或记事本打开它,查看分析结果。我们的第一个脚本已完成!

最终脚本

这是脚本最终应该看起来的样子:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8662_Chapter3Model1.py
# Created on: 2014-04-22 21:59:31.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import csv

# Local variables:
Bus_Stops = r"C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops"
CensusBlocks2010 = r"C:\Projects\PacktDB.gdb\SanFrancisco\CensusBlocks2010"
Inbound71 = r"C:\Projects\PacktDB.gdb\Chapter3Results\Inbound71"
Inbound71_400ft_buffer = r"C:\Projects\PacktDB.gdb\Chapter3Results\Inbound71_400ft_buffer"
Intersect71Census = r"C:\Projects\PacktDB.gdb\Chapter3Results\Intersect71Census"

# Process: Select
arcpy.Select_analysis(Bus_Stops,
 Inbound71,
 "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'")
# Process: Buffer
arcpy.Buffer_analysis(Inbound71,
 Inbound71_400ft_buffer,
 "400 Feet", "FULL", "ROUND", "NONE", "")
# Process: Intersect
arcpy.Intersect_analysis("{0} #;{1} #".format(Inbound71_400ft_buffer,CensusBlocks2010),
 Intersect71Census, "ALL", "", "INPUT")

dataDictionary = {}

with arcpy.da.SearchCursor(Intersect71Census, ["STOPID","POP10"]) as cursor:
 for row in cursor:
 busStopID = row[0]
 pop10 = row[1]
 if busStopID not in dataDictionary.keys():
 dataDictionary[busStopID] = [pop10]
 else:
 dataDictionary[busStopID].append(pop10)

with open(r'C:\Projects\Output\Averages2.csv', 'wb') as csvfile:
 spamwriter = csv.writer(csvfile, delimiter=',')
 for busStopID in dataDictionary.keys():
 popList = dataDictionary[busStopID]
 averagePop = sum(popList)/len(popList)
 data = [busStopID, averagePop]
 spamwriter.writerow(data)

print "Data Analysis Complete"

摘要

在本章中,我们介绍了如何构建分析模型并将其导出为脚本。在讨论脚本之后,我们对脚本进行了调整,以包含结果分析和总结,这些内容被输出到 CSV 文件中。特别是,我们讨论了如何使用 ModelBuilder 创建分析和将其导出为脚本,以及如何调整脚本使其更符合 Python 风格。我们还简要提到了搜索游标的使用,这将在第五章(Chapter 5)中更详细地介绍,即ArcPy 游标 – 搜索、插入和更新。此外,我们还看到了如何使用内置模块,如 CSV 模块,与 ArcPy 一起使用,以捕获格式化的电子表格中的分析输出。

在下一章中,我们将讨论如何创建更复杂的脚本并构建函数以避免代码重复。这些函数将使得代码一旦编写就可以永久使用。这种代码的重用将展示 Python 如何超越分析自动化,成为一套新的生产力工具集。

第四章。复杂 ArcPy 脚本和泛化函数

在本章中,我们将从基于 ModelBuilder 自动生成的简单脚本创建转向包含高级 Python 和 ArcPy 概念的复杂脚本,例如函数。函数可以在编写脚本时提高代码效率和节省时间。当创建模块或其他可重用代码时,它们也非常有用,允许将标准编程操作脚本化并准备好供将来使用。

本章将涵盖以下主题:

  • 创建函数以避免代码重复

  • 创建辅助函数以处理 ArcPy 的限制

  • 将函数泛化以使其可重用

Python 函数–避免代码重复

编程语言共享一个几十年来帮助程序员的观念:函数。从广义上讲,函数的概念是创建代码块,将对数据进行操作,根据程序员的所需对其进行转换,并将转换后的数据返回到代码的主体部分。在前几章中,我们已经介绍了 Python 的一些内置函数,例如 int 函数,它将字符串或浮点数转换为整数;现在是我们编写自己的函数的时候了。

函数之所以被使用,是因为它们在编程中解决了许多不同的需求。函数减少了编写重复代码的需要,从而减少了创建脚本所需的时间。它们可以用来创建数字范围(range() 函数),或者确定列表的最大值(max 函数),或者创建一个 SQL 语句来从要素类中选择一组行。它们甚至可以被复制并用于另一个脚本,或者作为可以导入到脚本中的模块的一部分。函数的重用不仅使编程更有用,而且减少了繁琐的工作。当脚本编写者开始编写函数时,这是将编程作为 GIS 工作流程一部分的重大步骤。

函数的技术定义

函数,在其他编程语言中也称为子程序或过程,是一段代码块,其设计目的是接受输入数据并对其进行转换,或者在没有输入要求的情况下调用时,向主程序提供数据。在理论上,函数只会转换提供给函数作为参数的数据;它不应更改函数中未包含的脚本的其他任何部分。为了实现这一点,引入了命名空间的概念。如第一章中所述,“ArcGIS 的 Python 简介”,命名空间用于隔离脚本中的变量;变量要么是全局的,可以在脚本的主体以及函数中使用,要么是局部的,仅在函数内部可用。

命名空间使得在函数内部使用变量名成为可能,并允许它表示一个值,同时也在脚本的另一部分使用相同的变量名。当从其他程序员那里导入模块时,这一点尤为重要;在该模块及其函数内部,它包含的变量可能具有与主脚本中变量名相同的名称。

在像 Python 这样的高级编程语言中,有内置对函数的支持,包括定义函数名称和数据输入(也称为参数)的能力。函数是通过使用def关键字加上函数名称以及可能包含或不包含参数的括号来创建的。参数也可以定义默认值,因此只有当参数与默认值不同时才需要传递给函数。函数返回的值也容易定义。

第一个函数

让我们创建一个函数,以了解编写函数时可能实现的功能。首先,我们需要通过提供def关键字和括号内的名称来调用函数。当调用firstFunction()时,它将返回一个字符串:

def firstFunction():
 'a simple function returning a string'
 return "My First Function"
>>>firstFunction()

输出如下:

'My First Function'

注意,这个函数有一个文档字符串或 doc 字符串(一个简单返回字符串的函数),它描述了函数的功能;这个字符串可以在以后被调用,以了解函数的功能,使用__doc__内部函数:

>>>print firstFunction.__doc__

输出如下:

'a simple function returning a string' 

函数被定义并赋予了一个名称,然后添加括号并跟一个冒号。接下来的行必须进行缩进(一个好的 IDE 会自动添加缩进)。该函数没有任何参数,因此括号是空的。然后函数使用return关键字从函数返回一个值,在这种情况下是一个字符串。

接下来,通过在函数名称后添加括号来调用函数。当它被调用时,它将执行它被指示执行的操作:返回一个字符串。

带参数的函数

现在,让我们创建一个接受参数并根据需要转换它们的函数。这个函数将接受一个数字并将其乘以 3:

def secondFunction(number):
 'this function multiples numbers by 3'
 return number *3
>>> secondFunction(4)

输出如下:

12

然而,该函数有一个缺点;无法保证传递给函数的值是一个数字。我们需要在函数中添加一个条件,以确保它不会抛出异常:

def secondFunction(number):
 'this function multiples numbers by 3'
 if type(number) == type(1) or type(number) == type(1.0):
 return number *3
>>> secondFunction(4.0)

输出如下:

12.0
>>>secondFunction(4)

输出如下:

12
>>>secondFunction("String")
>>> 

现在函数接受一个参数,检查它的数据类型,并返回参数的倍数,无论是整数还是函数。如果它是一个字符串或其他数据类型,如最后一个示例所示,则不返回任何值。

我们应该讨论的简单函数的另一个调整是参数默认值。通过在函数定义中包含默认值,我们避免了提供很少更改的参数。例如,如果我们想在简单函数中使用不同于 3 的乘数,我们可以这样定义它:

def thirdFunction(number, multiplier=3):
 'this function multiples numbers by 3'
 if type(number) == type(1) or type(number) == type(1.0):
 return number *multiplier
>>>thirdFunction(4)

输出如下:

12
>>>thirdFunction(4,5)

输出如下:

20

当只提供要乘以的数字时,该函数将正常工作,因为乘数默认值为 3。然而,如果我们需要另一个乘数,可以在调用函数时添加另一个值来调整其值。请注意,第二个值不必是数字,因为它没有类型检查。此外,函数中的默认值(s)必须跟在无默认值的参数之后(或者所有参数都可以有默认值,并且可以按顺序或按名称向函数提供参数)。

这些简单的函数结合了我们之前章节中讨论的许多概念,包括内置函数如 type条件语句参数参数默认值函数返回值。我们现在可以继续使用 ArcPy 创建函数。

使用函数替换重复代码

函数的主要用途之一是确保相同的代码不必反复编写。让我们回到上一章的例子,并将脚本中的代码转换为函数,以便能够对旧金山的任何公交线路进行相同的分析。

我们可以将脚本的第一部分转换为函数的是三个 ArcPy 函数。这样做将使脚本适用于 Bus Stop 特征类中的任何站点,并具有可调整的缓冲距离:

bufferDist = 400
buffDistUnit  = "Feet"
lineName = '71 IB'
busSignage = 'Ferry Plaza'
sqlStatement = "NAME = '{0}' AND BUS_SIGNAG = '{1}'"
def selectBufferIntersect(selectIn,selectOut,bufferOut,intersectIn, intersectOut, sqlStatement, bufferDist, buffDistUnit, lineName, busSignage): 
 'a function to perform a bus stop analysis'
 arcpy.Select_analysis(selectIn, selectOut, sqlStatement.format(lineName, busSignage))
 arcpy.Buffer_analysis(selectOut, bufferOut, "{0} {1}".format(bufferDist), "FULL", "ROUND", "NONE", "")
 arcpy.Intersect_analysis("{0} #;{1} #".format(bufferOut, intersectIn), intersectOut, "ALL", "", "INPUT")
 return intersectOut

此函数演示了如何调整分析以接受输入和输出特征类变量作为参数,以及一些新变量。

函数添加了一个变量来替换 SQL 语句,以及变量来调整公交站,并且调整了缓冲距离语句,以便可以调整距离和单位。在脚本中先前定义的特征类名称变量都已替换为局部变量名;虽然全局变量名可以保留,但这会降低函数的可移植性。

下一个函数将接受 selectBufferIntersect() 函数的结果,并使用搜索游标对其进行搜索,将结果传递到字典中。然后,该函数将返回字典供以后使用:

def createResultDic(resultFC):
 'search results of analysis and create results dictionary' 
 dataDictionary = {} 
 with arcpy.da.SearchCursor(resultFC, ["STOPID","POP10"]) as cursor:
 for row in cursor:
 busStopID = row[0]
 pop10 = row[1]
 if busStopID not in dataDictionary.keys():
 dataDictionary[busStopID] = [pop10]
 else:
 dataDictionary[busStopID].append(pop10)
 return dataDictionary

此函数仅需要一个参数:来自 searchBufferIntersect() 函数返回的特征类。首先创建一个持有字典的结果,然后通过搜索游标填充,使用 busStopid 属性作为键,并将人口普查区块人口属性添加到分配给该键的列表中。

字典在被填充了排序后的数据后,从函数中返回,用于最终函数createCSV()。这个函数接受字典和输出 CSV 文件的名称作为字符串:

def createCSV(dictionary, csvname):
  'a function takes a dictionary and creates a CSV file'
    with open(csvname, 'wb') as csvfile:
        csvwriter = csv.writer(csvfile, delimiter=',')
        for busStopID in dictionary.keys():
            popList = dictionary[busStopID]
            averagePop = sum(popList)/len(popList)
            data = [busStopID, averagePop]
            csvwriter.writerow(data)

最终的函数使用csv模块创建 CSV 文件。现在文件名,一个字符串,是一个可定制参数(这意味着脚本名称可以是任何有效的文件路径和具有.csv扩展名的文本文件)。csvfile参数传递给 CSV 模块的 writer 方法,并分配给变量csvwriter,然后访问和处理字典,并将其作为列表传递给csvwriter以写入CSV文件。csv.writer()方法将列表中的每个项目处理成 CSV 格式,并保存最终结果。使用 Excel 或记事本等文本编辑器打开CSV文件。

要运行这些函数,我们将在脚本中调用它们,在函数定义之后:

analysisResult = selectBufferIntersect(Bus_Stops,Inbound71, Inbound71_400ft_buffer,CensusBlocks2010, Intersect71Census, bufferDist, lineName, busSignage )
dictionary = createResultDic(analysisResult)
createCSV(dictionary,r'C:\Projects\Output\Averages.csv')

现在,脚本已经被分为三个函数,它们替换了第一个修改后的脚本中的代码。修改后的脚本看起来像这样:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8662_Chapter4Modified1.py
# Created on: 2014-04-22 21:59:31.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# Adjusted by Silas Toms
# 2014 05 05
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import csv

# Local variables:
Bus_Stops = r"C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops"
CensusBlocks2010 = r"C:\Projects\PacktDB.gdb\SanFrancisco\CensusBlocks2010"
Inbound71 = r"C:\Projects\PacktDB.gdb\Chapter3Results\Inbound71"
Inbound71_400ft_buffer = r"C:\Projects\PacktDB.gdb\Chapter3Results\Inbound71_400ft_buffer"
Intersect71Census = r"C:\Projects\PacktDB.gdb\Chapter3Results\Intersect71Census"
bufferDist = 400
lineName = '71 IB'
busSignage = 'Ferry Plaza'
def selectBufferIntersect(selectIn,selectOut,bufferOut,intersectIn,
 intersectOut, bufferDist,lineName, busSignage ):
 arcpy.Select_analysis(selectIn, 
 selectOut, 
 "NAME = '{0}' AND BUS_SIGNAG = '{1}'".format(lineName, busSignage))
 arcpy.Buffer_analysis(selectOut, 
 bufferOut, 
 "{0} Feet".format(bufferDist), 
 "FULL", "ROUND", "NONE", "")
 arcpy.Intersect_analysis("{0} #;{1} #".format(bufferOut,intersectIn), 
 intersectOut, "ALL", "", "INPUT")
 return intersectOut

def createResultDic(resultFC):
 dataDictionary = {}

 with arcpy.da.SearchCursor(resultFC, 
 ["STOPID","POP10"]) as cursor:
 for row in cursor:
 busStopID = row[0]
 pop10 = row[1]
 if busStopID not in dataDictionary.keys():
 dataDictionary[busStopID] = [pop10]
 else:
 dataDictionary[busStopID].append(pop10)
 return dataDictionary

def createCSV(dictionary, csvname):
 with open(csvname, 'wb') as csvfile:
 csvwriter = csv.writer(csvfile, delimiter=',')
 for busStopID in dictionary.keys():
 popList = dictionary[busStopID]
 averagePop = sum(popList)/len(popList)
 data = [busStopID, averagePop]
 csvwriter.writerow(data)
analysisResult = selectBufferIntersect(Bus_Stops,Inbound71, Inbound71_400ft_buffer,CensusBlocks2010,Intersect71Census, bufferDist,lineName, busSignage )
dictionary = createResultDic(analysisResult)
createCSV(dictionary,r'C:\Projects\Output\Averages.csv')
print "Data Analysis Complete"

函数的进一步泛化,虽然我们已经从原始脚本中创建了可以用来提取更多关于旧金山公交车站数据的函数,但我们的新函数仍然非常特定于它们被创建的数据集和分析。这对于那些不需要创建可重用函数的漫长而繁重分析来说非常有用。函数的第一个用途是消除重复代码的需要。接下来的目标是使代码可重用。让我们讨论一些方法,我们可以将这些函数从一次性函数转换为可重用函数甚至模块。

首先,让我们检查第一个函数:

def selectBufferIntersect(selectIn,selectOut,bufferOut,intersectIn,
 intersectOut, bufferDist,lineName, busSignage ):
 arcpy.Select_analysis(selectIn, 
 selectOut, 
 "NAME = '{0}' AND BUS_SIGNAG = '{1}'".format(lineName, busSignage))
 arcpy.Buffer_analysis(selectOut, 
 bufferOut, 
 "{0} Feet".format(bufferDist), 
 "FULL", "ROUND", "NONE", "")
 arcpy.Intersect_analysis("{0} #;{1} #".format(bufferOut,intersectIn), 
 intersectOut, "ALL", "", "INPUT")
 return intersectOut

这个函数似乎非常特定于公交车站分析。实际上,它如此特定,以至于尽管有几种方法可以调整它使其更通用(即,在其他可能不涉及相同步骤的脚本中也有用),我们不应该将其转换为单独的函数。当我们创建一个单独的函数时,我们引入了太多的变量到脚本中,试图简化它,这是一种事倍功半的努力。相反,让我们专注于泛化 ArcPy 工具本身的方法。

第一步将是将三个 ArcPy 工具分开,并检查每个工具可以调整什么。Select 工具应该调整以接受字符串作为 SQL 选择语句。SQL 语句可以由另一个函数生成,或者通过在运行时接受的参数生成(例如,通过 Script 工具传递给脚本,这将在下一章中讨论)。

例如,如果我们想使脚本接受每次运行脚本时多个公交车站(例如,每条线路的进站和出站车站),我们可以创建一个函数,该函数接受所需车站的列表和 SQL 模板,并返回一个 SQL 语句,可以插入到 Select 工具中。以下是一个示例:

def formatSQLIN(dataList, sqlTemplate):
 'a function to generate a SQL statement'
 sql = sqlTemplate #"OBJECTID IN "
 step = "("
 for data in dataList:
 step += str(data)
 sql += step + ")"
 return sql

def formatSQL(dataList, sqlTemplate):
 'a function to generate a SQL statement'
 sql = ''
 for count, data in enumerate(dataList):
 if count != len(dataList)-1:
 sql += sqlTemplate.format(data) + ' OR '
 else:
 sql += sqlTemplate.format(data)
 return sql

>>> dataVals = [1,2,3,4]
>>> sqlOID = "OBJECTID = {0}"
>>> sql = formatSQL(dataVals, sqlOID)
>>> print sql

输出如下:

OBJECTID = 1 OR OBJECTID = 2 OR OBJECTID = 3 OR OBJECTID = 4

这个新的函数formatSQL()是一个非常实用的函数。让我们通过将函数与其后的结果进行比较来回顾它所做的工作。该函数被定义为接受两个参数:一个值列表和一个 SQL 模板。第一个局部变量是空字符串sql,它将通过字符串添加来添加。该函数旨在将值插入到变量sql中,通过使用字符串格式化将它们添加到模板中,然后将模板添加到 SQL 语句字符串中(注意sql +=等同于sql = sql +)。此外,使用一个运算符(OR)使 SQL 语句包含所有匹配该模式的数据行。此函数使用内置的enumerate函数来计数列表的迭代次数;一旦它达到了列表中的最后一个值,运算符就不会添加到 SQL 语句中。

注意,我们还可以向函数添加一个参数,使其能够使用AND运算符而不是OR,同时仍然保留OR作为默认选项:

def formatSQL2(dataList, sqlTemplate, operator=" OR "):
 'a function to generate a SQL statement'
 sql = ''
 for count, data in enumerate(dataList):
 if count != len(dataList)-1:
 sql += sqlTemplate.format(data) + operator
 else:
 sql += sqlTemplate.format(data)
 return sql

>>> sql = formatSQL2(dataVals, sqlOID," AND ")
>>> print sql

输出如下:

OBJECTID = 1 AND OBJECTID = 2 AND OBJECTID = 3 AND OBJECTID = 4

虽然在 ObjectIDs 上使用AND运算符没有意义,但还有其他情况下使用它是有意义的,因此保留OR作为默认选项,同时允许使用AND。无论如何,这个函数现在可以用来生成我们的多站公交车站 SQL 语句(现在忽略公交标志字段):

>>> sqlTemplate = "NAME = '{0}'"
>>> lineNames = ['71 IB','71 OB']
>>> sql = formatSQL2(lineNames, sqlTemplate)
>>> print sql

输出如下:

NAME = '71 IB' OR NAME = '71 OB'

然而,我们不能忽视进站线路的公交标志字段,因为该线路有两个起点,所以我们需要调整函数以接受多个值:

def formatSQLMultiple(dataList, sqlTemplate, operator=" OR "):
    'a function to generate a SQL statement'
    sql = ''
    for count, data in enumerate(dataList):
        if count != len(dataList)-1:
            sql += sqlTemplate.format(*data) + operator
        else:
            sql += sqlTemplate.format(*data)
    return sql
>>> sqlTemplate = "(NAME = '{0}' AND BUS_SIGNAG = '{1}')"
>>> lineNames = [('71 IB', 'Ferry Plaza'),('71 OB','48th Avenue')]
>>> sql = formatSQLMultiple(lineNames, sqlTemplate)
>>> print sql

输出如下:

(NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza') OR (NAME = '71 OB' AND BUS_SIGNAG = '48th Avenue')

这个函数的细微差别在于数据变量前的星号,它允许数据变量内的值通过解包元组被正确地格式化到 SQL 模板中。注意,SQL 模板已被创建,通过使用括号来分隔每个条件。现在函数(们)可以重用,SQL 语句现在可以插入到选择工具中:

sql = formatSQLMultiple(lineNames, sqlTemplate)
arcpy.Select_analysis(Bus_Stops, Inbound71, sql)

接下来是缓冲工具。我们已经通过添加一个距离变量来使其通用化。在这种情况下,我们只需再添加一个变量,即单位变量,这将使缓冲单元可以从英尺调整到米或其他任何允许的单位。我们将保留其他默认设置不变。

这里是调整后的缓冲工具版本:

bufferDist = 400
bufferUnit = "Feet"
arcpy.Buffer_analysis(Inbound71, 
 Inbound71_400ft_buffer, 
 "{0} {1}".format(bufferDist, bufferUnit), 
 "FULL", "ROUND", "NONE", "")

现在,缓冲距离和缓冲单元都由前一个脚本中定义的变量控制,如果决定距离不足且可能需要调整变量,这将使其易于调整。

调整 ArcPy 工具的下一步是编写一个函数,该函数将允许使用交集工具将任意数量的特征类相互交集。这个新函数将与之前的formatSQL函数类似,因为它们将使用字符串格式化和附加来允许将特征类列表处理成交集工具可以接受的正确字符串格式。然而,由于这个函数将被构建成尽可能通用,它必须设计成可以接受任意数量的要交集的特征类:

def formatIntersect(features):
    'a function to generate an intersect string'
    formatString = ''
    for count, feature in enumerate(features):
        if count != len(features)-1:
            formatString += feature + " #;"
        else:
            formatString += feature + " #"
        return formatString
>>> shpNames = ["example.shp","example2.shp"]
>>> iString = formatIntersect(shpNames)
>>> print iString

输出结果如下:

example.shp #;example2.shp #

现在我们已经编写了formatIntersect()函数,需要创建的是要传递给函数的特征类列表。函数返回的字符串然后可以传递给交集工具:

intersected = [Inbound71_400ft_buffer, CensusBlocks2010]
iString = formatIntersect(intersected)
# Process: Intersect
arcpy.Intersect_analysis(iString, 
                         Intersect71Census, "ALL", "", "INPUT")

由于我们避免创建仅适用于此脚本或分析的函数,我们现在有两个(或更多)有用的函数可以在后续分析中使用,并且我们知道如何操作 ArcPy 工具以接受我们想要提供给它们的任何数据。

函数的更一般化

我们最初创建的用于搜索结果并生成结果电子表格的其他函数,也可以通过一些调整变得更加通用。

如果我们想要生成关于距离公交站点的每个普查区更详细的信息(例如,如果我们有一个包含收入数据和人口数据的普查区数据集),我们需要将一个属性列表传递给函数,以便从最终的特征类中提取这些属性。为了实现这一点,需要调整createResultDic()函数以接受这个属性列表:

def createResultDic(resultFC, key, values):
    dataDictionary = {}
    fields = [key]
    fields.extend(values)
    with arcpy.da.SearchCursor(resultFC, fields) as cursor:
        for row in cursor:
            busStopID = row[0]
            data = row[1:]
            if busStopID not in dataDictionary.keys():
                dataDictionary[busStopID] = [data]
            else:
                dataDictionary[busStopID].append(data)
    return dataDictionary

这个新的createResultDic()函数将为每个公交站生成一个列表的列表(即,每行的值包含在一个列表中,并添加到主列表中),然后可以通过知道列表中每个值的位位置来稍后解析。这种解决方案在需要将数据排序到字典中时很有用。

然而,这是一种不令人满意的结果排序方式。如果字段列表没有传递给字典,并且无法知道列表中数据的顺序怎么办?相反,我们希望能够使用 Python 字典的功能按字段名称对数据进行排序。在这种情况下,我们将使用嵌套字典来创建通过它们包含的数据类型(即,人口、收入或其他字段)可访问的结果列表:

def createResultDic(resultFC, key, values):
    dataDic = {}
    fields = []
   if type(key) == type((1,2)) or type(key) == type([1,2]):
            fields.extend(key)
            length = len(key)
    else:
        fields = [key]
        length = 1
    fields.extend(values)
    with arcpy.da.SearchCursor(resultFC, fields) as cursor:
        for row in cursor:
            busStopID = row[:length]
            data = row[length:]
            if busStopID not in dataDictionary.keys():

                dataDictionary[busStopID] = {}

            for counter,field in enumerate(values):
                if field not in dataDictionary[busStopID].keys():
                    dataDictionary[busStopID][field] = [data[counter]]
                else:
                    dataDictionary[busStopID][field].append(data[counter])
    return dataDictionary
>>> rFC = r'C:\Projects\PacktDB.gdb\Chapter3Results\Intersect71Census'
>>> key = 'STOPID'
>>> values = 'HOUSING10','POP10'
>>> dic = createResultDic(rFC, key, values)
>>> dic[1122023]

输出结果如下:

{'HOUSING10': [104, 62, 113, 81, 177, 0, 52, 113, 0, 104, 81, 177, 52], 'POP10': [140, 134, 241, 138, 329, 0, 118, 241, 0, 140, 138, 329, 118]}

在这个例子中,函数作为参数传递给一个要素类,STOPID以及要合并的字段。fields变量被创建出来,以便将所需的字段传递给搜索光标。光标返回每一行作为一个元组;元组的第一个成员是busStopID,其余的元组是该公交车站相关联的数据。然后,函数使用一个条件来评估该公交车站是否已经被分析过;如果没有,它将被添加到字典中,并分配一个第二级内部字典,该字典将用于存储与该站点相关的结果。通过使用字典,我们可以然后对结果进行排序,并将它们分配给正确的字段。

之前的例子展示了请求一个特定公交车站(1122023)的数据结果。由于这里传递了两个字段,数据已经被组织成两组,字段名称现在成为内部字典的键。正因为这种组织方式,我们现在可以为每个字段创建平均值,而不仅仅是单个值。

说到平均值,我们将搜索光标分析结果的平均值计算工作留给了createCSV()函数。这也应该避免,因为它通过添加额外的数据处理职责(这些职责应该由另一个函数负责)降低了createCSV()函数的有用性。让我们首先通过调整createCSV()函数来解决这个问题:

def createCSV(data, csvname, mode ='ab'):
    with open(csvname, mode) as csvfile:
        csvwriter = csv.writer(csvfile, delimiter=',')
        csvwriter.writerow(data)

这是一个简化版的函数,但它非常有用。通过这样调整函数,我们限制它只做两件事:打开 CSV 文件并向其中添加一行数据。因为我们使用了ab模式,如果 CSV 文件存在,我们只会向其中添加数据而不是覆盖它(如果不存在,它将被创建)。这种添加模式可以通过传递wb作为模式来覆盖,这将每次生成一个新的脚本。

现在我们可以对分析结果进行排序,计算平均值,并将它们传递给我们的新createCSV脚本。为此,我们将遍历由createResultDic()函数创建的字典:

csvname = r'C:\Projects\Output\Averages.csv'
dataKey = 'STOPID'
fields = 'HOUSING10','POP10'
dictionary = createResultDic(Intersect71Census, dataKey, fields)

header = [dataKey]
for field in fields:
    header.append(field)

createCSV(header,csvname, 'wb' )

for counter, busStop in enumerate(dictionary.keys()):
    datakeys  = dictionary[busStop]
    averages = [busStop]
    for key in datakeys:
        data = datakeys[key]
        average = sum(data)/len(data)
        averages.append(average)
    createCSV(averages,csvname)

最后一步展示了如何创建 CSV 文件:通过遍历字典中的数据,然后为每个公交车站的平均值。然后,这些平均值被添加到一个列表中,该列表包含每个公交车站的名称(以及在这个例子中它所属的线路),然后传递给createCSV()函数以写入到CSV文件中。

这里是最终的代码。请注意,我已经将许多自动生成的注释转换成了打印语句,以提供关于脚本状态的反馈:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8662_Chapter4Modified2.py
# Created on: 2014-04-22 21:59:31.00000
#   (generated by ArcGIS/ModelBuilder)
# Description: 
# Adjusted by Silas Toms
# 2014 04 23
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
import csv

Bus_Stops = r"C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops"
CensusBlocks2010 = r"C:\Projects\PacktDB.gdb\SanFrancisco\CensusBlocks2010"
Inbound71 = r"C:\Projects\PacktDB.gdb\Chapter4Results\Inbound71"
Inbound71_400ft_buffer = r"C:\Projects\PacktDB.gdb\Chapter4Results\Inbound71_400ft_buffer"
Intersect71Census = r"C:\Projects\PacktDB.gdb\Chapter4Results\Intersect71Census"
bufferDist = 400
bufferUnit = "Feet"
lineNames = [('71 IB', 'Ferry Plaza'),('71 OB','48th Avenue')]
sqlTemplate = "NAME = '{0}' AND BUS_SIGNAG = '{1}'"
intersected = [Inbound71_400ft_buffer, CensusBlocks2010]
dataKey = 'NAME','STOPID'
fields = 'HOUSING10','POP10'
csvname = r'C:\Projects\Output\Averages.csv'

def formatSQLMultiple(dataList, sqlTemplate, operator=" OR "):
    'a function to generate a SQL statement'
    sql = ''
    for count, data in enumerate(dataList):
        if count != len(dataList)-1:
            sql += sqlTemplate.format(*data) + operator
        else:
            sql += sqlTemplate.format(*data)
    return sql

def formatIntersect(features):
    'a function to generate an intersect string'
    formatString = ''
    for count, feature in enumerate(features):
        if count != len(features)-1:
            formatString += feature + " #;"
        else:
            formatString += feature + " #"
    return formatString

def createResultDic(resultFC, key, values):
    dataDictionary = {}
    fields = []
    if type(key) == type((1,2)) or type(key) == type([1,2]):
        fields.extend(key)
        length = len(key)
    else:
        fields = [key]
        length = 1
    fields.extend(values)
    with arcpy.da.SearchCursor(resultFC, fields) as cursor:
        for row in cursor:
            busStopID = row[:length]
            data = row[length:]
            if busStopID not in dataDictionary.keys():

                dataDictionary[busStopID] = {}

            for counter,field in enumerate(values):
                if field not in dataDictionary[busStopID].keys():
                    dataDictionary[busStopID][field] = [data[counter]]
                else:
                    dataDictionary[busStopID][field].append(data[counter])

    return dataDictionary

def createCSV(data, csvname, mode ='ab'):
    with open(csvname, mode) as csvfile:
        csvwriter = csv.writer(csvfile, delimiter=',')
        csvwriter.writerow(data)

sql = formatSQLMultiple(lineNames, sqlTemplate)

print 'Process: Select'
arcpy.Select_analysis(Bus_Stops, 
                      Inbound71, 
                      sql)

print 'Process: Buffer'
arcpy.Buffer_analysis(Inbound71, 
                      Inbound71_400ft_buffer, 
                      "{0} {1}".format(bufferDist, bufferUnit), 
                      "FULL", "ROUND", "NONE", "")

iString = formatIntersect(intersected)
print iString

print 'Process: Intersect'
arcpy.Intersect_analysis(iString, 
                          Intersect71Census, "ALL", "", "INPUT")

print 'Process Results'
dictionary = createResultDic(Intersect71Census, dataKey, fields)

print 'Create CSV'
header = [dataKey]
for field in fields:
    header.append(field)
createCSV(header,csvname, 'wb' )

for counter, busStop in enumerate(dictionary.keys()):
    datakeys  = dictionary[busStop]
    averages = [busStop]

    for key in datakeys:
        data = datakeys[key]
        average = sum(data)/len(data)
        averages.append(average)
    createCSV(averages,csvname)

print "Data Analysis Complete"

摘要

在本章中,我们讨论了如何将自动生成的代码进行泛化,同时添加可以在其他脚本中重用的功能,这将使生成必要的代码组件,如 SQL 语句,变得更加容易。我们还讨论了何时最好不过度创建函数,以避免使它们过于特定。

在下一章中,我们将探讨强大的数据访问模块及其搜索游标、更新游标和插入游标。

第五章。ArcPy 游标 – 搜索、插入和更新

现在我们已经了解了如何使用 ArcPy 与 ArcToolbox 工具交互,并且我们也已经介绍了如何使用 Python 创建函数和导入模块,我们对如何使用 Python 提高 GIS 工作流程有了基本的了解。在本章中,我们将介绍数据游标和数据访问模块,这些模块是在 10.1 版本中引入的。这些数据访问游标在 arcgisscripting 模块(ArcPy 的前身)和 ArcPy 的早期版本中使用的游标上有了很大的改进。游标不仅可以像我们所看到的那样搜索数据,还可以使用更新游标更新数据,并可以使用插入游标添加新的数据行。

数据游标用于通过逐行迭代方法访问数据表中的数据记录。这个概念是从关系数据库中借用的,在关系数据库中,数据游标用于从 SQL 表达式返回的表中提取数据。游标用于搜索数据,但也用于更新数据或添加新数据。

当我们讨论使用 ArcPy 游标创建数据搜索时,我们不仅仅是在谈论属性信息。新的数据访问模型游标可以直接与形状字段交互,并且当与 ArcPy 几何对象结合使用时,可以执行地理空间函数,从而取代将数据传递给 ArcToolbox 工具的需要。数据访问游标代表了在 Python 自动化 GIS 领域中最有用的创新。

在本章中,我们将介绍:

  • 使用搜索游标访问属性和空间数据

  • 使用更新游标调整行内的值

  • 使用插入游标向数据集中添加新数据

  • 使用游标和 ArcPy 几何对象类型在内存中执行地理空间分析

数据访问模块

随着 ArcGIS 10.1 的发布,新的数据访问模块 arcpy.da 使得数据交互比之前的数据游标更加容易和快速。通过允许以各种形式(形状对象、X 值、Y 值、质心、面积、长度等)直接访问形状字段,以及多种格式(JavaScript 对象表示法 (JSON)、Keyhole 标记语言 (KML)、已知二进制 (WKB)、已知文本 (WKT)),数据访问模块大大提高了 GIS 分析师提取和控制形状字段数据的能力。

数据访问游标接受一系列必需和可选参数。必需参数包括作为字符串(或表示路径的变量)的特征类路径以及要返回的字段。如果需要所有字段,可以使用星号表示法,并提供一个包含星号作为字符串的字段参数列表([ * ])。如果只需要少数几个字段,请提供这些字段的字符串字段名(例如 [ "NAME", "DATE"])。

其他参数是可选的,但对于搜索和更新游标都非常重要。可以提供下一个以 SQL 表达式形式的where子句;此子句将限制从数据集中返回的行数(如最后一章中的脚本中演示的 SQL 表达式所示)。搜索和更新游标使用的 SQL 表达式不是完整的 SQL 表达式,因为游标的选择或更新命令是由游标的选择自动提供的。此参数只需要 SQL 表达式的where子句。

可以在 ArcPy 空间参考格式中提供下一个空间参考;如果数据格式正确,则这不是必需的,但可以用于在运行时将数据转换到另一个投影。然而,无法指定使用的空间变换。第三个可选参数是一个布尔值(或 True/False),声明数据是否应以展开的点(即单个顶点的列表)或原始几何格式返回。最后一个可选参数是另一个列表,可以用来组织游标返回的数据;此列表将包括 SQL 关键字,如 DISTINCT、ORDER BY 或 GROUP BY。然而,此最终参数仅在处理地理数据库时可用。

让我们看看如何使用arcpy.da.SearchCursor进行形状字段交互。如果我们需要生成一个包含特定路线沿线所有公交车站的电子表格,并包含以 X/Y 格式显示的数据位置,我们可以使用 ArcToolbox 中的添加 XY 工具。然而,这会在我们的数据中添加两个新字段,这并不总是允许的,尤其是在数据存储在具有固定模式的企业地理数据库中时。相反,我们将使用数据访问模块中内置的 SHAPE@XY 令牌来轻松提取数据并将其传递给第四章中的createCSV()函数,复杂的 ArcPy 脚本和泛化函数,以及限制结果只到感兴趣站点的 SQL 表达式:

csvname = "C:\Projects\Output\StationLocations.csv"
headers = 'Bus Line Name','Bus Stop ID', 'X','Y'
createCSV(headers, csvname, 'wb') 
sql = "(NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza') OR (NAME = '71 OB' AND BUS_SIGNAG = '48th Avenue')"
with arcpy.da.SearchCursor(Bus_Stops,['NAME', 'STOPID', 'SHAPE@XY'], sql) as cursor:
 for row in cursor:
 linename = row[0]
 stopid = row[1]
 locationX = row[2][0]
 locationY = row[2][1]
 locationY = row[2][1]
 data = linename, stopid, locationX, locationY
 createCSV(data, csvname)

注意,每行数据都作为元组返回;这是有意义的,因为搜索游标不允许任何数据操作,一旦创建,元组就是不可变的。相比之下,来自更新游标的数据以列表格式返回,因为列表可以更新。两者都可以使用之前显示的索引进行访问。

光标返回的每一行都是一个包含三个对象的元组:公交车站的名称、公交车站 ID,以及最后包含该站点 X/Y 位置的另一个元组。变量row中元组内的对象可以通过索引访问:公交车站名称位于索引 0,ID 位于索引 1,位置元组位于索引 2。

在位置元组中,X 值位于索引 0,Y 值位于索引 1;这使得通过传递一个值来访问位置元组中的数据变得很容易,如下所示:

 locationX = row[2][0]

能够将列表、元组甚至字典添加到另一个列表、元组或字典中是 Python 的一个强大组件,这使得数据访问逻辑性和数据组织变得容易。

然而,前一段代码返回的电子表格有几个问题:位置以要素类的本地投影(在本例中为州平面投影)返回,并且有一些数据行是重复的。如果能提供电子表格中的纬度和经度值并删除重复值,那就更有帮助了。在我们将数据传递给 createCSV() 函数之前,使用可选的空间参考参数和列表对数据进行排序:

spatialReference = arcpy.SpatialReference(4326)
sql = "(NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza') OR (NAME = '71 OB' AND BUS_SIGNAG = '48th Avenue')"
dataList = []
with arcpy.da.SearchCursor(Bus_Stops, ['NAME','STOPID','SHAPE@XY'], sql, spatialReference) as cursor:
 for row in cursor:
 linename = row[0]
 stopid = row[1]
 locationX = row[2][0]
 locationY = row[2][1]
 data = linename, stopid, locationX, locationY
 if data not in dataList:
 dataList.append(data)

csvname = "C:\Projects\Output\StationLocations.csv"
headers = 'Bus Line Name','Bus Stop ID', 'X','Y'
createCSV(headers, csvname, 'wb') 
for data in dataList: 

空间参考是通过传递代表所需投影系统的代码来创建的。在这种情况下,WGS 1984 纬度和经度地理系统的代码是 4326,并将其传递给 arcpy.SpatialReference() 方法来创建一个空间参考对象,该对象可以传递给搜索光标。此外,if 条件用于过滤数据,只接受一个列表作为每个停止点进入名为 dataList 的列表。这个新版本的代码将生成一个包含所需数据的 CSV 文件。这个 CSV 文件可以随后通过 www.convertcsv.com/csv-to-kml.htm 提供的服务转换为 KML,或者甚至更好地使用 Python。使用字符串格式化和循环将数据插入到预构建的 KML 字符串中。

属性字段交互

除了形状字段交互之外,数据访问模块光标提供的另一项改进是,可以使用列表调用要素类中的字段,如前所述。早期数据光标需要使用效率较低的 get value 函数调用,或者需要将字段调用得像函数可用的方法一样。新方法允许通过传递一个星号来调用所有字段,这对于访问之前未检查过的要素类中的字段来说是一个有价值的方法。

更有价值的一项改进是,无需知道数据集是要素类还是形状文件,就能访问唯一标识字段。因为形状文件有一个要素 ID 或 FID,而要素类有一个对象 ID,所以编写脚本工具来访问唯一标识字段更困难。数据访问模块光标允许使用 OID@ 字符串从任何类型的输入中请求唯一标识。这使得了解唯一标识的类型变得无关紧要。

如前所述,其他属性字段是通过列表中的字符串请求的。字段名称必须与字段的真正名称匹配;别名名称不能传递给游标。字段可以在列表中以任何所需的顺序排列,并将按请求的顺序返回。列表中只需要包含所需的字段。

这里是一个请求字段信息的演示:

sql = "OBJECTID = 1"
with arcpy.da.SearchCursor(Bus_Stops,
 ['STOPID','NAME', 'OID@'],
 sql) as cursor:
for row in cursor:

如果字段列表中的字段被调整,结果行中的数据将反映这些调整。此外,游标返回的元组的所有成员都可以通过零索引访问。

更新游标

更新游标用于调整现有数据行中的数据。在计算数据或将空值转换为非空值时,更新变得非常重要。结合特定的 SQL 表达式,可以使用新收集或计算的数据值来针对数据进行更新。

注意,运行包含更新游标的代码将更改或更新其操作的数据。在将代码运行在原始数据上之前,制作数据的副本以测试代码是一个好主意。

所有之前讨论过的数据访问模块搜索游标参数对更新游标同样有效。主要区别在于更新游标返回的数据行是以列表形式返回的。因为列表是可变的,所以可以使用列表值赋值来调整它们。

例如,让我们假设 71 路公交车线路将被重命名为 75 路。这将影响所有往返线路,因此必须包含一个 SQL 表达式来获取与该线路相关的所有数据行。一旦创建数据游标,返回的行必须调整名称,重新添加到列表中,并调用更新游标的updateRow方法。以下是这种情况在代码中的样子:

sql = "NAME LIKE '71%'"
with arcpy.da.UpdateCursor(Bus_Stops, ['NAME'],sql),) as cursor:
 for row in cursor:
 lineName = row[0]
 newName = lineName.replace('71','75')
 row[0] = newName

SQL 表达式将返回所有以71开头的数据行;这包括71 IB71 OB。请注意,SQL 表达式必须用双引号括起来,因为属性值需要用单引号。

对于每一行数据,返回的行中位置零的名称被分配给变量lineName。这个变量,一个字符串,使用replace()方法将字符71替换为字符75。这也可以只是将1替换为5,但我想要明确指出正在替换的内容。

生成新的字符串后,它被分配给变量newName。然后使用列表赋值将此变量添加到游标返回的列表中;这将替换列表中最初占据零位置的原始数据值。一旦行值被分配,它随后被传递到游标的updateRow()方法。此方法接受行并更新特定行的要素类中的值。

更新形状字段

对于每一行,光标返回的列表中包含的所有值都可以进行更新,除了唯一标识符(尽管不会抛出异常,UID 值将不会被更新)。即使是形状字段也可以进行调整,但有一些注意事项。主要注意事项是更新的形状字段必须与原始行的几何类型相同,一个点可以被一个点替换,一条线可以被一条线替换,一个多边形可以被另一个多边形替换。

调整点位置

如果一个公交车站从当前位置沿街道向下移动,则需要使用更新光标来更新它。此操作需要一个 X/Y 格式的新的位置,最好与要素类相同的投影,以避免在空间变换中丢失位置精度。根据访问数据的方法,我们有两种创建新点位置的方法。第一种方法是在使用SHAPE@令牌请求位置数据时使用,需要使用 ArcPy 几何类型,在这种情况下是点类型。ArcPy 几何类型将在下一章中详细讨论。

sql = 'OBJECTID < 5'
with arcpy.da.UpdateCursor(Bus_Stops, [ 'OID@', 'SHAPE@'],sql) as cursor:
 for row in cursor:
 row[1] = arcpy.Point(5999783.78657, 2088532.563956)

通过向 ArcPy 点几何传递 X 和 Y 值,创建一个点形状对象并将其传递给光标,在光标返回的更新列表中。将新位置分配给形状字段,然后使用光标的updateRow()方法允许调整形状字段值到新位置。因为前四个公交车站位于同一位置,它们都被移动到新位置。

第二种方法适用于所有其他形状字段交互形式,包括SHAPE@XYSHAPE@JSONSHAPE@KMLSHAPE@WKTSHAPE@WKB令牌。这些通过将新位置以请求的格式返回给光标并更新列表来更新:

sql = 'OBJECTID < 5'
with arcpy.da.UpdateCursor(Bus_Stops, [ 'OID@', 'SHAPE@XY'],sql) as cursor:
 for row in cursor:
 row[1] =(5999783.786500007, 2088532.5639999956)

这里是使用SHAPE@JSON关键字和数据 JSON 表示的相同代码:

sql = 'OBJECTID < 5'
with arcpy.da.UpdateCursor(Bus_Stops, [ 'OID@', 'SHAPE@JSON'],sql) as cursor:
 for row in cursor:
 print row
 row[1] = u'{"x":5999783.7865000069, "y":2088532.5639999956,
 "spatialReference":{"wkid":102643}}'

只要关键字、数据格式和几何类型匹配,位置就会更新到新坐标。关键字方法在更新点时非常有用,然而,SHAPE@XY关键字不能与线或多边形一起使用,因为返回的位置代表请求的几何体的质心。

使用更新光标删除行

如果我们需要删除一行数据,UpdateCursor有一个deleteRow方法可以用来删除行。请注意,这将完全删除数据行,使其无法恢复。此方法不需要传递任何参数;相反,它将删除当前行:

sql = 'OBJECTID < 2'
Bus_Stops = r'C:\Projects\PacktDB.gdb\Bus_Stops'
with arcpy.da.UpdateCursor(Bus_Stops,
 ['OID@',
 'SHAPE@XY'],sql) as cursor:
 for row in cursor:

使用插入光标

现在我们已经掌握了如何更新现有数据的方法,让我们来研究如何使用插入游标创建新数据并将其添加到要素类中。涉及的方法与其他数据访问游标的使用非常相似,只是我们不需要创建一个可迭代的游标来提取数据行;相反,我们将创建一个具有特殊 insertRow 方法的游标,该方法能够逐行将数据添加到要素类行中。

可以使用相同的 with..as 语法来调用插入游标,但通常它是在脚本流程中作为一个变量创建的。

注意

注意,一次只能调用一个游标;如果没有首先使用 Python 的 del 关键字删除初始游标来从内存中移除游标变量,创建两个插入(或更新)游标将生成异常(Python 错误)。这就是为什么许多人都更喜欢使用 with..as 语法。

数据访问模块的插入游标需要与其他游标相同的某些参数。需要写入的要素类和将要插入数据的字段列表(包括形状字段)是必需的。空间参考将不会被使用,因为新的形状数据必须与要素类具有相同的空间参考。不允许在插入游标中使用 SQL 表达式。

要添加到要素类中的数据将以元组或列表的形式存在,其顺序与字段列表参数中列出的字段顺序相同。只需要包含感兴趣的字段在字段列表中,这意味着不是每个字段都需要在添加到列表中的值。当向要素类添加新行数据时,唯一 ID 将自动生成,因此不需要在添加字段的列表中显式包含唯一 ID(以 OID@ 关键字的形式)。

让我们探索可以用来生成新的公交车站的代码。我们将写入一个名为 TestBusStops 的测试数据集。我们只对名称和车站 ID 字段感兴趣,因此这些字段以及形状字段(位于州平面投影系统中)将被包含在要添加的数据列表中:

Bus_Stops = r'C:\Projects\PacktDB.gdb\TestBusStops'
insertCursor = arcpy.da.InsertCursor(Bus_Stops, ['SHAPE@', 'NAME','STOPID'])
coordinatePair = (6001672.5869999975, 2091447.0435000062)
newPoint = arcpy.Point(*coordinatePair)
dataList = [newPoint,'NewStop1',112121]
insertCursor.insertRow(dataList)
del insertCursor

如果有一个可迭代的列表数据需要插入到要素类中,在进入迭代之前创建插入游标变量,一旦数据迭代完成,就删除插入游标变量,或者使用 with..as 方法在迭代完成后自动删除插入游标变量:

Bus_Stops = r'C:\Projects\PacktDB.gdb\TestBusStops'
listOfLists = [[(6002672.58675, 2092447.04362),'NewStop2',112122],
 [(6003672.58675, 2093447.04362),'NewStop3',112123],
 [(6004672.58675, 2094447.04362),'NewStop4',112124]
 ]

with arcpy.da.InsertCursor(Bus_Stops,
 ['SHAPE@',
 'NAME',
 'STOPID']) as iCursor:
 for dataList in listOfLists:
 newPoint = arcpy.Point(*dataList[0])
 dataList[0] = newPoint

作为列表,listOfLists 变量是可迭代的。其中的每个列表在迭代时被视为 dataListdataList 中的第一个值(坐标对)被传递给 arcpy.Point() 函数以创建一个 Point 对象。arcpy.Point() 函数需要两个参数,XY;这些参数通过星号从坐标对元组中提取,星号会将元组“展开”,并将包含的值传递给函数。然后,使用基于索引的列表赋值将 Point 对象添加回 dataList,如果 dataList 变量是一个元组(我们则必须创建一个新的列表,并将 Point 对象和其他数据值添加进去),则我们无法使用这种方法。

插入折线几何

要从一系列点创建并插入多边形类型的形状字段,最好使用 SHAPE@ 关键字。我们还将进一步探讨 ArcPy 几何类型,这些内容将在下一章中讨论。当使用 SHAPE@ 关键字时,我们必须使用 ESRI 的空间二进制格式中的数据,并且必须使用 ArcPy 几何类型以相同的格式将数据写回字段。

要创建折线,有一个要求,至少由两个坐标对组成的两个有效点。当使用 SHAPE@ 关键字时,有一种将坐标对转换为 ArcPy 点并将其添加到 ArcPy 数组中的方法,然后该数组被转换为 ArcPy 折线,并写回到形状字段:

listOfPoints = [(6002672.58675, 2092447.04362),
 (6003672.58675, 2093447.04362),
 (6004672.58675, 2094447.04362)
 ]
line = 'New Bus Line'
lineID = 12345
busLine = r'C:\Projects\PacktDB.gdb\TestBusLine'
insertCursor = arcpy.da.InsertCursor(busLine, ['SHAPE@',                   'LINE', 'LINEID'])
lineArray = arcpy.Array()
for pointsPair in listOfPoints:
 newPoint = arcpy.Point(*pointsPair)
 lineArray.add(newPoint)
newLine = arcpy.Polyline(lineArray)
insertData = newLine, line, lineID

元组中的三个坐标对被迭代并转换为 Point 对象,这些对象随后被添加到名为 lineArray 的数组对象中。然后,该数组对象被添加到名为 newLine 的折线对象中,该对象随后与其它数据属性一起添加到一个元组中,并通过 InsertCursor 插入到要素类中。

插入多边形几何

多边形也是通过游标插入或更新的。ArcPy 多边形几何类型不需要坐标对包含第一个点两次(即作为第一个点和最后一个点)。多边形会自动通过 arcpy.Polygon() 函数关闭:

listOfPoints = [(6002672.58675, 2092447.04362),
 (6003672.58675, 2093447.04362),
 (6004672.58675, 2093447.04362),
 (6004672.58675, 2091447.04362)
 ]
polyName = 'New Polygon'
polyID = 54321
blockPoly = r'C:\Projects\PacktDB.gdb\Chapter5Results\TestPolygon'
insertCursor = arcpy.da.InsertCursor(blockPoly, ['SHAPE@', 'BLOCK', 'BLOCKID'])
polyArray = arcpy.Array()
for pointsPair in listOfPoints:
 newPoint = arcpy.Point(*pointsPair)
 polyArray.add(newPoint)
newPoly = arcpy.Polygon(polyArray)
insertData = newPoly, polyName, polyID
insertCursor.insertRow(insertData)

下面是插入操作结果的可视化:

插入多边形几何

摘要

在本章中,我们介绍了数据访问模块游标的基本用法。探讨了搜索、更新和插入游标,并特别关注了这些游标用于从形状字段提取形状数据的使用。还介绍了游标参数,包括空间参考参数和 SQL 表达式 where 子句参数。在下一章中,我们将进一步探讨游标的使用,特别是与 ArcPy 几何类型的使用。

第六章:使用 ArcPy 几何对象

地理空间分析的本质是使用几何形状——点、线和多边形——来模拟现实世界对象及其基于位置的关系。简单的形状及其位置、长度和面积等几何属性通过地理空间操作进行处理,以生成分析结果。正是模型化地理数据和相关的属性信息将地理信息系统与其他所有信息系统区分开来。

在 ArcPy 之前,使用地理空间操作处理要素类几何主要依赖于 ArcToolbox 中的预建工具。ArcPy 使得直接访问存储在要素类形状字段中的数学表示的几何形状成为可能。一旦访问,这些几何数据就被加载到 ArcPy 几何对象中,以便在 ArcPy 脚本中进行分析。正因为这一进步,编写访问几何字段并使用它们进行分析的脚本已经改变了 ArcGIS 的地理空间分析。在本章中,我们将探讨如何生成和使用 ArcPy 几何对象来执行地理空间操作,并将它们应用于公交车站分析。

在本章中,我们将介绍:PointArray构造函数对象以及PointGeometryPolylinePolygon几何对象

  • 如何使用几何对象执行地理空间操作

  • 如何将几何对象集成到脚本中

  • 如何使用几何对象执行常见的地理空间操作

  • 如何在脚本中使用几何对象方法替换 ArcToolbox 工具的使用

ArcPy 几何对象类

在设计几何对象时,ArcPy 的作者使得在内存中执行地理空间操作成为可能,减少了使用 ArcToolbox 中工具进行这些操作的需求。这将导致速度提升,因为在分析的每一步中都不需要将计算结果写入磁盘。相反,步骤的结果可以在脚本中的函数之间传递。分析的最后结果可以写入硬盘上的要素类,或者可以写入电子表格或传递给另一个程序。

几何对象被编写为 Python 类——包含内部函数的特殊代码块。内部函数是几何对象的方法和属性;当被调用时,它们允许对象执行操作(方法)或揭示有关几何对象的信息(属性)。Python 类使用一个包含共享方法和属性的父类编写,以及引用父类但还具有特定方法和不共享属性的子类。在这里,父类是 ArcPy 的Geometry对象,而子类是PointGeometryMultipointPolylinePolygon对象。

几何对象可以通过三种方式生成。第一种方法需要使用数据游标来读取现有的要素类,并传递一个特殊的关键字作为字段名。游标返回的形状数据是一个几何对象。第二种方法是通过将原始坐标传递给构造函数对象(无论是Point还是Array对象)来创建新数据,然后将其传递给几何对象。第三种方法是通过使用 ArcToolbox 中的复制要素工具从要素类中读取数据。

每个几何对象都有允许进行读取访问和写入访问的方法。读取访问方法对于访问构成点、线和多边形的坐标点是重要的。写入访问方法在生成新的数据对象时很重要,这些对象可以进行分析或写入磁盘。

PointGeometryMultipointPolylinePolygon几何对象用于对其各自的几何类型进行分析。通用几何对象可以接受任何几何类型,并在不需要区分几何类型时执行地理空间操作。

将使用两个其他 ArcPy 类在内存中执行地理空间操作:Array对象和Point对象。它们是构造函数对象,因为它们不是从几何类派生的子类,而是用于构建几何对象。Point对象用于从原始坐标创建坐标点。Array对象是一系列坐标点,可以传递给PolylinePolygon对象,因为普通的 Python 列表中的 ArcPy Point对象不能用来生成这些几何对象。

ArcPy Point对象

Point对象是用于生成几何对象的基本构建块。此外,所有几何对象在读取访问方法时都将返回组件坐标作为Point对象。Point对象允许通过其XYZ属性进行简单的几何访问,以及有限数量的地理空间方法,例如containsoverlapswithintouchescrossesequalsdisjoint。让我们使用 IDLE 来探索一些这些方法,使用具有相同坐标的两个Point几何对象:

>>> Point = arcpy.Point(4,5)
>>> point1  = arcpy.Point(4,5)
>>> Point.equals(point1)
True
>>> Point.contains(point1)
True
>>> Point. crosses(point1)
False
>>> Point.overlaps(point1)
False
>>> Point.disjoint(point1)
False
>>> Point.within(point1)
True
>>> point.X, Point.Y
(4.0, 5.0)

在这些示例中,我们可以看到Point对象的一些特性。对于具有相同坐标的两个点,equals方法和disjoint方法的结果符合预期。当两个对象不共享坐标时,disjoint方法将返回True,而equals方法则相反。contains方法将适用于两个Point对象并返回Truecrosses方法和overlaps方法的结果有些令人惊讶,因为两个Point对象在位置上确实重叠,并且可以被认为是交叉的;然而,这些方法并没有返回预期的结果,因为它们并不是用来比较两个点的。

ArcPy Array对象

在我们继续学习PolylinePolygon对象之前,我们需要了解 ArcPy 的Array对象。它是Point对象和需要多个坐标点的几何对象之间的桥梁。Array对象接受Point对象作为参数,而Array对象反过来又作为参数传递给要创建的几何对象。让我们使用Point对象和Array对象一起,更好地理解它们是如何协同工作的。

Array对象类似于 Python 列表,具有extendappendreplace方法,并且还具有独特的addclone方法。add方法将用于单独添加Point对象:

>>> Point = arcpy.Point(4,5)
>>> point1  = arcpy.Point(7,9)
>>> Array = arcpy.Array()
>>> Array.add(point)
>>> Array.add(point1)

extend()方法会一次性添加一个Point对象的列表:

>>> Point = arcpy.Point(4,5)
>>> point1 = arcpy.Point(7,9)
>>> pList = [Point,point1]
>>> Array = arcpy.Array()
>>> Array.extend(pList)

insert方法将在 Array 中的特定索引处放置一个Point对象,而replace方法通过传递索引和一个新的Point对象来替换 Array 中的Point对象:

>>> Point  = arcpy.Point(4,5)
>>> point1  = arcpy.Point(7,9)
>>> point2  = arcpy.Point(11,13)
>>> pList = [Point,point1]
>>> Array = arcpy.Array()
 >>> Array.extend(pList)
>>> Array.replace(1,point2)
>>> point3  = arcpy.Point(17,15)
>>> Array.insert(2,point3)

Array对象加载了Point对象后,可以用来生成其他几何对象。

ArcPy Polyline 对象

Polyline对象是由至少包含两个Point对象的Array对象生成的。如下面的 IDLE 示例所示,一旦生成了Array对象并加载了Point对象,它就可以作为参数传递给Polyline对象:

>>> Point  = arcpy.Point(4,5)
>>> point1  = arcpy.Point(7,9)
>>> pList = [Point,point1]
>>> Array = arcpy.Array()

>>> Array.extend(pList)
>>> pLine = arcpy.Polyline(Array)

现在 Polyline 对象已经创建,我们可以访问其方法。这包括揭示多段线中的构成坐标点和其他相关信息的方法:

>>> pLine.firstPoint
<Point (4.0, 5.0, #, #)>
>>> pLine.lastPoint
<Point (7.0, 9.0, #, #)>
pLine.getPart()
<Array [<Array [<Point (4.0, 5.0, #, #)>, <Point (7.0, 9.0, #, #)>]>]>
>>> pLine.trueCentroid
<Point (5.5, 7.0, #, #)>
>>> pLine.length
5.0
>>> pLine.pointCount
2

这个示例Polyline对象尚未分配空间参考系统,因此长度是无单位的。当一个几何对象确实有空间参考系统时,线性单位和面积单位将以系统的线性单位返回。

Polyline对象也是我们第一个可以调用执行地理空间操作(如缓冲区、距离分析和裁剪)的几何对象:

>>> bufferOfLine = pLine.buffer(10)
>>> bufferOfLine.area
413.93744395
>>> bufferOfLine.contains(pLine)
True
>>> newPoint = arcpy.Point(25,19)
>>> pLine.distanceTo(newPoint)
20.591260281974

Polyline对象的另一个有用方法是positionAlongLine方法。它用于在特定位置返回一个PointGeometry对象,如下文所述。这个位置可以是沿线的数值距离,也可以是百分比(以 0-1 之间的浮点数表示),当使用可选的第二个参数时:

>>> nPoint = pLine.positionAlongLine(3)
>>> nPoint.firstPoint.X, nPoint.firstPoint.Y
(5.8, 7.4)>>> pPoint = pLine.positionAlongLine(.5,True)
 >>> pPoint.firstPoint.X,pPoint.firstPoint.Y
(5.5, 7.0)

Polyline对象还有许多其他方法可用。更多信息请参阅:resources.arcgis.com/en/help/main/10.2/index.html#//018z00000008000000

ArcPy Polygon 对象

要创建一个 Polygon 对象,必须使用 Array 对象加载 Point 对象,然后将它作为参数传递给 Polygon 对象。一旦生成了 Polygon 对象,它所拥有的方法对于执行地理空间操作非常有用。几何对象也可以使用 ArcToolbox 的 CopyFeatures 工具保存到磁盘。以下 IDLE 示例演示了如何通过将 Polygon 对象和原始字符串文件名传递给工具来生成 shapefile

>>> import arcpy
>>> point1 = arcpy.Point(12,16)
>>> point2 = arcpy.Point(14, 18)
>>> point3 = arcpy.Point(11, 20)
>>> Array = arcpy.Array()
>>> Points = [point1,point2,point3]
>>> Array.extend(points)
>>> Polygon = arcpy.Polygon(array)
>>> arcpy.CopyFeatures_management(polygon, r'C:\Projects\Polygon.shp')
<Result 'C:\\Projects\\Polygon.shp'>

Polygon 对象缓冲区

Polygon 对象,就像 Polyline 对象一样,拥有使执行地理空间操作(如缓冲区)变得容易的方法。通过将一个数字传递给缓冲区方法作为参数,可以在内存中生成一个缓冲区。数字的单位由 SpatialReference 系统确定。可以通过提供负缓冲区数字来生成内部缓冲区;生成的缓冲区是 Polygon 对象内部在指定距离处的区域。裁剪、并集、对称差集以及更多操作都作为方法提供,还有在内部或包含操作;只要传递了 SpatialReference 系统对象作为参数,就可以使用 Polygon 对象方法执行投影。以下是一个脚本,它将创建两个具有两个不同的 SpatialReference 系统的 shapefile,每个系统由 EPSG 编码系统中的一个数字代码(2227 和 4326)标识:

import arcpyPoint  = arcpy.Point(6004548.231,2099946.033)
point1  = arcpy.Point(6008673.935,2105522.068)
point2  = arcpy.Point(6003351.355,2100424.783)Array = arcpy.Array()
Array.add(point1)
Array.add(point)
array.add(point2)
Polygon = arcpy.Polygon(array, 2227)
buffPoly = Polygon.buffer(50)
features = [Polygon,buffPoly]
arcpy.CopyFeatures_management(features,
 r'C:\Projects\Polygons.shp')
spatialRef = arcpy.SpatialReference(4326)
polygon4326 = Polygon.projectAs(spatialRef)
arcpy.CopyFeatures_management(polygon4326,
 r'C:\Projects\polygon4326.shp')

下面是第二个 shapefile 在 ArcCatalog 预览窗口中的样子:

Polygon 对象缓冲区

其他 Polygon 对象方法

与 ArcToolbox 中的裁剪工具不同,该工具可以使用另一个多边形裁剪多边形,裁剪方法需要一个范围对象(另一个 ArcPy 类)并且仅限于裁剪区域周围的矩形边界。要从多边形中移除区域,可以使用差集方法,就像在 ArcToolbox 中的裁剪或擦除工具一样:

buffPoly = Polygon.buffer(500)
donutHole =buffPoly.difference(Polygon)
features = [Polygon,donutHole]
arcpy.CopyFeatures_management(features,
                              r"C:\Projects\Polygons2.shp")

下面是缓冲区和差集操作类似甜甜圈孔的结果。具有甜甜圈孔的缓冲区包围了原始的 Polygon 对象:

其他 Polygon 对象方法

ArcPy 几何对象

通用几何对象对于在内存中创建特征类的几何形状的副本非常有用,而无需首先知道特征类包含哪种类型的几何形状。与所有 ArcPy 几何对象一样,它的读取方法包括以多种格式提取数据,如 JSON、WKT 和 WKB。每个几何形状的面积(如果它是多边形)、质心、范围和构成点也都可以使用,如前所述。

下面是使用 CopyFeatures 工具将特征类的几何形状读入内存的示例:

import arcpy
cen2010 = r'C:\Projects\ArcPy.gdb\SanFrancisco\CensusBlocks2010'
blockPolys = arcpy.CopyFeatures_management(cen2010,
 arcpy.Geometry())

变量 blockPolys 是一个包含所有加载到其中的几何形状的 Python 列表;在这种情况下,它是人口普查区块。然后可以遍历该列表以进行分析。

ArcPy PointGeometry 对象

PointGeometry对象在执行与Point对象不可用的相同地理空间操作时非常有用。当使用光标从具有PointGeometry类型的要素类中检索形状数据时,形状数据以PointGeometry对象的形式返回。当不使用光标从要素类中检索数据时,需要Point对象来构建所有其他几何对象,但执行点地理空间操作时使用的是PointGeometry对象。

让我们探索从数据访问模块SearchCursor获取PointGeometry对象,并使用返回的数据行创建缓冲点。在我们的公交车站分析中,这将取代使用 ArcToolbox 缓冲工具在每个车站周围创建 400 英尺缓冲区的需求。以下脚本使用字典收集缓冲区对象,然后使用另一个搜索光标搜索人口普查块。要使用SearchCursor()方法访问形状字段,将SHAPE@令牌作为其中一个字段传递。然后,脚本将遍历公交车站,找到与每个车站相交的所有人口普查块:

# Generate 400 foot buffers around each bus stop
import arcpy,csv
busStops = r"C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops"
censusBlocks2010 = r"C:\Projects\PacktDB.gdb\SanFrancisco\CensusBlocks2010"

sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
dataDic = {}
with arcpy.da.SearchCursor(busStops, ['NAME','STOPID','SHAPE@'], sql) as cursor:
 for row in cursor:
 linename = row[0]
 stopid = row[1]
 shape = row[2]
 dataDic[stopid] = shape.buffer(400), linename

现在数据已经检索,并且使用PointGeometry对象的缓冲方法生成了缓冲区,可以使用迭代和搜索光标将缓冲区与人口普查块几何体进行比较。在此分析中,将使用两种地理空间方法:overlapintersect。重叠方法是一个布尔运算,当比较两个几何体时返回 true 或 false。intersect方法用于获取相交的实际面积以及识别每个块的人口。使用intersect需要两个参数:第二个几何对象和一个整数,表示要返回的几何类型(1 为点,2 为线,4 为多边形)。我们希望返回的多边形相交面积具有相交面积和人口数据:

# Intersect census blocks and bus stop buffers
processedDataDic = {} = {}
for stopid in dataDic.keys():
 values = dataDic[stopid]
 busStopBuffer = values[0]
 linename = values[1]
 blocksIntersected = []
 with arcpy.da.SearchCursor(censusBlocks2010, ['BLOCKID10','POP10','SHAPE@']) as cursor:
for row in cursor:
 block = row[2]
 population = row[1]
 blockid = row[0] 
 if busStopBuffer.overlaps(block) ==True:
 interPoly = busStopBuffer.intersect(block,4)
 data = row[0],row[1],interPoly, block
 blocksIntersected.append(data)
 processedDataDic[stopid] = values, blocksIntersected

此部分脚本遍历块并与缓冲的公交车站相交。现在我们可以识别接触每个车站周围缓冲区的块,并且感兴趣的数据已经收集到字典中,可以对其进行处理,并计算所有被缓冲区接触到的块的总体人口:

# Create an average population for each bus stop
dataList = []
for stopid in processedDataDic.keys():
 allValues = processedDataDic[stopid]
 popValues = []
 blocksIntersected = allValues[1]
 for blocks in blocksIntersected:
 popValues.append(blocks[1])
 averagePop = sum(popValues)/len(popValues)
 busStopLine = allValues[0][1]
 busStopID = stopid
 finalData = busStopLine, busStopID, averagePop
 dataList.append(finalData)

现在数据已经创建并添加到列表中,可以使用我们在第四章中创建的createCSV模块将其输出到电子表格中,复杂的 ArcPy 脚本和泛化函数

# Generate a spreadsheet with the analysis results
def createCSV(data, csvname, mode ='ab'):
 with open(csvname, mode) as csvfile:
 csvwriter = csv.writer(csvfile, delimiter=',')
 csvwriter.writerow(data)

csvname = "C:\Projects\Output\StationPopulations.csv"
headers = 'Bus Line Name','Bus Stop ID', 'Average Population'
createCSV(headers, csvname, 'wb') 
for data in dataList:
 createCSV(data, csvname)

数据已经处理并写入电子表格。我们可以对数据进行的一个额外步骤是使用相交面积为每个缓冲区创建一个成比例的人口值。让我们重新处理数据以包括成比例的面积:

dataList = []
for stopid in processedDataDic.keys():
 allValues = processedDataDic[stopid]
 popValues = []
 blocksIntersected = allValues[1]
 for blocks in blocksIntersected:
 pop = blocks[1]
 totalArea = blocks[-1].area
 interArea = blocks[-2].area
 finalPop = pop * (interArea/totalArea)
 popValues.append(finalPop)
 averagePop = round(sum(popValues)/len(popValues),2)
 busStopLine = allValues[0][1]
 busStopID = stopid
 finalData = busStopLine, busStopID, averagePop
 dataList.append(finalData)

现在,脚本正在充分利用 ArcPy 几何对象的力量,并且脚本完全在内存中运行,这避免了产生任何中间数据集。

摘要

在本章中,我们详细讨论了 ArcPy 几何对象的使用。这些多样化的对象具有相似的方法,实际上是从同一个 Python 类派生出来的子类。它们对于执行内存中的地理空间分析非常有用,这避免了从硬盘读取和写入数据,同时也跳过了创建任何中间数据。

ArcPy 几何对象将成为自动化地理空间工作流程的重要组成部分。将它们与搜索游标结合使用,使得 ArcPy 比任何早期的 ArcGIS Python 脚本工具实现都更有用。接下来,我们将原始脚本转换为可以直接从 ArcToolbox 或地理数据库中的个人工具箱中执行的脚本工具。

第七章:创建脚本工具

现在我们已经涵盖了创建和执行 ArcPy 脚本的基础,我们需要采取下一步,创建可重用的脚本工具。创建脚本工具将允许代码的重用性更高,并使为其他 GIS 分析师和客户创建自定义工具变得容易。通过 Python 脚本后端或代码和熟悉的 ArcGIS 工具前端或用户界面,代码的细节对用户来说是隐藏的;它变成了另一个工具,尽管是一个可以节省数天和数周工作的工具。

本章将涵盖以下主题:

  • 向脚本添加参数以接受输入并按用户需求生成输出

  • 创建自定义工具前端和自定义工具箱

  • 将工具前端参数设置以允许其传递参数到代码后端

向脚本添加动态参数

我们在前几章中生成的脚本都使用了硬编码的输入。输入值以字符串或数字的形式写入脚本中,并分配给变量。虽然它们可以被手动更新以替换输入和输出文件路径以及 SQL 语句,但程序员应致力于创建每次使用时无需编辑的脚本。相反,脚本应设计为动态的,并接受文件路径和其他输入作为参数或参数,这与我们创建的函数接受参数的方式非常相似。

Python 的设计考虑了这一点,sys模块中有一个名为sys.argv的方法,它接受在脚本执行时传递给脚本的输入。虽然 ArcPy 及其前身arcgisscripting 模块的设计者在最初利用了sys.argv方法,但随着时间的推移,他们为接受脚本参数设计了 ArcPy 方法。由于在编写 ArcPy 脚本时可以使用这两种方法,并且两者都可在网络上的示例脚本中找到,因此识别sys.argv方法和arcpy.GetParameterAsText()之间的细微差别很重要。这两种方法的主要区别在于sys.argv将动态参数作为一个列表接受。列表的成员通过列表索引访问并分配给变量。Arcpy.GetParameterAsText()是一个接受索引号参数的函数。传递的索引号反映了参数在工具前端中的顺序;第一个参数是零,下一个是一,依此类推。

注意

如果在工具前端调整了参数的顺序,这种调整必须在代码后端得到反映。

使用 arcpy.AddMessage 显示脚本消息

从脚本中获取反馈以评估脚本在执行分析时的进度非常重要。尽管这看起来很简单,但 Python 脚本和编程语言通常默认不提供任何反馈,除了错误和脚本的终止。这对新手程序员来说可能有点令人沮丧,因为所有内置的反馈都是负面的。

为了缓解这种缺乏反馈的情况,使用打印语句可以让脚本在运行时报告分析进度。然而,当使用脚本工具时,print 语句没有任何效果。它们不会在任何地方显示,并且被 Python 可执行文件忽略。要在执行脚本工具时在脚本控制台中显示消息,ArcPy 有一个 arcpy.AddMessage() 方法。

Arcpy.AddMessage 语句被添加到脚本中,以便程序员需要反馈的地方。所需的反馈作为参数传递给方法并显示;无论是列表、字符串、浮点数还是整数。Arcpy.AddMessage 使得检查分析计算的结果变得容易,以确保使用了正确的输入并产生了正确的输出。由于脚本中的这种反馈可以是强大的调试工具,因此每当需要从脚本工具获取反馈时,都应使用 arcpy.AddMessage

注意

注意,传递给 arcpy.AddMessage 的语句仅在脚本作为脚本工具运行时才会显示。

向脚本添加动态组件

要开始将脚本转换为脚本工具,我们应首先将我们在 第六章 中创建的脚本 使用 ArcPy 几何对象 Chapter6_1.py 复制到名为 Chapter7 的新文件夹中,并将其重命名为 Chapter7_1.py。然后,我们可以开始使用 arcpy.GetParameterAsText 替换硬编码的变量为动态变量。还有一个名为 GetParameter 的 ArcPy 方法,它接受对象作为输入,但就我们的目的而言,GetParameterAsText 是要使用的方法。

通过将 arcpy.GetParameterAsTextarcpy.AddMessage 添加到脚本中,我们将迈出创建脚本工具的第一步。必须小心确保从动态参数创建的变量顺序正确,因为重新排序可能会很耗时。一旦将参数添加到脚本中,并将脚本的硬编码部分替换为变量,脚本就准备好成为脚本工具了。

首先,将所有硬编码的变量移动到脚本的顶部。然后,将所有分配的值替换为 arcpy.GetParameterAsText,使它们成为动态值。每个参数都使用零基索引进行引用;然而,它们是单独传递给函数的,而不是作为列表的成员:

#Chapter 7.py
import arcpy, csv
busStops = arcpy.GetParameterAsText(0)
censusBlocks2010 = arcpy.GetParameterAsText(1)
censusBlockField = arcpy.GetParameterAsText(2)
csvname = arcpy.GetParameterAsText(3)
headers = arcpy.GetParameterAsText(4).split(',')
sql = arcpy.GetParameterAsText(5)
keyfields = arcpy.GetParameterAsText(6).split(';')
dataDic = {}

censusFields = ['BLOCKID10',censusBlockField, 'SHAPE@']
if "SHAPE@" not in keyfields:
 keyfields.append("SHAPE@")

arcpy.AddMessage(busStops)
arcpy.AddMessage(censusBlocks2010)
arcpy.AddMessage(censusBlockField)
arcpy.AddMessage(csvname)
arcpy.AddMessage(sql)
arcpy.AddMessage(keyfields)

如您从变量keyfields和变量标题中看到的那样,必须对某些变量应用进一步的处理,因为并非所有变量都将作为字符串使用。在这种情况下,通过使用string函数 split 和在每个分号上拆分传递给变量keyfields的字符串来创建一个列表,而headers变量是通过逗号拆分创建的。对于其他变量,如censusBlockField变量和keyfields变量,添加了SHAPE@关键字,因为每次运行分析时都需要它。如果需要为每次数据运行指定特定的字段,例如BLOCKID10字段,则它可以保留为硬编码在脚本中,或者可以选择成为脚本工具中可选择的字段参数。

然后将变量添加到脚本的正确位置,使脚本准备好成为地理数据库或 ArcToolbox 中自定义工具箱的一部分的脚本工具。然而,我们必须首先创建脚本工具的工具部分,以便收集值并将其传递到脚本中。

创建脚本工具

创建脚本工具是 ArcPy 的强大功能和 ArcToolbox 中工具易用性的强大组合。

第一步是创建一个自定义工具箱来存放脚本工具。为了实现这一点,请执行以下操作:

  1. 打开ArcCatalog,在SanFrancisco.gdb文件地理数据库上右键单击。

  2. 从菜单中选择新建,然后选择工具箱

  3. 将工具箱命名为Chapter8Tools

  4. Chapter8Tools上右键单击,选择添加,然后选择脚本

以下菜单将出现,允许您设置脚本工具。添加一个无空格的标题和标签,以及描述。我更喜欢在前台运行脚本工具以查看其传递的消息,但这不是必需的,并且在需要启动脚本并继续其他任务时可能会令人烦恼。一旦菜单填写完毕,点击下一步

创建脚本工具

下一菜单包含一个输入字段和一个文件对话框按钮,允许用户找到将传递收集到的参数的脚本。使用文件对话框导航到并选择脚本,并确保已勾选在进程中运行 Python 脚本

创建脚本工具

现在,一旦脚本被识别,点击下一步

创建脚本工具

标记和定义参数

下一个对话框是最重要的一个。在这里,将要传递的参数被标记并且定义了它们的数据类型。必须仔细选择每个参数的正确数据类型,因为某些参数可以提供多种数据类型。此外,每个参数的属性也将被定义;这些信息将表征要收集的数据,并有助于确保数据以正确的格式以及正确的数据类型存在。

首先添加要收集的每个参数的显示名。显示名应该解释所需的输入类型。例如,第一个参数将是公交车站的要素类,所以显示名可以是公交车站要素类

注意

确保按照它们在脚本中传递给变量的顺序添加显示名。

标记和定义参数

添加数据类型

接下来,为每个参数添加数据类型。这很复杂,因为将有一个很长的数据类型列表可供选择,并且通常有几个类型可以适用于许多参数。例如,如果创建了一个 shapefile 参数,它将允许用户选择一个 shapefile,正如预期的那样。然而,使用要素类数据类型可能更好,因为这样就可以在分析中使用 shapefiles 和要素类。

将公交车站要素类添加为参数

第一个参数是公交车站要素类,它应该是一个要素类数据类型。点击显示名旁边的数据类型字段,将出现一个下拉列表。一旦选择了数据类型,就可以查看参数列表下方的参数属性。对于公交车站要素类,默认值是可以接受的,因为要素类是必需的,不是多值,没有默认或环境设置,并且不是从任何其他参数获取的。

将公交车站要素类添加为参数

注意

一些参数将需要先选择另一个参数,因为它们是从第一个参数派生出来的值。人口普查区块字段参数和 SQL 语句参数分别从人口普查区块要素类和公交车站要素类参数派生值。

将人口普查区块要素类添加为参数

人口普查区块要素类与公交车站要素类相似。它将是一个要素类数据类型,允许选择 shapefiles 和要素类,并且不需要调整默认参数属性。一旦设置了数据类型,人口普查区块参数就准备好使用了。

将人口普查区块字段添加为参数

人口普查区块字段参数有一个新的变化;它从人口普查区块要素类参数中获取,并且只有在创建第一个参数后才会填充。为了实现这一点,必须设置获取自参数属性。选择字段作为数据类型,然后点击获取自参数属性旁边的空白区域,并从下拉列表中选择Census_Block_Feature_Class。这将创建一个包含在人口普查区块要素类中的字段列表。

将人口普查区块字段作为参数添加

将输出电子表格作为参数添加

由于脚本工具通过分析运行产生的电子表格是一个逗号分隔值CSV)文件,因此选择文本文件作为数据类型。将默认参数属性设置为文件名可以节省时间,并将所需的文件扩展名更容易识别。此外,由于电子表格将由脚本工具作为输出生成,因此方向参数属性应该是输出

将输出电子表格作为参数添加

将电子表格字段名称作为参数添加

作为输出电子表格的表头选择的字段名称应该是字符串数据类型,字段名称之间用逗号分隔,没有空格。脚本使用字符串数据类型的split方法来分隔字段名称。将逗号传递给split方法通过在逗号上分割输入字符串来分隔参数,从而创建一个字段名称列表。该字段名称列表将被csv模块用作创建电子表格的表头。

将电子表格字段名称作为参数添加

将 SQL 语句作为参数添加

SQL 语句参数将需要有用的 SQL 表达式构建器菜单,因此应该是SQL 表达式数据类型。SQL 表达式构建器将使用从公交车站要素类获取的字段。通过点击该属性并从下拉列表中选择Bus_Stop_Feature_Class来将Obtained from parameter属性设置为公交车站要素类。

将 SQL 语句作为参数添加

将公交车站字段作为参数添加

最后一个参数是公交车站字段参数,它是一个字段数据类型,将从公交车站要素类中获取。将多值参数属性从更改为,以允许用户选择多个字段。同时,请记住将获取自参数属性设置为Bus_Stop_Feature_Class,以便字段从公交车站要素类参数中填充。

将公交车站字段作为参数添加

现在所有参数都已描述,并且它们的属性已设置,脚本工具就绪。点击完成以关闭菜单。

检查最终脚本

一旦收集了所有参数,这些变量随后被用来替换硬编码的字段列表或其他静态脚本元素,用从脚本工具收集的新动态参数替换:这样,脚本就变成了一种宝贵的工具,可以用于多次数据分析,因为现在要分析的字段是动态的:

import arcpy, csv
busStops = arcpy.GetParameterAsText(0)
censusBlocks2010 = arcpy.GetParameterAsText(1)
censusBlockField = arcpy.GetParameterAsText(2)
csvname = arcpy.GetParameterAsText(3)
headers = arcpy.GetParameterAsText(4).split(',')
sql = arcpy.GetParameterAsText(5)
keyfields = arcpy.GetParameterAsText(6).split(';')
dataDic = {}
censusFields = [ 'BLOCKID10',censusBlockField,'SHAPE@']
if "SHAPE@" not in keyfields:
 keyfields.append("SHAPE@")

arcpy.AddMessage(busStops)
arcpy.AddMessage(censusBlocks2010)
arcpy.AddMessage(censusBlockField)
arcpy.AddMessage(csvname)
arcpy.AddMessage(sql)
arcpy.AddMessage(keyfields)
x = 0
with arcpy.da.SearchCursor(busStops, keyfields, sql) as cursor:
 for row in cursor:
 stopid = x
 shape = row[-1]
 dataDic[stopid] = []
 dataDic[stopid].append(shape.buffer(400))
 dataDic[stopid].extend(row[:-1])
 x+=1

processedDataDic = {}
for stopid in dataDic.keys():
 values = dataDic[stopid]
 busStopBuffer = values[0]
 blocksIntersected = []
 with arcpy.da.SearchCursor(censusBlocks2010, censusFields) as cursor:
 for row in cursor:
 block = row[-1]
 population = row[1]
 blockid = row[0] 

 if busStopBuffer.overlaps(block) ==True:
 interPoly = busStopBuffer.intersect(block,4)
 data = population,interPoly, block
 blocksIntersected.append(data)
 processedDataDic[stopid] = values, blocksIntersected

dataList = []
for stopid in processedDataDic.keys():
 allValues = processedDataDic[stopid]
 popValues = []
 blocksIntersected = allValues[-1]
 for blocks in blocksIntersected:
 pop = blocks[0]
 totalArea = blocks[-1].area
 interArea = blocks[-2].area
 finalPop = pop * (interArea/totalArea)
 popValues.append(finalPop)
 averagePop = round(sum(popValues)/len(popValues),2)
 busStopLine = allValues[0][1]
 busStopID = stopid
 finalData = busStopLine, busStopID, averagePop
 dataList.append(finalData)

def createCSV(data, csvname, mode ='ab'):
 with open(csvname, mode) as csvfile:
 csvwriter = csv.writer(csvfile, delimiter=',')
 csvwriter.writerow(data)

headers.insert(0,"ID")
createCSV(headers, csvname, 'wb') 
for data in dataList:
 createCSV(data, csvname)

变量x被添加来跟踪字典dataDic的成员,在第六章使用 ArcPy 几何对象的脚本中,它依赖于STOPID字段。为了消除这种依赖性,引入了x

运行脚本工具

现在前端已经设计好,可以接受用户的参数,后端脚本也准备好接受来自前端的数据,工具现在可以执行了。在工具箱中双击脚本工具以打开工具对话框。

运行脚本工具

选择参数,就像使用任何 ArcToolbox 工具一样(例如,使用文件对话框导航文件树到公交车站特征类)。一旦设置了参数,点击确定以执行脚本。

一个可选的调整是在计算平均人口的地方添加一个arcpy.AddMessage行。通过这样做,可以显示单个区块的人口,并且脚本控制台会提供关于脚本进度的反馈。

在脚本中变量finalData定义的下一行插入:

arcpy.AddMessage(finalData)

这一行提供的反馈将清楚地表明脚本正在运行,这在脚本执行长时间分析时很有用。在执行长时间分析时,向用户提供反馈是一个好习惯,这样他们就可以看到脚本按预期工作。当有大量数据传递给arcpy.AddMessage时,将换行符(\n)作为参数传递。这将把数据分成离散的块,使其更容易阅读。

摘要

在本章中,我们学习了如何将脚本转换成一个永久性和可共享的脚本工具,它可以被没有任何编程经验的 ArcGIS 用户使用。通过使用熟悉的 ArcGIS 工具界面创建前端,然后将参数传递给使用 ArcPy 构建的自定义工具,GIS 程序员可以将 ArcToolbox 的易用性与 Python 的强大功能结合起来。

在下一章中,我们将探讨如何使用 ArcPy 控制从地图文档导出地图。通过调整地图元素,如标题和图例,我们可以创建动态的地图输出,以反映地图分析产生的数据特性。在第九章更多 Arcpy 映射技术中,我们将把地图的输出添加到本章创建的脚本工具中。

第八章。ArcPy.Mapping 简介

制作地图是一种艺术,需要通过多年的专注研究地图学才能掌握。信息的视觉展示既令人兴奋又具有挑战性,可以成为地理空间专业人士日常工作流程中的一项有价值的部分。一旦掌握了基础知识,地图学输出就变成了一个持续不断的战斗,以更快的速度制作出更多的地图。ArcPy 再次提供了一个强大的解决方案:arcpy.mapping 模块。

通过允许自动操作所有地图组件,包括地图窗口、图层、图例以及所有文本元素,arcpy.mapping 使得创建、修改和输出多张地图变得快速简单。地图集创建——对于地理空间专业人士来说的另一项重要技能,也通过该模块变得容易。在本章中,我们将讨论通过 arcpy.mapping 可用的基本功能,并使用它来输出公交车站及其周边人口普查区的地图集。

本章将涵盖以下主题:

  • 检查和更新地图文档(MXD)图层数据源

  • 将 MXDs 导出为 PDF 或其他图像格式

  • 调整地图文档元素

使用 ArcPy 与地图文档

认识到先前 arcgisscripting 模块的局限性,ESRI 设计了 ArcPy 模块,不仅能够处理数据,还包含了 arcpy.mapping 模块,以便直接与地图文档(MXDs)及其包含的图层进行交互。这个新模块为地图自动化提供了多种可能性。一个脚本可能有助于识别损坏的图层链接,更新这些图层的数据源,并应用新的配色方案到图层上。另一个脚本可能使用地图模板创建一系列地图,每个地图对应于一个要素数据集中的要素类。第三个脚本可以创建一个地图集,从 MXD 中移动到网格图层中的每个单元格,以输出书籍的页面,甚至可以实时计算坐标。基于最新分析数据动态创建的地图可以在数据生成的同时输出。Arcpy.mapping 将 ArcPy 模块从有用的工具转变为任何地理空间工作流程中的工具。

为了调查 arcpy.mapping 模块的功效,我们需要一个 MXD 模板的帮助。我已经准备了一个包含数据和 MXD 的地图包,我们将用它来完成本章的练习。它包括我们旧金山公交车站分析的数据,我们将继续扩展它以包括地图。

检查和替换图层源

arcpy.mapping 模块的第一种也是最重要的用途是识别和修复地图文档中图层与它们数据源之间的损坏链接。图层符号和 GIS 数据存储是分开的,这意味着图层数据源经常被移动。Arcpy.mapping 提供了一个快速的解决方案,尽管并不完美。

此解决方案依赖于 arcpy.mapping 模块中包含的多个方法。首先,我们需要识别损坏的链接,然后修复它们。为了识别损坏的链接,我们将使用 arcpy.mapping 中包含的 ListBrokenDataSources() 方法。

ListBrokenDataSources() 方法需要将 MXD 路径传递给 arcpy.mapping.MapDocument() 方法。一旦创建地图文档对象,它就会被传递给 ListBrokenDataSources() 方法,并生成一个包含图层对象列表,每个损坏链接的图层都有一个图层对象。图层对象有多个属性可供使用。使用这些属性,让我们使用每个对象的名称和数据源属性打印出每个图层的名称和数据源:

import arcpy
mxdPath = 'C:\Projects\MXDs\Chapter8\BrokenLinks.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
brokenLinks = arcpy.mapping.ListBrokenDataSources(mxdObject)
for link in brokenLinks:
 print link.name, link.dataSource

修复损坏的链接

现在我们已经识别出损坏的链接,下一步是修复它们。在这种情况下,发现数据源应该在名为 Data 的文件夹中,但它们并不包含在该文件夹内。然后脚本必须升级以替换每个图层的数据源,使它们指向实际的数据源位置。

图层对象和地图文档对象中包含的方法可以完成此下一步。如果 MXD 的所有数据源都已移动,则最好使用 MXD 对象及其方法来修复数据源。在示例 MXD 中,数据源已全部移动到一个名为 NewData 的新文件夹中,因此我们将使用 findAndReplaceWorkspacePaths() 方法来修复链接:

oldPath = r'C:\Projects\MXDs\Data'
newPath = r'C:\Projects'
mxdObject.findAndReplaceWorkspacePaths(oldPath,newPath)
mxdObject.save() 

只要数据源仍然以相同的格式存在(例如,形状文件仍然是形状文件或要素类仍然是要素类),findAndReplaceWorkspacePaths() 方法就会起作用。如果数据源类型已更改(例如,形状文件已导入到文件地理数据库中),则必须使用 replaceWorkspaces() 方法代替,因为它需要一个工作空间类型作为参数:

oldPath = r'C:\Projects\MXDs\Data'
oldType = 'SHAPEFILE_WORKSPACE'
newPath = r'C:\Projects'
newType = 'FILEGDB_WORKSPACE'
mxdObject.replaceWorkspaces(oldPath,oldType,newPath,newType)
mxdObject.save()

修复单个图层的链接

如果各个图层不共享数据源,则需要使用图层对象可用的 findAndReplaceWorkspacePath() 方法进行调整。此方法与之前使用的方法类似,但它只会替换应用到的图层对象的数据源,而不是所有图层。当与字典结合使用时,可以通过图层名称属性更新图层数据源:

import arcpy
layerDic = {'Bus_Stops':[r'C:\Projects\OldDataPath', r'C:\Projects'],
 'stclines_streets': [r'C:\Projects\OtherPath', r'C:\Projects']}
mxdPath = r'C:\Projects\MXDs\Chapter8\BrokenLinks.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
brokenLinks = arcpy.mapping.ListBrokenDataSources(mxdObject)
for layer in brokenLinks:
 oldPath, newPath = layerDic[layer.name]
 layer.findAndReplaceWorkspacePath(oldPath, newPath )
 mxdObject.save()

这些解决方案适用于单个地图文档和图层。通过使用内置的 glob 模块的 glob.glob() 方法(它有助于生成匹配特定文件扩展名的文件列表)和 os 模块的 os.path.join() 方法,也可以扩展到包含 MXD 文件的文件夹:

import arcpy, glob, os
oldPath = r'C:\Projects\MXDs\Data'
newPath = r'C:\Projects'
folderPath = r'C:\Projects\MXDs\Chapter8'
mxdPathList = glob.glob(os.path.join(folderPath, '*.mxd'))
for path in mxdPathList: 
 mxdObject = arcpy.mapping.MapDocument(mxdPath)
 mxdObject.findAndReplaceWorkspacePaths(oldPath,newPath)
 mxdObject.save()

从 MXD 导出为 PDF

arcpy.mapping 的下一个最重要的用途是自动导出 MXDs。以下代码将突出显示 PDF 的导出,但请注意,该模块还支持 JPEG 和其他图像格式的导出。使用 arcpy.mapping 进行此过程是一种乐趣,因为通常打开和导出 MXDs 的过程涉及很多等待 ArcMap 启动和地图加载,这可能会浪费很多时间:

import arcpy, glob, os
mxdFolder = r'C:\Projects\MXDs\Chapter8'
pdfFolder = r'C:\Projects\PDFs\Chapter8'
mxdPathList = glob.glob(os.path.join(mxdFolder, '*.mxd'))
for mxdPath in mxdPathList:
 mxdObject = arcpy.mapping.MapDocument(mxdPath)
 arcpy.mapping.ExportToPDF(mxdObject,
 os.path.join(pdfFolder,
 basepath( 
 mxdPath.replace('mxd','pdf')
 )))

注意

注意,输出文件夹必须存在,此代码才能正确运行。虽然 os 模块有方法来检查路径是否存在(os.path.exists)和创建文件夹(os.mkdir),但这些方法不包括在此代码片段中,并且如果输入或输出路径不存在,arcpy.mapping.ExportToPDF() 方法将抛出异常。

此示例代码非常有用,可以将其转换为接受文件夹路径作为参数的函数。然后,可以将该函数添加到脚本工具中,正如上一章所讨论的那样。

调整地图文档元素

Arcpy.mapping 包含了一些重要的方法,这些方法将有助于地图文档操作的自动化。这些方法包括在 MXDs 中添加新图层或开启/关闭图层的能力,调整数据帧的比例或移动数据帧以聚焦特定区域的能力,以及调整地图文本组件(如标题或副标题)的能力。这些方法将在我们继续进行公交车站分析时进行讨论。

打开名为 MapAdjust.mxd 的 MXD。这代表我们的基础地图文档,其中包含我们将调整以满足我们需求的图层和元素。它包含我们从分析中生成的图层和填充地图的基础图层。还有许多文本元素,脚本将自动替换以适应每个地图的具体需求。然而,它并不能很好地表示分析结果,因为与公交车站缓冲区相交的普查区重叠,使得地图解读变得困难。

该脚本将替换人口普查区图层和公交车站图层的源数据,以便只为每个公交车站制作一张地图,以及与每个公交车站周围缓冲区相交的普查区。

调整地图文档元素

要实现这一点,我们需要创建两个空的特征类:一个包含所有人口普查区的属性,另一个包含公交车站的属性。这将允许数据源被分析产生的数据所替换。

打开SanFrancisco.gdb文件地理数据库,在Chapter8Results要素数据集上右键单击。从下拉菜单中选择新建,然后选择要素类。将第一个要素类命名为SelectedCensusBlocks并使其成为多边形。在下一个菜单中选择默认关键字,然后在接下来的菜单中,点击导入按钮。从 SanFrancisco 要素数据集中选择CensusBlocks要素类;这将把字段加载到新的要素类中。对第二个名为SelectedBusStops的要素类重复此操作,但确保它是点几何类型,并从BusStops要素类导入模式。对第三个名为SelectedStopBuffers的要素类重复相同的过程,但确保它是点几何类型,并从Buffers要素类导入模式。

一旦创建了要素类,现在就可以使用它们来加载分析的结果。我们将重新在内存中执行分析,并将结果写入新创建的要素类,这样就可以捕捉整个普查区,而不仅仅是与缓冲区相交的部分,这将更好地说明分析的结果。

MapAdjust.mxd地图文档的初始状态包含我们现在熟悉的几个要素类:下载的要素类Bus_Stops、生成的要素类Buffers、相交和裁剪的Census Blocks,以及用于制图的四个要素类,即街道要素类、公园要素类、一个邻里要素类,以及旧金山的轮廓。有两个数据框,一个默认命名为图层,另一个命名为插入框,用于创建将显示图层数据框在旧金山周围移动位置的小插入框。描述图层数据框范围的矩形是一个在插入框数据框属性中创建的范围框架。

下面是地图文档初始状态的导出视图:

调整地图文档元素

这里提出的想法是,利用我们分析初始结果来生成人口层、公交车站层和缓冲层的符号。一旦设置并保存,它们就可以作为我们将从这个基本地图文档中生成的单个地图页面的基础。

注意

注意构成标题和副标题的文本元素,以及右窗格底部的图例和归属文本。这些元素可以通过使用arcpy.mapping.ListElements()方法与地图文档中构成图层和数据源的其他元素一起进行调整。

自动地图文档调整

现在我们已经了解了地图文档的初始配置,我们将介绍一个将自动调整的脚本。此脚本将包括我们在本章和前几章中介绍的一些概念,并将介绍一些用于地图文档调整的新方法,我们将在以下内容中详细说明:

import arcpy, os
dirpath = os.path.dirname
basepath = os.path.basename
Bus_Stops = r"C:\Projects\SanFrancisco.gdb\Bus_Stops"
selectedBusStop = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedBusStop'
selectedStopBuffer = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedStopBuffer'
CensusBlocks2010 = r"C:\Projects\SanFrancisco.gdb\CensusBlocks2010"
selectedBlock = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedCensusData'
pdfFolder = r'C:\Projects\PDFs\Chapter8\Map_{0}'
bufferDist = 400
sql = "(NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza')"
mxdObject = arcpy.mapping.MapDocument("CURRENT")
dataFrame = arcpy.mapping.ListDataFrames(mxdObject, "Layers")[0]
elements = arcpy.mapping.ListLayoutElements(mxdObject)
for el in elements:
 if el.type =="TEXT_ELEMENT":
 if el.text == 'Title Element':
 titleText = el
 elif el.text == 'Subtitle Element':
 subTitleText = el
arcpy.MakeFeatureLayer_management(CensusBlocks2010, 'blocks_lyr') 
layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
layerStops = layersList[0]
layerCensus = layersList[1]
layerBuffer = layersList[2]
layerBlocks = layersList[3] 
if layerBlocks.dataSource != selectedBlock:
 layerBlocks.replaceDataSource(dirpath(dirpath(layerBlocks.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedBlock))
if layerStops.dataSource != selectedBusStop:
 layerStops.replaceDataSource(dirpath(dirpath(layerStops.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedBusStop))
if layerBuffer.dataSource != selectedStopBuffer:
 layerBuffer.replaceDataSource(dirpath(dirpath(layerBuffer.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedStopBuffer))
layerStops.visible = True
layerBuffer.visible = True
layerCensus.visible = False
with arcpy.da.SearchCursor(Bus_Stops,['SHAPE@','STOPID','NAME',
 'BUS_SIGNAG' ,'OID@','SHAPE@XY'],sql) as cursor:
 for row in cursor:
 stopPointGeometry = row[0]
 stopBuffer = stopPointGeometry.buffer(bufferDist)
 with arcpy.da.UpdateCursor(layerBlocks,['OID@']) as dcursor:
 for drow in dcursor:
 dcursor.deleteRow()
 arcpy.SelectLayerByLocation_management('blocks_lyr', 'intersect', stopBuffer, "", "NEW_SELECTION")
 with arcpy.da.SearchCursor('blocks_lyr',['SHAPE@','POP10','OID@']) as bcursor:
 inCursor = arcpy.da.InsertCursor(selectedBlock,['SHAPE@','POP10'])
 for drow in bcursor: 
 data = drow[0],drow[1]
 inCursor.insertRow(data)
 del inCursor
 with arcpy.da.UpdateCursor(selectedBusStop,['OID@']) as dcursor:
 for drow in dcursor:
 dcursor.deleteRow()
 inBusStopCursor = arcpy.da.InsertCursor(selectedBusStop,['SHAPE@'])
 data = [row[0]]
 inBusStopCursor.insertRow(data)
 del inBusStopCursor
 with arcpy.da.UpdateCursor(selectedStopBuffer,['OID@']) as dcursor:
 for drow in dcursor:
 dcursor.deleteRow()
 inBufferCursor = arcpy.da.InsertCursor(selectedStopBuffer,['SHAPE@'])
 data = [stopBuffer]
 inBufferCursor.insertRow(data)
 del inBufferCursor
 layerStops.name = "Stop #{0}".format(row[1])
 arcpy.RefreshActiveView()
 dataFrame.extent = arcpy.Extent(row[-1][0]-1200,row[-1][1]-1200,
 row[-1][0]+1200,row[-1][1]-1200)
 subTitleText.text = "Route {0}".format(row[2])
 titleText.text = "Bus Stop {0}".format(row[1])
 outPath  = pdfFolder.format( str(row[1])+ "_" + str(row[-2])) + '.pdf'
 print outPath
 arcpy.mapping.ExportToPDF(mxdObject,outPath)
 titleText.text = 'Title Element'
 subTitleText.text = 'Subtitle Element'
 arcpy.RefreshActiveView()

哇!代码真的很多。让我们逐节回顾,以了解脚本的每一部分都在做什么。

此代码将在 MXD 的 Python 窗口中运行,因此请确保打开 MXD。一旦打开,请打开Python窗口,并在其中右键单击,然后从右键菜单中选择加载。使用文件导航浏览器,找到名为Chapter8_6_AdjustmapCURRENT.py的脚本,并单击它以选择它。点击确定,它将在 Python 窗口中加载。按Enter键将执行脚本,或使用滚动条浏览加载的行。

变量

在脚本中,首先创建了一些变量来保存string文件路径、integer缓冲距离和用于识别感兴趣公交线路的sql语句:

import arcpy, os
Bus_Stops = r"C:\Projects\SanFrancisco.gdb\Bus_Stops"
selectedBusStop = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedBusStop'
selectedStopBuffer = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedStopBuffer'
CensusBlocks2010 = r"C:\Projects\SanFrancisco.gdb\CensusBlocks2010"
selectedBlock = r'C:\Projects\SanFrancisco.gdb\Chapter8Results\SelectedCensusData'
pdfFolder = r'C:\Projects\PDFs\Chapter8\Map_{0}'
bufferDist = 400
sql = "(NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza')"

这些将在以后被用来允许我们搜索图层并对它们进行分析。

地图文档对象和文本元素

因为此代码将在打开的地图文档中执行,所以我们不需要将 MXD 文件路径传递给arcpy.mapping.MapDocument()方法。相反,我们将使用关键字CURRENT来表示我们正在引用打开的地图文档:

mxdObject = arcpy.mapping.MapDocument("CURRENT")
dataFrame = arcpy.mapping.ListDataFrames(mxdObject, "Layers")[0]
elements = arcpy.mapping.ListLayoutElements(mxdObject)
for el in elements:
 if el.type =="TEXT_ELEMENT":
 if el.text == 'Title Element':
 titleText = el
 elif el.text == 'Subtitle Element':
 subTitleText = el

一旦创建了地图文档对象,就使用ListDataFrames()方法从数据帧列表中选择图层数据帧,并将其传递给名为 dataFrame 的变量。

接下来,使用ListLayoutElements()方法将布局元素作为列表传递给 elements 变量。布局元素包括地图文档布局视图的各种元素:图例、整洁线、指向北方的箭头、比例尺和用作标题和描述的文本元素。不幸的是,返回的列表没有很好的顺序,因为它们在布局中的位置是不确定的。要访问我们希望分配给变量以供以后使用的文本元素,必须使用元素对象的两个属性:类型和文本。我们想要调整标题和副标题元素,因此使用for循环遍历元素列表,并使用属性来找到感兴趣的元素。

图层对象

Make Feature Layer 工具是数据管理工具集的一部分,用于将数据从磁盘复制到内存中作为一个图层。ArcGIS 需要生成图层来对数据进行选择和操作,而不是直接在要素类上操作。通过使用图层来执行这些操作,可以保护源要素类。

使用 ArcPy 的MakeFeatureLayer_management()方法访问 Make Feature Layer 工具。当在 Python 窗口中使用此工具时,结果作为将在目录表中可见的图层添加到地图文档中。当在 ArcMap 的 Python 窗口中未使用此工具时,生成的图层仅存在于内存中,并不会添加到地图文档中。

在以下代码的部分中,通过传递人口普查区块要素类的文件路径,在内存中生成了一个名为blocks_lyr的图层。然后使用arcpy.mapping()模块的ListLayers()方法访问初始 MXD 中包含的图层对象。它们按照在地图文档的目录表中列出的顺序返回,并使用列表索引分配给变量,包括新创建的blocks_lyr

arcpy.MakeFeatureLayer_management(CensusBlocks2010, 'blocks_lyr') 
layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
layerStops = layersList[0]
layerCensus = layersList[1]
layerBuffer = layersList[2]
layerBlocks = layersList[3] 

替换数据源

现在我们已经将图层对象分配给了变量,我们将检查它们的数据源是否是我们用于地图生产的正确要素类。使用每个图层对象的dataSource属性,我们将它们与我们要用作数据源的文件路径变量进行比较:

if layerBlocks.dataSource != selectedBlock:
 layerBlocks.replaceDataSource(dirpath(dirpath(layerBlocks.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedBlock))
if layerStops.dataSource != selectedBusStop:
 layerStops.replaceDataSource(dirpath(dirpath (layerStops.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedBusStop))
if layerBuffer.dataSource != selectedStopBuffer:
 layerBuffer.replaceDataSource(dirpath( dirpath(layerBuffer.dataSource)),
 'FILEGDB_WORKSPACE',basepath(selectedStopBuffer))

使用If语句检查数据源是否正确。如果不正确,将使用replaceDataSource()图层方法将它们替换为正确的数据源。此方法需要三个参数:工作空间(在这种情况下,文件地理数据库)、工作空间类型以及新要素类数据源名称,该名称必须位于同一工作空间中,以便replaceDataSource()方法能够工作(尽管它不需要位于同一要素数据集中)。

调整图层可见性

图层对象有一个属性允许我们调整它们的可见性。将此布尔属性设置为TrueFalse将调整图层的可见性(开启为 True,关闭为 False):

layerStops.visible = True
layerBuffer.visible = True
layerCensus.visible = False

我们希望图层变量layerCensus,即新的blocks_lyr对象,被关闭,因此将其设置为False,但公交车站和缓冲区图层对象需要可见,因此将它们设置为True

从公交车站要素类生成缓冲区

所有变量都已生成或分配,所以下一步是使用SearchCursor搜索选定的公交车站。对于每个公交车站,将生成缓冲区对象以找到与这些单独公交车站相交的人口普查区块:

with arcpy.da.SearchCursor(Bus_Stops,['SHAPE@','STOPID','NAME',
 'BUS_SIGNAG' ,'OID@','SHAPE@XY'],sql) as cursor:
 for row in cursor:
 stopPointGeometry = row[0]
 stopBuffer = stopPointGeometry.buffer(bufferDist)
 with arcpy.da.UpdateCursor(layerBlocks,['OID@']) as                   dcursor:
 for drow in dcursor:
 dcursor.deleteRow()

对于从公交车站特征类检索到的每一行数据,返回一系列属性,包含在一个元组中。其中第一个,row[0],是一个PointGeometry对象。此对象有一个缓冲方法,用于在内存中生成一个缓冲Polygon对象,然后将其分配给stopBuffer变量。一旦创建了缓冲区对象,就使用数据访问 UpdateCursor 的deleteRow()方法来擦除人口普查区图层中的行。一旦删除了行,该图层就可以用下一节中将要识别的新选定的人口普查区重新填充。

交集公交车站缓冲区和人口普查区

为了识别与每个公交车站周围的缓冲区相交的人口普查区,使用 ArcToolbox 工具SelectLayerByLocation通过 ArcPy 方法SelectLayerByLocation_management()调用:

arcpy.SelectLayerByLocation_management('blocks_lyr', 'intersect', stopBuffer, "", "NEW_SELECTION")
 with arcpy.da.SearchCursor('blocks_lyr', ['SHAPE@', 'POP10','OID@']) as bcursor:
 inCursor = arcpy.da.InsertCursor(selectedBlock,['SHAPE@', 'POP10'])
 for drow in bcursor: 
 data = drow[0],drow[1]
 inCursor.insertRow(data)
 del inCursor

此方法需要内存中的blocks_lyr图层对象和分配给变量stopBuffer的新创建的缓冲区对象。它还需要选择类型intersect以及另一个参数,该参数控制选择是否添加到现有选择中,或者将是一个新选择。在这种情况下,我们想要一个新选择,因为只需要与当前公交车站相交的人口普查区。

一旦选择了人口普查区并进行了识别,形状数据和人口数据通过InsertCursor传递到由变量selectedBlock表示的特征类。InsertCursor必须使用 del 关键字删除,因为一次只能有一个InsertCursorUpdateCursor在内存中。

填充选定的公交车站和缓冲特征类

以类似的方式,下一步是填充将在地图制作中使用的公交车站和缓冲特征类。首先使用deleteRow()方法将公交车站特征类清空,然后将选定的公交车站形状字段数据插入到特征类中。然后对公交车站缓冲特征类和缓冲几何对象执行相同的步骤:

 with arcpy.da.UpdateCursor(selectedBusStop,['OID@']) as dcursor:
 for drow in dcursor:
 dcursor.deleteRow()
 inBusStopCursor = arcpy.da.InsertCursor(selectedBusStop,['SHAPE@'])
 data = [row[0]]
 inBusStopCursor.insertRow(data)
 del inBusStopCursor
 with arcpy.da.UpdateCursor(selectedStopBuffer,['OID@']) as dcursor:
 for drow in dcursor:
 dcursor.deleteRow()
 inBufferCursor = arcpy.da.InsertCursor(selectedStopBuffer,['SHAPE@'])
 data = [stopBuffer]
 inBufferCursor.insertRow(data)
 del inBufferCursor

更新文本元素

现在数据已经生成并写入创建来存储它们的特征类后,下一步是更新布局元素。这包括将影响图例、数据框架范围和文本元素的图层属性:

layerStops.name = "Stop #{0}".format(row[1])
dataFrame.extent = arcpy.Extent(row[-1][0]-1200,row[-1][1]-1200,
 row[-1][0]+1200,row[-1][1]-1200)
subTitleText.text = "Route {0}".format(row[2])
titleText.text = "Bus Stop {0}".format(row[1])
arcpy.RefreshActiveView()

使用其名称属性调整公交车站图层名称,以反映当前公交车站。通过创建一个arcpy.Extent对象并传递四个参数:XminYminXmaxYmax来调整数据框架范围。为了生成这些值,我使用了 1200 英尺这个相对任意的值,在公交车站周围创建一个正方形。使用它们的文本属性更新文本元素。最后,使用RefreshActiveView()方法确保地图文档窗口正确更新到新的范围。

将调整后的地图导出为 PDF

最后一步是将新调整的地图文档对象传递给 ArcPy 的 ExportToPDF 方法。此方法需要两个参数,即地图文档对象和表示 PDF 文件路径的字符串:

outPath = pdfFolder.format(str(row[1])+"_"+ str(row[-2]))+'.pdf'
arcpy.mapping.ExportToPDF(mxdObject,outPath)
titleText.text = 'Title Element'
subTitleText.text = 'Subtitle Element'
arcpy.RefreshActiveView()

PDF 文件路径字符串是由 pdfFolder 字符串模板和公交站点的 ID 以及对象 ID 和文件扩展名 .pdf 生成的。一旦这些和变量 mxdObject 所表示的地图文档对象传递给 ExportToPDF 方法,就会生成 PDF。然后重置文本元素并刷新视图,以确保地图文档在下次使用脚本时准备就绪。

在 Python 窗口中运行脚本

如果尚未打开,请打开名为 MapAdjust.mxd 的地图文档。打开 Python 窗口,并在窗口中右键单击。从菜单中选择 Load。当文件对话框打开时,找到名为 Chapter8_6_AdjustmapCURRENT.py 的脚本并选择它,确保其中的文件路径正确。点击 OK,它将在 Python 窗口中加载。脚本加载后,按一次 Enter 运行脚本。脚本运行可能需要几秒钟或更长时间才能明显看出它在运行。

注意

注意,在大多数情况下,Python 窗口并不是执行 ArcPy 脚本的好地方,因为它与 IDE 相比有些局限。使用它来加载和执行执行这些地图文档调整的脚本,是 Python 窗口最好的用途之一。

一旦脚本开始运行,地图文档的调整将开始出现并重复。这是一个令人着迷的过程,因为运行脚本的效果以一种在运行 Python 脚本时不太容易看到的方式显现出来。一旦开始生成 PDF 文件,打开一个查看输出。脚本将为所选公交线上的每个公交站生成一个地图,因此生成一定数量的 PDF 文件后,您可以自由关闭地图文档。

这里是输出示例:

在 Python 窗口中运行脚本

脚本生成的地图显示每个公交站位于中心,周围是缓冲区和与缓冲区相交的符号化人口普查区块。标题、副标题和图例已调整,以指示地图中描绘的公交站。使用 ArcPy,我们现在控制着地理空间分析的各个方面:分析本身以及描绘输出结果的制图生产。

摘要

在本章中介绍了 arcpy.mapping 并使用它来控制需要调整以创建自定义地图的地图文档元素。通过将地理空间分析和地图制作结合起来,我们更接近于利用 ArcPy 的全部功能。

在下一章中,我们将进一步探讨 arcpy.mapping,并创建一个可以添加到 ArcToolbox 中的脚本工具,该工具将运行分析并从结果数据生成地图。我们还将完善脚本,并介绍数据驱动页面,以讨论如何在该 ArcPy 脚本中使用这个强大的工具。

第九章。更多 ArcPy 映射技术

能够控制地图文档的制图,同时运行地理空间分析,增加了 ArcPy 的功能和实用性。arcpy.mapping 的属性和方法可以用来操作图层对象、地图比例尺和数据帧范围,甚至可以设置定义查询。通过将自动地理空间分析与动态地图制作相结合,可以实现脚本化映射系统。本章将涵盖以下主题:

  • Arcpy.mapping 图层对象

  • 图层对象定义查询和范围

  • Arcpy.mapping 数据帧对象

  • 创建动态缩放的地图

使用 arcpy.mapping 控制图层对象

Arcpy.mapping 图层对象用于控制地图文档数据帧中图层的属性。通过图层对象属性可以开启和关闭图层可见性、添加新图层以及调整图层顺序。

创建图层对象涉及向 arcpy.mapping.ListLayers() 方法传递参数。如第八章所述,在引用 arcpy.mapping.MapDocument 对象时,可以使用基于零的索引访问地图文档中的图层。此代码将打印出在 MXD 中名为“Layers”的数据帧中包含的图层对象列表:

import arcpy
mxdPath = r'C:\Projects\MXDs\Chapter9\MapDocument1.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
dataFrame = arcpy.mapping.ListDataFrames(mxdObject, "Layers")[0]
layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
print layersList

被称为“Layers”的数据帧中的图层已使用 ListLayers() 方法分配给变量 layersListlayersList 中的每个图层都可以使用基于零的索引进行访问。一旦在列表中访问了图层并将其分配给变量或放入 for 循环中,就可以利用图层对象的属性和方法。

注意

ListLayers 方法的第二个参数在这里为空,但不必如此。它是一个通配符参数,将限制返回的图层对象仅限于与通配符模式匹配的对象。例如,Stops将返回所有以Stops*结尾的图层。可以使用多个星号来查找图层名中包含单词的位置在开头、中间或结尾。

图层对象方法和属性

图层对象属性和方法可以是只读的,这意味着可以检查但不能调整,或者它们是读写,这意味着可以在脚本中进行调整。让我们探索一些这些属性和方法,看看它们如何用来控制从地图文档生成的地图的外观和感觉,以及来自脚本分析的数据。

定义查询

图层对象的一个重要属性是能够动态设置定义查询。定义查询是一个 SQL 语句的where子句,它将可用于显示、查询或其他数据操作(缓冲区、交集等)的数据限制为仅匹配where子句的行。定义查询可以通过在 MXD 中打开图层属性菜单并使用定义查询选项卡来设置,但在这里我们关注的是如何以编程方式添加它们。以下是如何做到这一点的示例:

layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
busStops = layersList[0]
busStops.definitionQuery = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"

这个有价值的属性可以用来重新格式化第八章 ArcPy 映射简介 中的代码。记得Chapter8_6.py脚本的复杂第二部分,其中选择了沿71 Inbound线路的每个公交车站,并将它的几何形状写入另一个要素类?相反,我们可以使用图层对象和定义查询来执行相同类型的几何操作。让我们检查当使用定义查询时,该操作的第一部分(选择公交车站几何形状并在其周围创建缓冲区)看起来是什么样的:

import arcpy
bufferDist = 400
mxdPath = r'C:\Projects\MXDs\Chapter9\MapDocument1.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
dataFrame= arcpy.mapping.ListDataFrames(mxdObject, "Layers")[0]
layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
busStops = layersList[0]
defQuery = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
busStops.definitionQuery = defQuery
idList =[]
with arcpy.da.SearchCursor(busStops,['OID@']) as cursor:
 for row in cursor:
 idList.append(row[0])
for oid in idList:
 newQuery = "OBJECTID = {0}".format(oid)
 print newQuery
 busStops.definitionQuery = newQuery
 with arcpy.da.SearchCursor(busStops,['SHAPE@','STOPID','NAME','BUS_SIGNAG','OID@','SHAPE@XY']) as cursor:
 for row in cursor:
 stopPointGeometry = row[0]
 stopBuffer = stopPointGeometry.buffer(bufferDist)

在这个例子中,定义查询被用来将SearchCursor的潜在结果限制为查询指定的公交车站。然而,这过于繁琐,并且定义查询并没有增加太多价值,因为首先还需要另一个SearchCursor来从busStops图层中提取ObjectID信息。这使代码变得更加复杂,而实际上只需要一个SearchCursor

定义查询应该用来选择与缓冲区相交的块,因为这将消除使用复杂搜索光标和插入光标设置的需要,这种设置在第八章 ArcPy 映射简介 中被采用。让我们重新编写代码,以便在人口普查块图层对象上正确使用定义查询。

第一步是添加一些代码来生成将用作定义查询的 SQL 语句:

import arcpy
bufferDist = 400
mxdPath = r'C:\Projects\MXDs\Chapter9\MapDocument1.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
dataFrame = arcpy.mapping.ListDataFrames(mxdObject, 
 "Layers")[0]
layersList = arcpy.mapping.ListLayers(mxdObject,"",dataFrame)
busStops = layersList[0]
censusBlocks = layersList[3]
sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops,['SHAPE@', 'STOPID', 'NAME', 'BUS_SIGNAG','OID@'],sql) as cursor: 
 for row in cursor:
 bus Query = 'OBJECTID = {0}'.format(row[-1])
 busStops.definitionQuery = bus Query
 stopPointGeometry = row[0]
 stop Buffer = stopPointGeometry. Buffer(bufferDist)
 arcpy.SelectLayerByLocation_management(censusBlocks,'intersect',stopBuffer,"","NEW_SELECTION")
 blockList = []
 with arcpy.da.SearchCursor(censusBlocks,
 ['OID@']) as bcursor:
 for brow in bcursor:
 blockList.append(brow[0])
 newQuery = 'OBJECTID IN ('for COUNTER, oid in enumerate(blockList):
 if COUNTER < len(blockList)-1:
 newQuery += str(oid) + ','
 else:
 newQuery += str(oid)+ ')'
 print newQuery

在本节中,代码将 MXD 中的人口普查块图层分配给变量censusBlocks。然后创建公交车站SearchCursor,并为每一行生成 400 英尺的缓冲区以选择围绕公交车站的普查块。一旦选择了正确的块,就在censusBlocks图层对象上使用第二个SearchCursor来找到所选块的ObjectID(使用OID@令牌)。然后将ObjectIDs附加到名为blockList的列表中。

此列表随后在for循环中进行迭代,以生成一个字符串 SQL 语句。使用分配给变量newQuery的初始字符串,for循环将每个选择区块的ObjectIDs添加到字符串中,以创建一个有效的 SQL 语句。for循环使用enumerate函数来计数for循环执行的循环次数;这允许使用if/then语句。if/then语句确定字符串中ObjectID之后的内容,因为每个ObjectID必须由逗号分隔,除了最后的ObjectID,它必须后跟一个闭括号。for循环生成一个类似于以下示例的 SQL 语句:

OBJECTID IN (910,1664,1812,1813,2725,6382)

结尾的print语句用于展示此段代码的结果,并带来看到代码工作结果的温暖舒适感。一旦我们确信代码正在生成有效的 SQL 语句(闭括号和逗号分隔的ObjectIDs),下一步是将定义查询分配给censusBlocks层对象,并使用结果生成该区域的地图。

控制数据帧窗口的范围和比例

在第八章 ArcPy 映射简介中,我们开始探索数据帧的属性和方法。使用arcpy.Extent对象,我们能够将数据帧的范围设置为脚本中硬编码的范围。然而,这并不总是能够捕捉到大型普查区块的整个范围。通过结合定义查询和数据帧的范围和比例属性,我们可以避免这些不希望的结果。

有两种数据帧对象方法用于将数据帧窗口移动到感兴趣的区域,在这种情况下是选定的普查区块。第一个,我们在这里没有使用,是dataFrame.zoomToSelectedFeatures。第二个是将数据帧的范围属性分配给定义查询后分配给它的普查区块层范围。

我更喜欢第二种方法,因为它即使在没有选定的普查区块的情况下也能工作。此外,由于此脚本生成的地图不应该显示区块的选择,我们不得不添加代码来明确清除已识别的正确普查区块的选择:

 censusBlocks.definitionQuery = newQuery
 dataFrame.extent = censusBlocks.getExtent()
 arcpy.SelectLayerByAttribute_management(censusBlocks,
 "CLEAR_SELECTION")

定义查询使得将数据帧窗口移动到感兴趣的区域变得容易,因为层的范围矩形(或边界)现在仅围绕指定的区块,并且可以将dataFrame的范围属性设置为范围矩形。然而,这并不总是地图学上所希望的,因为看起来将数据帧窗口从范围矩形移回更好。为了做到这一点,我们将访问数据帧对象的缩放属性。

缩放属性可以设置为当前缩放的倍数,以避免在调整数据框架范围时硬编码任何特定距离。当使用缩放属性时,重要的是要记住使用arcpy.RefreshActiveView()方法,因为它将数据框架窗口刷新到新的缩放。

dataFrame.scale = dataFrame.scale * 1.1
arcpy.RefreshActiveView()

由于数据框架范围在之前的几行中已设置,当前缩放表示所选普查区块的范围。要调整它,评估属性并应用一个乘数。在这种情况下,乘数是 1.1,但它可以是任何值。这使得生成的地图看起来更好,因为它为分析结果提供了一些背景上下文。

添加图层对象

在导出地图之前,最后一步是将上面创建的 400 英尺缓冲区作为一个图层添加到数据框架对象中。为了完成这个任务,我们需要提前创建一个符号化图层并复制其符号化,以确保其看起来符合预期。这将作为一个占位符图层添加到MXD中,并分配给脚本中的bufferLayer变量。

  1. 打开一个MXD并添加公交车站要素类

  2. ArcToolbox分析工具集中的邻近工具集中运行缓冲工具,将公交车站要素类作为输入,并将缓冲区大小设置为400 英尺。工具运行完成后,打开缓冲图层属性并按需符号化图层。

  3. 一旦图层被符号化,右键单击图层并选择另存为图层文件

  4. 将图层保存在文件夹中,并关闭MXD

  5. 打开MapDocument1.mxd地图文档,并使用添加数据按钮添加图层。

  6. 确保将其名称更改为400 Foot Buffer,并将其添加到人口部分上面的图例中。

  7. 在脚本中,将缓冲区图层分配给变量bufferLayer

  8. 在脚本较低的部分,在公交车站的SearchCursor中,在生成缓冲区周围的代码下方添加以下行:

    arcpy.CopyFeatures_management(stopBuffer, r"C:\Projects\Output\400Buffer.shp") 
    bufferLayer.replaceDataSource(r"C:\Projects\Output","SHAPEFILE_WORKSPACE","400Buffer")
    
    

这两条线将生成的缓冲区作为 shapefile 复制到磁盘上,然后替换bufferLayer图层对象的源数据。注意,shapefile 的名称不包括.shp扩展名;SHAPEFILE_WORKSPACE参数使得这一点变得不必要。

备注

为了确保每个新的缓冲区 shapefile 可以覆盖现有的 shapefile,在import arcpy行下方添加以下行,以确保文件可以被覆盖:

 arcpy.env.overwriteOutput = 1

导出地图

此脚本的最后一步是导出每个公交车站周围的地图。为此,我们将从Chapter8_6_AdjustMap.py脚本中借用一些代码,并将整个脚本添加到名为Chapter9.py的文件中。此代码将识别并调整标题和副标题元素,使得可以自定义每个生成的 PDF:

import arcpy
arcpy.env.overwriteOutput = 1
bufferDist = 400
pdfFolder = r'C:\Projects\PDFs\Chapter9\Map_{0}'
mxdPath = r'C:\Projects\MXDs\Chapter9\MapDocument1.mxd'
mxdObject = arcpy.mapping.MapDocument(mxdPath)
dataFrame = arcpy.mapping.ListDataFrames(mxdObject,"Layers")[0]
elements = arcpy.mapping.ListLayoutElements(mxdObject)
for el in elements:
 if el.type =="TEXT_ELEMENT":
 if el.text == 'Title Element':
 titleText = el
 elif el.text == 'Subtitle Element':
 subTitleText = el 
layersList = arcpy.mapping.ListLayers(mxdObject,
 "",dataFrame)

busStops = layersList[0]
bufferLayer = layersList[2]
censusBlocks = layersList[4]
sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops,['SHAPE@',
 'STOPID',
 'NAME',
 'BUS_SIGNAG',
 'OID@'],sql) as cursor: 
 for row in cursor:
 busQuery = 'OBJECTID = {0}'.format(row[-1])
 busStops.definitionQuery = busQuery
 stopPointGeometry = row[0]
 stopBuffer = stopPointGeometry.buffer(bufferDist)
 arcpy.CopyFeatures_management(stopBuffer,r"C:\Projects\Output\400Buffer.shp")
 bufferLayer.replaceDataSource(r"C:\Projects\Output",
 "SHAPEFILE_WORKSPACE",
 "400Buffer")
 arcpy.SelectLayerByLocation_management(censusBlocks,
 'intersect',
 stopBuffer,
 "",
 "NEW_SELECTION")
 blockList = []
 with arcpy.da.SearchCursor(censusBlocks,
 ['OID@']) as bcursor:
 for brow in bcursor:
 blockList.append(brow[0])
 newQuery = 'OBJECTID IN ('
 for COUNTER, oid in enumerate(blockList):
 if COUNTER < len(blockList)-1:
 newQuery += str(oid) + ','
 else:
 newQuery += str(oid)+ ')'
 print newQuery
 censusBlocks.definitionQuery = newQuery
 dataFrame.extent = censusBlocks.getExtent()
 arcpy.SelectLayerByAttribute_management(censusBlocks,
 "CLEAR_SELECTION")
 dataFrame.scale = dataFrame.scale * 1.1
 arcpy.RefreshActiveView()
 subTitleText.text = "Route {0}".format(row[2])
 titleText.text = "Bus Stop {0}".format(row[1])
 outPath  = pdfFolder.format( str(row[1])) + '.pdf'
 print outPath
 arcpy.mapping.ExportToPDF(mxdObject,outPath)
 titleText.text = 'Title Element'
 subTitleText.text = 'Subtitle Element'
 censusBlocks.definitionQuery = ''
 busStops.definitionQuery = ''

摘要

在本章中,我们介绍了使用图层定义查询、数据框架范围和比例以及图层源替换来简化地图制作的方法。通过使用定义查询,可以修改图层到新的范围,从而更容易地缩放到图层范围并设置数据框架的比例。定义查询还限制了在数据框架内显示图层的哪些成员。图层源替换被用作制图控制,使我们能够预先生成图层的样式并动态调整它所代表的数据。

在下一章中,我们将结合前三章的教训,使我们能够创建一个脚本工具,该工具可以从分析结果中运行分析并生成电子表格和地图。

第十章。高级几何对象方法

在本章中,我们将讨论之前在第六章中讨论过的高级几何对象方法,即使用 ArcPy 几何对象。本书的目标是介绍 ArcPy 及其模块,同时展示如何在创建持久的 GIS 工作流程时应用这些工具。进行一次分析是好的,但反复进行,只需点击一下按钮,就更好了。使分析结果以行业标准格式共享也是理想的选择。在 ArcGIS 世界中,实现这一目标最佳的方式是使用 ArcPy 和利用几何对象方法的脚本工具。

本章将涵盖以下主题:

  • 将常用函数添加到 Python 路径中的模块

  • 通过添加点生成使分析更高级

  • 高级多边形对象方法

  • 使用 XLWT 创建 Excel 电子表格

创建 Python 模块

创建可重用代码的重要一步是将组件函数打包成一个模块,这样任何脚本都可以从 Python 路径中调用它。首先,我们需要在site-packages文件夹中创建一个文件夹,Python 模块在下载并使用 Python 模块过程提取时放置在这个文件夹中,或者当运行包含在共享模块中的setup.py脚本时。在 Windows 资源管理器中打开site-packages文件夹,导航到C:\Python27\ArcGIS10.2\Lib\site-packages(如果你使用的是标准 Python 2.7 安装,则为C:\Python27\Lib\site-packages)。一旦进入文件夹,创建一个名为common的新文件夹,如图所示:

模块将一个或多个脚本中的函数打包到一个文件夹中,以便与他人共享(尽管它们通常依赖于其他模块才能运行)。我们已经使用了一些内置模块,如csv模块和第三方模块如 ArcPy。让我们探索它们的构建,以了解模块是如何打包以供使用和共享的。

注意

许多模块并没有放置在site-packages文件夹中,但它们需要修改 Python 路径才能使它们可导入。将模块放置在site-packages文件夹中可以消除这一要求。

在 Windows 资源管理器中打开site-packages文件夹,导航到C:\Python27\ArcGIS10.2\Lib\site-packages(如果你使用的是标准 Python 2.7 安装,则为C:\Python27\Lib\site-packages)文件夹。一旦进入文件夹,创建一个名为common的新文件夹,如图所示:

创建 Python 模块

init.py 文件

在这个文件夹中,需要添加一个特殊文件,以便 Python 能够识别该文件夹为模块。这个文件称为__init__.py,它利用 Python 的一个特殊属性,即magic对象或属性,这些属性是内置到 Python 中的。这些magic对象使用前导和尾随的双下划线来避免与自定义函数混淆。

注意

注意,这些是双下划线;单下划线通常用于自定义 Python 类中的所谓私有函数。

__init__.py 文件用于指示该文件夹是一个模块(使其可以通过 import 关键字导入),并通过调用它可能依赖的任何模块来初始化该模块。然而,没有要求在 __init__.py 文件中添加导入命令;它可以是空文件,并且仍然会执行我们所需的模块识别功能。

  1. 打开 IDLEAptana 或您喜欢的 IDE,在名为 common 的文件夹中添加一个新的 Python 文件,并将其命名为 __init__.py。这个文件现在将保持为空。__init__.py 文件

  2. 现在我们已经初始化了模块,我们需要创建一个脚本,用于存放我们的公共函数。让我们称它为 useful.py,因为这些函数将在这个分析和其他分析中非常有用。

  3. 下一步是将我们在早期章节中创建的函数转移过来。这些宝贵的函数被锁定在这些脚本中,所以通过将它们添加到 useful.py,我们将使它们对所有其他我们编写的脚本可用。

    注意

    一个重要的函数是来自 第四章 的 formatSQLMultiple 函数,它使用模板和数据列表生成 SQL 语句。通过将其添加到 useful.py,我们将在需要 SQL 语句时随时可以调用该函数。

  4. 打开脚本 Chapter4Modified2.py,复制函数,然后将其粘贴到 useful.py 中。它没有依赖项,因此不需要修改。

那个脚本中另一个有用的函数是 formatIntersect 函数,它生成在运行 ArcToolbox Intersect 工具时使用的文件路径字符串。虽然自从该函数设计以来我们已经深入到 ArcPy 中,并且不再需要在我们的公交车站分析中调用 Intersect 工具,但这并不意味着我们将来永远不会需要调用它。它仍然很有用,应该添加到 useful.py 中。

我们可以访问的最后一个函数是 createCSV() 函数。从 Chapter4Modified.py 复制并粘贴到 useful.py 中。但是,为了避免需要单独导入 CSV 模块,我们需要稍微修改该函数。以下是它应该看起来像什么:

def createCSV(data, csvname, mode ='ab'):
 'creates a csv file'
 import csv
 with open(csvname, mode) as csvfile:
 csvwriter = csv.writer(csvfile, delimiter=',')
 csvwriter.writerow(data)
 del csv

通过导入然后删除 csv 模块,我们能够使用它来生成 csv 文件,并使用 del 关键字从内存中移除该模块。

现在我们已经将将要重用的函数保存在 useful.py 脚本中,在公共模块内部,让我们探索如何使用 Python 的 import 方法调用它们。打开 Python 可执行文件,使用 Python.exe 或 IDLE,或者在 Aptana 中的内置终端。在三个箭头提示符(>>>)中,输入以下行:

>>> from common.useful import createCSV>>>

如果出现第二个三重箭头形状的提示,则表示函数已正确从模块中导入。要在脚本中导入此模块中的函数,请使用相同的导入结构,并列出所需的函数,使用逗号分隔:

from common.useful import createCSV, formatSQLMultiple

脚本useful.py中的函数使用 Python 点符号调用。这是由于__init__.py文件告诉 Python,文件夹common现在是一个模块,并且它应该期望存在一个名为useful的方法,其中包含createCSVformatSQLMultiple函数。

添加高级分析组件

我们用来介绍 ArcPy 的公交车站分析可以进一步扩展以生成更精确的结果。为了更好地估计每个公交车站服务的真实人数,让我们添加一个函数,该函数将在考虑的区块内生成随机点,同时消除公园和其他不含住宅的区域。

要做到这一点,我们需要从San Francisco地理数据库中引入一个新的数据集,即RPD_Parks要素类。通过使用这个要素类来减少我们分析考虑的区域,我们可以为每个公交车站生成一个更现实的客运服务区人口评估。

当运行空间分析时,使用ArcToolbox 擦除工具擦除RPD_Parks多边形所表示的区域通常是一个常规步骤,但这个选项存在一些缺点。第一个缺点是擦除工具仅在 ArcGIS for Desktop 高级许可证级别下可用,这使得它仅对某些用户可用。第二个缺点是工具会产生一个中间数据集,在可能的情况下应尽量避免。

使用 ArcPy 将使我们能够避免这两个缺点。我们可以创建一个脚本,该脚本将仅在未与RPD_Parks要素类相交的普查区块多边形的一部分内生成随机点。为此,我们将深入到 ArcPy Polygon对象的方法中。

高级多边形对象方法

在第六章使用 ArcPy 几何对象中,我们开始探索 ArcPy 几何对象以及如何使用它们的方法来执行内存中的空间分析。介绍了缓冲区交集方法,并使用这些方法生成分析结果。接下来,我们将讨论更多这些方法,并展示它们如何有助于改进内存中的空间分析。

Polygon 对象有一个名为 Difference 的方法,允许我们找到两个多边形相交时的非相交区域。将 census block polygonpark polygon 作为参数传递将返回(作为多边形对象)第一个参数中不发生重叠的部分。另一个重要的方法是 Overlaps,它用于测试两个几何对象(点、线或多边形)是否相交。如果存在重叠,则 Overlaps 方法将返回 True,如果没有重叠,则返回 FalseUnion 也是一个重要的方法,将在本章中使用,它允许将两个几何对象 合并 成一个对象。

让我们探索这些重要的方法。为了找到两个多边形对象的非相交区域,以下函数结合了 OverlapsDifference 方法:

def nonIntersect(poly1,poly2):
 'returns area of non-intersect between two polygons'
 if poly1.overlaps(poly2) == True:
 return poly1.difference(poly2)

函数 nonIntersect 接受两个 Polygon 对象作为参数。第一个参数,poly1,是相交的多边形(人口普查区块多边形),第二个参数,poly2,是要检查重叠的多边形。if 条件语句使用 Overlaps 方法,如果两个参数之间存在重叠,则返回 True。如果存在任何重叠,则 difference() 方法将返回非相交区域作为多边形对象。然而,这个函数应该扩展以涵盖 Overlaps() 方法返回 False 的情况:

def nonIntersect(poly1,poly2):
 'returns area of non-intersect between two polygons'
 if poly1.overlaps(poly2) == True:
 return poly1.difference(poly2)
 else:
 return poly1

Overlaps 方法返回 False 时,函数现在将返回第一个参数,这表明两个多边形对象之间没有重叠。现在这个函数已经完成,可以在分析中使用。因为 nonIntersect() 是一个可以在其他空间分析中使用的函数,所以复制它并将其添加到 useful.py 中。

生成随机点以表示人口

改进公交车站分析的下一步是生成点以表示每个人口普查区块的人口。虽然随机点不能完美地表示人口,但它将作为良好的人口模型,并允许我们避免对每个由公交车站服务的区块进行面积平均以找到粗略的人口。ArcToolbox 数据管理工具集中的 CreateRandomPoints 工具使生成点变得简单。

工具 CreateRandomPoints 接受一系列必需和可选参数。由于该工具生成要素类,必需参数是要素类将被放置的工作空间和要素类的名称。感兴趣的可选参数是约束要素类和要生成的点数。由于我们希望在分析的中途步骤中避免创建新的要素类,我们可以利用 in_memory 工作空间,这意味着要素类将在内存中生成,这意味着它们不会被写入硬盘。

由于需要为每个普查区生成特定数量的随机点,我们应该创建一个函数,该函数将接受一个约束多边形和表示每个普查区的总人口数。然而,in_memory工作空间并不适用于所有情况,因此我们将提供工作空间参数的默认值:

def generatePoints(fc, pop,constrant, workspace='in_memory'):
 'generate random points'
 import os, arcpy
 arcpy.CreateRandomPoints_management(workspace, fc,                  constrant, "", pop, "")
 return os.path.join(workspace, fc)

该函数将在所需的工作空间中创建要素类,并将返回要素类的路径(使用os模块连接),以便在脚本的其他部分中使用。此函数也可重用,应将其复制到useful.py文件中。

在脚本中使用函数

现在我们已经创建了帮助我们在脚本中运行更高级空间分析的函数,让我们将它们添加到脚本中,并添加一些SearchCursors来遍历数据:

# Import the necessary modules
import arcpy, os
from common.useful import nonIntersect, generatePoints,createCSV

# Add an overwrite statement
arcpy.env.overwriteOutput = True

# Define the data inputs
busStops = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\Bus_Stops'
parks = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\RPD_Parks'
censusBlocks = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\CensusBlocks2010'
csvName = r'C:\Projects\Output\Chapter10Analysis.csv'

# Create the spreadsheet in memory and add field headers
headers = 'Line Name','Stop ID', 'Total Population Served'
createCSV(headers,csvName,mode='wb')

# Copy the census block data into a feature layer
arcpy.MakeFeatureLayer_management(censusBlocks,'census_lyr')

# Copy the park data geometries into a list and union them allparkGeoms = arcpy.CopyFeatures_management(parks,arcpy.Geometry())
parkUnion = parkGeoms[0]
for park in parkGeoms[1:]:
 parkUnion = parkUnion.union(park)

# Create a search cursor to iterate the bus stop data
sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops, ['NAME','STOPID','SHAPE@'],sql) as cursor:
 for row in cursor:

lineName = row[0]
 stopID = row[1]
 stop = row[2]
 busBuf = stop.buffer(400)
 # Select census blocks that intersect the bus buffer
 arcpy.SelectLayerByLocation_management("census_lyr","intersect", busBuf,'','NEW_SELECTION')
 # Use a second Cursor to find the selected population
 totalPopulation = 0
 with arcpy.da.SearchCursor("census_lyr",['SHAPE@','POP10',
 'BLOCKID10']) as ncursor:
 for nrow in ncursor:
 block = nrow[0]
 checkedBlock = nonIntersect(block, parkUnion)

blockName = nrow[2]
 population = nrow[1]
 if population != 0:
 points = generatePoints("PopPoints",
 population,checkedBlock)
 pointsGeoms = arcpy.CopyFeatures_management(points,arcpy.Geometry())
 pointsUnion = pointsGeoms[0]
 for point in pointsGeoms[1:]:
 pointsUnion = pointsUnion.union(point)
 pointsInBuffer=busBuf.intersect(pointsUnion, 1)
 intersectedPoints = pointsInBuffer.pointCount
 totalPopulation += intersectedPoints
 # Add the tallied data to the spreadsheet

data = lineName, stopID, totalPopulation
 print 'data written', data
 createCSV(data, csvName)

#Start the spreadsheet to see the results
os.startfile(csvName)

让我们逐节回顾代码,因为一开始要吸收的内容很多。

导入部分是我们调用常规模块,如 arcpy 和 os,以及common模块中的自定义函数:

import arcpy, os
from common.useful import nonIntersect
from common.useful import generatePoints
from common.useful import formatSQLMultiple
from common.useful import nonIntersectcreateCSV

如前所述,common模块的useful方法中的函数使用 Python 点符号和from … import ...导入风格调用,使它们直接可用。许多函数可以在一行中导入,用逗号分隔,或者像这里一样单独导入。

下一个设置 ArcPy 环境overwrite属性为True的行非常重要,因为它允许我们覆盖创建随机点操作的结果。如果不覆盖,函数结果将使用所有可用内存并导致脚本失败:

arcpy.env.overwriteOutput = True

注意

在此覆盖设置上要小心,因为它将允许覆盖任何要素类。我们所有的输出都在内存中,并且只为分析生成,所以这里几乎不需要担心,但运行脚本时要确保没有重要内容被覆盖。

下一个部分是将在脚本中使用的一组变量,并将初始化用于收集分析结果的电子表格:

busStops = r'C:\PacktDB.gdb\SanFrancisco\Bus_Stops'
parks = r'C:\PacktDB.gdb\SanFrancisco\RPD_Parks'
censusBlocks = r'C:\PacktDB.gdb\SanFrancisco\CensusBlocks2010'
csvName = r'C:\Projects\Output\Chapter10Analysis.csv'
headers = 'Line Name','Stop ID', 'Total Population Served'
createCSV(headers,csvName,mode='wb')

这里分配给变量的文件路径如果要将此转换为脚本工具,可以替换为 ArcPy 参数,但到目前为止,硬编码的路径是可行的。在变量下方,创建结果电子表格并添加列字段标题。

值得注意的是,电子表格是使用wb模式创建的。这种以二进制文件打开的模式称为wb(写二进制),用于创建新文件。必须显式传递给createCSV()函数作为默认模式参数是ab(追加二进制),如果不存在将创建新文件,或者添加到已存在的文件中(第三种二进制模式是rbread binary,用于打开现有文件)。

接下来的几行代码使特征类中的数据在内存中可用。人口普查区块数据被转换为Feature Layer,而RPD_Parks数据则被读取到内存中,作为一系列Polygon对象,然后这些对象被合并成一个单一的、统一的Polygon对象,称为parkUnion

arcpy.MakeFeatureLayer_management(censusBlocks,'census_lyr')parkGeoms = arcpy.CopyFeatures_management(parks,arcpy.Geometry())
parkUnion = parkGeoms[0]
for park in parkGeoms[1:]:
 parkUnion = parkUnion.union(park)

通过使用数据管理工具集中的CopyFeatures工具,parkGeoms变量接收RPD_Parks特征类中每行数据的几何形状列表。然而,我们不想必须遍历公园几何形状来与每个区块进行比较,因此调用Union方法从整个列表创建一个Polygon对象。通过将列表的第一个成员分配给parkUnion变量,然后遍历parkGeoms列表逐个合并其他几何形状,最终得到一个代表RPD_Parks数据集中所有公园的Polygon对象。

一旦所有模块都已导入并且变量已分配,我们就可以进入数据访问SearchCursorfor循环以开始分析。然而,我们不想对所有公交车站进行分析,因此我们将使用 SQL 语句的where子句来限制分析范围到单条公交线:

sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops, ['NAME','STOPID','SHAPE@'],sql) as cursor:
 for row in cursor:
 lineName = row[0]
 stopID = row[1]
 stop = row[2]
 busBuf = stop.buffer(400)
 arcpy.SelectLayerByLocation_management("census_lyr","intersect,busBuf,'','NEW_SELECTION')
 totalPopulation = 0

迭代的第一部分涉及进入for循环并将每行的值分配给一个变量。在SearchCursor返回的PointGeometry对象周围创建一个400英尺的Polygon对象缓冲区。然后,使用这个缓冲区与人口普查区块的Feature Layer相交,以找到并选择所有与缓冲区相交的区块。为了统计每个缓冲区服务的人口总数,创建了变量totalPopulation

一旦执行了选择操作,就可以使用第二个SearchCursor遍历所选区块,以检索它们的 人口值和用于随机点生成的Polygon对象:

with arcpy.da.SearchCursor("census_lyr",['SHAPE@','POP10',
 'BLOCKID10']) as ncursor:
 for nrow in ncursor:

block = nrow[0]
 checkedBlock = nonIntersect(block, parkUnion)
 blockName = nrow[2]
 population = nrow[1]

在这个迭代过程中,一旦检索到每个区块(以Polygon对象的形式),就使用之前创建的nonIntersect函数将该区块与合并的公园几何形状进行比对。这确保了点只会在非公园区域内创建,即更有可能代表人们居住的地方。同时也会检索人口值。

一旦评估了约束多边形(例如人口普查区块),移除了任何潜在的公园部分,并且获得了人口值,就可以使用generatePoints()函数生成随机点:

if population != 0:
 points = generatePoints("PopPoints",population,checkedBlock)
 pointsGeoms = arcpy.CopyFeatures_management(points,arcpy.Geometry())
 pointsUnion = pointsGeoms[0]
 for point in pointsGeoms[1:]:
 pointsUnion = pointsUnion.union(point)
 pointsInBuffer = busBuf.intersect(pointsUnion,1)
 intersectedPoints = pointsInBuffer.pointCount
 totalPopulation += intersectedPoints

generatePoints()函数需要三个参数。第一个是要生成的特征类的名称;每次生成时都会覆盖,从而避免为每个区块创建一个in_memory特征类,以避免过度使用内存。其他两个参数是人口值和约束多边形对象。

一旦这些参数传递给函数,它将返回新创建的特征类的文件路径,并将该文件路径分配给变量 points。然后使用CopyFeatures工具从 points 中提取几何形状,并将其分配给变量points。再次使用Union方法创建一个单一的、统一的PointGeometry对象,该对象将与公交车站缓冲区相交。一旦运行了此交集,结果几何形状将被分配给pointsInBuffer变量,并使用pointCount方法找到在缓冲区域内生成的点的数量。这是我们估计的人口普查区块内的人口,该值被添加到totalPopulation变量中,最终得出距离公交车站 400 英尺内的总估计人口。

脚本的最后几行演示了如何将数据收集到一个元组中,并将其传递给createCSV()模块以写入我们的最终电子表格:

 data = lineName, stopID,totalPopulation
 print 'data written', data
 createCSV(data, csvName)
os.startfile(csvName)

最后一行,os.startfile(csvName),使用os模块的startfile方法在分析完成后自动打开电子表格。在这种情况下,电子表格C:\Projects\Output\Chapter10Analysis.csv已填充了分析结果,并打开以显示这些结果。然而,用户可能需要指示这些行是逗号分隔值,才能打开脚本。

在脚本中使用函数

而不是创建逗号分隔值,我们可以利用 ArcGIS 10.2 和 ArcPy 安装时安装的另一个 Python 模块。这个模块被称为XLWT,用于生成 Excel 电子表格,并且与 Excel 电子表格读取模块XLRD一起,是 Python 用户可用最有用的模块之一。

使用 XLWT 创建 XLS

XLWT 是一个功能强大的模块,允许多种样式选项。然而,对于我们的目的,我们可以忽略这些选项,创建一个将生成我们的空间分析结果的电子表格的函数。当然,此函数可以添加到common.useful中:

def generateXLS(indatas, sheetName, fileName):
 import xlwt
 workbook = xlwt.Workbook()
 sheet = workbook.add_sheet(sheetName)
 for YCOUNTER, data in enumerate(indatas):
 for XCOUNTER, value in enumerate(data):
 sheet.write(YCOUNTER, XCOUNTER, value)
 workbook.save(fileName)

此函数需要三个参数,indatas-一个包含可迭代数据行的列表,一个字符串工作表名称,以及一个以.xls扩展名结尾的字符串文件名。

要使用此函数,将其添加到common.useful中。一旦添加,复制并重命名较旧的脚本分析,以便进行调整:

import arcpy, os
from common.useful import nonIntersect, generatePoints, generateXLS

arcpy.env.overwriteOutput = True

busStops = r'C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops'
parks = r'C:\Projects\PacktDB.gdb\SanFrancisco\RPD_Parks'
censusBlocks = r'C:\Projects\PacktDB.gdb\SanFrancisco\CensusBlocks2010'
xlsName = r'C:\Projects\Output\Chapter10Analysis.xls'

headers = 'Line Name','Stop ID', 'Total Population Served'
indatas = [headers]

arcpy.MakeFeatureLayer_management(censusBlocks,'census_lyr')parkGeoms = arcpy.CopyFeatures_management(parks,arcpy.Geometry())
parkUnion = parkGeoms[0]
for park in parkGeoms[1:]:

parkUnion = parkUnion.union(park)

sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops, ['NAME','STOPID',
 'SHAPE@'],sql) as cursor:
 for row in cursor:
 lineName = row[0]
 stopID = row[1]
 stop = row[2]
 busBuf = stop.buffer(400)
 arcpy.SelectLayerByLocation_management("census_lyr","intersect",busBuf,'','NEW_SELECTION')
 totalPopulation = 0
 with arcpy.da.SearchCursor("census_lyr", ['SHAPE@','POP10',
 'BLOCKID10']) as ncursor:
 for nrow in ncursor:

block = nrow[0]
 checkedBlock = nonIntersect(block, parkUnion)
 blockName = nrow[2]
 population = nrow[1]
 if population != 0:
 points = generatePoints("PopPoints",population,checkedBlock)

 pointsGeoms = arcpy.CopyFeatures_management(points,arcpy.Geometry())
 pointsUnion = pointsGeoms[0]
 for point in pointsGeoms[1:]:
 pointsUnion = pointsUnion.union(point)
 pointsInBuffer = busBuf.intersect(pointsUnion,1)

intersectedPoints = pointsInBuffer.pointCount
 totalPopulation += intersectedPoints
 data = lineName, stopID, totalPopulation
 indatas.append(data)
generateXLS(indatas, "Results", xlsName)
os.startfile(xlsName)

我们现在可以像生成CSV文件一样轻松地生成 Excel 电子表格,同时使用可重用的函数。我们现在能够快速执行可重复的空间分析,并以行业标准格式生成结果。

摘要

在本章中,我们探讨了如何创建模块和可重用的函数,这些函数将节省未来的脚本编写时间,使我们能够避免重写这些有用的函数。我们还进一步探讨了通过 ArcPy 几何对象可用的方法,包括 Intersect(交集)、Overlaps(重叠)和 Union(并集)方法。我们创建了一个不将要素类写入磁盘的空间分析,这样分析时间就减少了,并且避免了不必要的文件。最后,我们探讨了如何使用 XLWT 模块生成 Excel 电子表格,以便分析结果可以以行业标准格式共享。

在下一章中,我们将探讨如何使用 ArcPy 与 ArcGIS for Desktop 的扩展程序如网络分析师和空间分析师进行交互。通过在脚本中整合它们的功能,我们进一步增强了创建快速且可重复的空间分析工作流程的能力。

第十一章:使用 ArcPy 的网络分析师和空间分析师

使用 ArcGIS for Desktop 扩展也受益于 Python 和 ArcPy 的强大功能。使用街道数据集或公交车路线数据集通过 ArcPy 建模路线的能力将帮助我们把整个工作流程转换为脚本工具。网络分析师和空间分析师都内置了 ArcPy 的访问模块,以增强对可用工具、方法和属性的掌控。

本章将涵盖以下主题:

  • 创建简单的网络数据集

  • 检查扩展

  • ArcPy 网络分析师模块

  • ArcPy 空间分析师模块

网络分析师扩展

ESRI 的网络分析师扩展是一个强大的工具,可以在 ArcGIS 中启用路由和网络连接功能。当用于街道路由时,该扩展允许用户在道路网络中找到两点之间的最短路径。路线可以由多个因素限制,例如交通或左转,以更好地模拟道路旅行。使用其他类型的网络,如水管网络或电网,也可以运行类似的分析。

使用网络分析师

要使用网络分析师扩展,需要 ArcGIS for Desktop 高级许可证。在 ArcCatalog 或 ArcMap 中,点击自定义菜单并选择扩展。一旦扩展菜单打开,点击旁边的复选框以启用网络分析师扩展

创建要素数据集

使用网络数据集的第一步是在要素数据集中创建一个。为此,我们将生成一个要素数据集来保存感兴趣的数据。在包含公交站数据的文件地理数据库上右键单击,然后从新建菜单中选择要素数据集。将其命名为Chapter11Results并点击下一步

创建要素数据集

接下来,选择空间参考系统SRS)。在这种情况下,我们将使用旧金山的本地州平面区域的空间参考系统。它是一个投影坐标系统,因此选择该文件夹,然后点击State Plane文件夹。一旦打开,选择名为NAD 1983(US Feet)的文件夹。从可用的参考系统中,选择名为NAD 1983 StatePlane California III FIPS 0403 (US Feet)的选项。点击下一步进入下一个菜单。

注意

该系统也被称为 2227,在已知 IDWKID)或欧洲石油调查组EPSG)系统中。有关这些代码的更多信息可在spatialreference.org网站找到,该网站用于查找全球使用的数千个空间参考系统。

点击垂直坐标系统文件夹,然后选择北美洲文件夹。选择北美洲 1988 年垂直基准(英尺)NAVD 1988 US survey feet)。这将使得垂直和水平线性单位处于相同的测量系统。点击下一步进入下一个菜单。

下一页上的容限也非常重要,但在此处我们不会详细说明。接受默认设置并点击完成以最终确定要素数据集

导入数据集

将公交车站、街道和公交路线要素类导入到Chapter 11 Results Feature Dataset中。右键单击数据集并选择导入,然后要素类(单个)。逐个添加要素类,并给它们一个新名字,以保持它们与SanFrancisco Feature Dataset中包含的版本分开。导入它们将确保它们位于正确的 SRS 中,并且可以创建网络数据集。

创建网络数据集

现在我们有了数据容器,我们可以从街道要素类创建一个网络数据集。右键单击Chapter11Results要素数据集并选择新建,然后选择网络数据集

创建网络数据集

网络数据集命名为Street_Network并点击下一步。选择街道要素类作为将参与网络数据集的类,并点击下一步进入下一个菜单。选择全局转弯来模拟网络内的转弯。在下一个菜单中,使用默认的连通性设置。然后,接受使用几何中的 Z 坐标值设置。接受默认的成本限制和驾驶方向设置,最后点击完成以生成网络数据集。然后,使用最终菜单构建网络数据集。网络数据集已准备好使用。

使用 ArcPy 访问网络数据集

现在必要的设置已经完成,可以将 street_network 网络数据集添加到脚本中以生成路线。因为这是一个简单的分析,将使用的唯一阻抗值将是街道段落的长度。通过使用SearchCursor,可以访问并添加来自公交车站的PointGeometry对象作为要搜索的位置:

import arcpy
arcpy.CheckOutExtension("Network")
busStops = r'C:\Projects\PacktDB.gdb\Chapter11Results\BusStops'
networkDataset = r'C:\Projects\PacktDB.gdb\Chapter11Results\street_network'
networkLayer = "streetRoute"
impedance = "Length"
routeFile = "C:\Projects\Layer\{0}.lyr".format(networkLayer)
arcpy.MakeRouteLayer_na(networkDataset,
 networkLayer, impedance)
print 'layer created'
sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops,['SHAPE@', 'STOPID'],sql) as cursor:
 for row in cursor:
 stopShape = row[0]
 print row[1]
 arcpy.AddLocations_na(networkLayer,'Stops',stopShape, "", "") 
 arcpy.Solve_na(networkLayer,"SKIP")
arcpy.SaveToLayerFile_management(networkLayer,routeLayerFile,"RELATIVE")
print 'finished'

脚本分解

让我们剖析这个脚本,一旦完成,它将生成一个包含添加的停靠点和街道上的路线的图层文件,以最佳方式从起点停靠站到达目的地停靠站。

脚本首先导入 arcPy 模块。下一行允许我们使用网络分析扩展:

arcpy.CheckOutExtension("Network")

使用arcpy.CheckOutExtension()方法调用网络分析扩展时,需要将正确的关键字作为参数传递给方法。一旦调用,扩展的工具就可以在脚本中调用和执行。

将公交车站要素类和街道网络网络数据集分配给变量后,它们可以传递给 ArcPy 的 MakeRouteLayer_na() 方法,以及一个表示阻抗值的变量:

arcpy.MakeRouteLayer_na(networkDataset,
 networkLayer, impedance)

MakeRouteLayer_na 工具在内存中生成一个 RouteLayer。这个空白层需要填充停靠点以生成它们之间的路线(s)。为此,我们需要一个 SearchCursor 来访问 PointGeometry 对象,以及一个 SQL 语句,该语句将返回的结果限制在感兴趣的线上:

sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops,['SHAPE@', 'STOPID'],sql) as cursor:
 for row in cursor:
 stopShape = row[0]
 print row[1]
 arcpy.AddLocations_na(networkLayer,'Stops',stopShape,"", "")

搜索游标将允许当与 AddLocations 工具结合使用时,填充由 MakeRouteLayer 工具生成的层的 Stops 子层。一旦填充,RouteLayer 可以传递给 Solve 工具以找到感兴趣点之间的路线。再次强调,路线是通过在两点之间找到最低的 阻抗 来解决的。在这个例子中,唯一的阻抗是段长度,但如果有数据可用,它可能是交通、海拔或其他限制类型:

 arcpy.Solve_na(networkLayer,"SKIP")
arcpy.SaveToLayerFile_management(networkLayer,routeLayerFile,"RELATIVE")

最终结果是使用 SaveToLayerFile 工具写入磁盘的图层文件。

分解脚本

网络分析模块

为了使网络分析扩展的使用更加 Pythonic,较新的网络分析器 (na) 模块调整了访问与 ArcToolbox 网络分析工具相对应的方法的方式。现在,不是直接从 ArcPy 调用工具,而是将工具作为 na 模块的方法。移除网络分析工具集的首字母缩写词也减少了混淆,并使得记住方法名称更容易。如下所示,可以看到差异:

import arcpy
arcpy.CheckOutExtension("Network")
busStops = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\Bus_Stops
networkDataset = r'C:\Projects\SanFrancisco.gdb\Chapter11Results\street_network'
networkLayer = "streetRoute"
impedance = "Length"
routeLayerFile = "C:\Projects\Layer\{0}_2.lyr".format(networkLayer)arcpy.na.MakeRouteLayer(networkDataset, networkLayer,impedance)
print 'layer created'
sql = "NAME = '71 IB' AND BUS_SIGNAG = 'Ferry Plaza'"
with arcpy.da.SearchCursor(busStops,['SHAPE@','STOPID'],sql) as cursor:
 for row in cursor:
 stopShape = row[0]
 print row[1]
 arcpy.na.AddLocations(networkLayer,'Stops', stopShape, "", "")
arcpy.na.Solve(networkLayer,"SKIP")
arcpy.management.SaveToLayerFile(networkLayer,routeLayerFile,"RELATIVE")
print 'finished'

工具将产生与原始脚本相同的层输出,但将网络分析工具重新组织到 na 模块中使得代码更加逻辑。例如,使用 arcpy.na.Solve() 调用 Solve 比使用 arcpy.Solve_na() 更有意义,因为这强调了 Solve 是网络分析器 (na) 模块的方法。随着 ArcPy 的不断发展,我预计将发生更多的 Pythonic 代码重组。

访问空间分析扩展

空间分析扩展对于对栅格和矢量数据集进行分析非常重要,但它通常用于执行地表分析和栅格数学。通过使用 ArcPy,这些操作变得更加容易,因为空间分析工具箱中所有可用的工具都通过空间分析访问模块公开。这包括栅格计算器工具,通过使用简单表达式中的工具和运算符,使得地图代数变得简单。

将海拔添加到公交车站

高程栅格 "sf_elevation" 已从 NOAA 下载并添加到文件地理数据库中。然而,它覆盖了整个湾区,我们应该编写一个脚本来仅提取旧金山城市的某个区域,这将减少我们脚本运行所需的时间。我们将使用 SQL 语句作为 where 子句来限制结果仅限于南市场(SoMa)街区。为此,让我们利用搜索游标和空间分析访问模块的 Extract by Polygon 属性:

import arcpy
arcpy.CheckOutExtension("Spatial")
busStops = r'C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops'
sanFranciscoHoods = r'C:\Projects\PacktDB.gdb\SanFrancisco\SFFind_Neighborhoods'
sfElevation = r'C:\Projects\PacktDB.gdb\sf_elevation'
somaGeometry = []
sql = "name = 'South of Market'"
with arcpy.da.SearchCursor(sanFranciscoHoods,['SHAPE@XY'],sql,None, True) as cursor:
 for row in cursor:
 X = row[0][0]
 Y = row[0][1]
 somaGeometry.append(arcpy.Point(X,Y))
somaElev = arcpy.sa.ExtractByPolygon(sfElevation,somaGeometry,"INSIDE")
somaOutPath = sfElevation.replace('sf_elevation','SOMA_elev')
somaElev.save(somaOutPath)
print 'extraction finished'

ExtractByPolygon() 方法有些误导,因为它不接受 Polygon 对象作为参数。相反,它需要一个表示我们想要提取的区域顶点的 Point 对象列表。当 SearchCursor 遍历街区数据集时,游标会返回一个 Polygon 对象。幸运的是,SearchCursor 有一个我们尚未探索的最终参数,允许我们提取构成 Soma 街区多边形的单个点或顶点。通过将搜索游标的可选 Explode to Points 参数(它将多边形对象转换为每个顶点的坐标对)设置为 True,可以通过将每个返回顶点的 XY 值传递给 arcpy.Point 方法来生成点对象。这些 Point() 对象被追加到 somaGeometry 列表中,然后传递给空间分析访问模块的 ExtractByPolygon 方法。

注意

使用多边形对象而不是点对象将返回一个错误。

使用地图代数生成英尺高程

我们现在有一个栅格可以用来提取高程值。然而,原始栅格和生成的 SoMa 街区栅格都包含以米为单位的高程值,最好将它们转换为英尺,以保持与公交车站投影的一致性。让我们使用栅格数学和 Times() 方法将值从米转换为英尺:

somaOutPath = sfElevation.replace('sf_elevation','SOMA_elev')
outTimes = arcpy.sa.Times(somaOutPath, 3.28084)
somaFeetOutPath = sfElevation.replace('sf_elevation','SOMA_feet')
outTimes.save(somaFeetOutPath)

Times() 方法生成一个新的栅格,以获取我们需要的感兴趣公交车站的地面高程值。

添加公交车站和获取高程值

现在我们已经生成了一个可以用来查找英尺高程值的栅格,我们需要添加一个新的 arcpy.sa() 方法来生成点。ExtractValuesToPoints() 方法将生成一个新的公交车站要素类,其中包含一个新的字段,用于存储高程值:

with arcpy.da.SearchCursor(sanFranciscoHoods,['SHAPE@'],sql) as cursor:
 for row in cursor:
 somaPoly = row[0]
arcpy.MakeFeatureLayer_management(busStops, 'soma_stops')
arcpy.SelectLayerByLocation_management("soma_stops", "INTERSECT",somaPoly)
outStops = r'C:\Projects\PacktDB.gdb\Chapter11Results\SoMaStops'
arcpy.sa.ExtractValuesToPoints("soma_stops", somaOutFeet,outStops,"INTERPOLATE","VALUE_ONLY")
print 'points generated'

最终结果

我们生成了一个包含高程值字段的公交车站子要素类。这个过程可以重复应用于整个城市,一次一个街区,或者可以在原始高程栅格的整个公交车站要素类上执行,为每个车站生成一个值:

import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
busStops = r'C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops'
sanFranciscoHoods = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\SFFind_Neighborhoods'
sfElevation = r'C:\Projects\SanFrancisco.gdb\sf_elevation'

somaGeometry = []
sql = "name = 'South of Market'"
with arcpy.da.SearchCursor(sanFranciscoHoods,['SHAPE@XY'],sql,None, True) as cursor:
 for row in cursor:
 somaGeometry.append(arcpy.Point(row[0][0],row[0][1]))
somaElev = arcpy.sa.ExtractByPolygon(sfElevation, somaGeometry,"INSIDE")
somaOutput = sfElevation.replace('sf_elevation','SOMA_elev')
somaElev.save(somaOutput)
print 'extraction finished'

somaOutput = sfElevation.replace('sf_elevation','SOMA_elev')
outTimes = arcpy.sa.Times(somaOutput, 3.28084)
somaOutFeet = sfElevation.replace('sf_elevation','SOMA_feet')
outTimes.save(somaOutFeet)
print 'conversion complete'

with arcpy.da.SearchCursor(sanFranciscoHoods,['SHAPE@'],sql) as cursor:
 for row in cursor:
 somaPoly = row[0]

arcpy.MakeFeatureLayer_management(busStops, 'soma_stops')
arcpy.SelectLayerByLocation_management("soma_stops", "INTERSECT",somaPoly)

outStops = r'C:\Projects\SanFrancisco.gdb\Chapter11Results\SoMaStops'
arcpy.sa.ExtractValuesToPoints("soma_stops", somaOutFeet,outStops,"INTERPOLATE","VALUE_ONLY")
print 'points generated'

此脚本很好地展示了访问 ArcPy 中高级扩展的价值,并将它们与 SearchCursorsGeometry 对象结合使用。该脚本可以通过添加一个 SearchCursor 来进一步扩展,以遍历 outstops 数据集并导出结果到电子表格,或者甚至向原始公交车站数据集添加一个新字段以填充海拔值。它甚至可以用作阻抗值输入到网络分析师扩展分析中——这是一个有趣的编码任务,我希望你们会尝试。

摘要

在本章中,我们介绍了在 ArcPy 中使用常见的 ArcGIS for Desktop Advanced 扩展的基础知识,重点关注网络分析师访问模块和空间分析师访问模块。我们探讨了如何生成网络以及如何使用 ArcPy 创建网络路径。我们还探讨了如何访问空间分析师工具,并使用它们与 SearchCursors 结合来处理栅格和矢量数据以进行空间分析。

在下一章中,我们将探讨 ArcPy 的一些最终拼图,这将允许创建高级脚本和脚本工具。

第十二章。开端的结束

本书几乎完成,但关于在 Python 和 ArcPy 中编写代码还有很多东西要了解。不幸的是,我无法将其全部放入一本书中,但这也意味着你可以享受探索 ArcPy 的所有方法和属性。作为本书的结论,我们将涵盖一些在编写 ArcPy 脚本时可能出现的其他重要主题。结合前面章节的教训,我希望你很快就能在工作中、在学校里或只是为了乐趣(为什么不呢?)使用 ArcPy。

本章将涵盖以下主题:

  • 处理字段信息 - 类型、别名、域、空间类型等

  • 访问描述要素类的信息

  • 自动生成要素类并填充字段

  • 自动创建文件地理数据库和要素数据集

  • 创建一个脚本工具,该工具将运行公交车站分析,并在自动生成的文件地理数据库、要素数据集和要素类中生成结果

从要素类获取字段信息

在创建脚本工具或仅运行脚本时,有时需要从要素类(或形状文件)中提取字段信息。这些信息可以包括字段名称和别名、字段类型和长度、比例、域或子类型。这些都是通过 arcpy.ListFields 方法可用的属性。我们将探讨许多属性,如何提取它们,以及如何在脚本中使用它们。

通过将 ArcPy 方法组织成函数,数据将以我们偏好的形式组织,而不是依赖于 ArcPy 设计者使用的默认组织方式。重要的是要记住,你创建的脚本应该反映你的需求,而创建这些函数包装器是向将原始 ArcPy 工具精炼以适应你的工作流程迈出的第一步。

访问 ListFields 的属性

列表字段工具作为 ArcPy 方法提供。Arcpy.ListFields 只接受一个参数,即要素类或形状文件。一旦参数传递成功,就可以使用点符号访问一系列重要属性。为了进一步利用这些属性,我们将创建函数,使其能够轻松获取所需的信息,并按照所需的格式进行。

列表推导式

在这些字段信息函数中,我们将利用一种名为列表推导式的 Python 数据结构。它们简化了 for 循环结构,使其更容易填充所需值(在这种情况下为字段信息)的列表。

要创建列表推导式,需要在括号内生成一个 for 循环,并将列表填充生成的值。以下是一个示例,它创建了一个包含从 1 到 10 的数字的平方值的列表,在 Python 解释器中运行:

>>>originalList = range(1,11)
>>>print originalList
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
>>>newList =  [x**2 for x in originalList]
>>>print newList
[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

使用列表推导式是因为它们更快、更容易编写,尽管可能需要一些时间来习惯语法。通过实验来更好地理解它们的用途和限制,并查阅网上可用的许多资源。

创建字段信息函数

每个函数都将是一个独立的实体,但它们都将具有相似的结构。每个函数都将接受一个参数,即感兴趣的要素类。ArcPy 将被导入,并在稍后从内存中删除,以确保ListFields()方法可以无错误地调用。一旦要素类被传递给ListFields()方法,所需值将填充在列表推导式内的一个列表中。一旦填充完成,它将通过使用return关键字从函数中返回。

这里是字段名称的函数集:

def returnfieldnames(fc):
 import arcpy
 fieldnames = [f.name for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldnames

def returnfieldalias(fc):
 import arcpy
 fieldalias = [f.aliasName for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldalias

def returnfieldbasename(fc):
 import arcpy
 fieldtypes = [f.baseName for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldtypes

这些名称函数在基于另一个要素类创建新要素类时非常有用。有时需要保留原始要素类的确切名称和别名,使用这些函数可以使这一点成为可能。在这种情况下,还需要提供其他字段信息。以下是与字段类型、长度、精度和比例相关的函数:

def returnfieldtypes(fc):
 import arcpy
 fieldtypes = [f.type for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldtypes

def returnfieldlength(fc):
 import arcpy
 fieldlengths = [f.length for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldlengths

def returnfieldprecision(fc):
 import arcpy
 fieldprecise = [f.precision for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldprecise

def returnfieldscale(fc):
 import arcpy
 fieldscales = [f.scale for f in arcpy.ListFields(fc)]
 del arcpy
 return fieldscales

甚至有一个属性用于请求域信息:

def returnfielddomain(fc):
 import arcpy
 fielddomains = [f.domain for f in arcpy.ListFields(fc)]
 del arcpy
 return fielddomains

这些函数都具有前面讨论的结构,并且具有使用简单、易于搜索的优点。因为要素类中的字段有特定的顺序,所以每个函数返回的列表都将有一个信息返回的顺序,可以通过特定的索引号访问。

fieldsubtypes也通过数据访问模块提供。因为它们与字段相关,所以作为字典返回:

def returnfieldsubtypes(fc):
 import arcpy
 fieldsubdic = {}
 subtypes = arcpy.da.ListSubtypes(fc)
 for stcode, stdict in subtypes.iteritems():
 for stkey in stdict.iterkeys():
 if stkey == 'FieldValues':
 fields = stdict[stkey]
 for field, fieldvals in fields.iteritems():
 sub = fieldvals[0]
 desc = fieldvals[1]
 fieldsubdic[field] = sub, desc

 del arcpy
 return fieldsubdic

注意

将这些函数添加到公共模块中的useful.py脚本中,将使它们对任何脚本或脚本工具可用。使用import关键字将它们添加到任何新脚本中。它们是自包含的函数,只需要提供感兴趣要素类的文件路径。

查询要素类信息

使用ListFields()方法无法访问有关传入要素类的一些重要信息。相反,将使用许多不同的方法来查找每个要素类的几何类型、空间参考或字段子类型。其中一些是通过 ArcPy 的Describe方法发现的,该方法旨在提供

对于几何类型,我们将使用Describe()方法的shapeType属性:

def returngeometrytype(fc):
 import arcpy
 arcInfo = arcpy.Describe(fc)
 geomtype = arcInfo.shapeType
 del arcpy
 return str(geomtype)

Shape 字段的名称(通常默认为 Shape)也可以通过Describe方法请求,并返回一个字符串数据类型:

def returngeometryname(fc):
 import arcpy
 arcInfo = arcpy.Describe(fc)
 geomname = arcInfo.shapeFieldName
 del arcpy
 return str(geomname)

通过Describe方法,要素类的spatial_reference也是可用的。数据以spatial_reference对象的形式返回:

def returnspatialreference(fc):
 import arcpy
 spatial_reference = arcpy.Describe(fc).spatialReference
 del arcpy
 return spatial_reference

spatial_reference对象具有许多重要的属性。projectionnameprojectioncode是其中重要的属性

def returnprojectioncode(fc):
 import arcpy
 spatial_reference = arcpy.Describe(fc).spatialReference
 proj_code = spatial_reference.projectionCode
 del arcpy
 return  proj_code

def returnprojectionname(fc):
 import arcpy
 spatial_reference = arcpy.Describe(fc).spatialReference
 proj_name = spatial_reference.name
 del arcpy
 return  proj_name

许多其他属性和方法可以类似地利用,使它们在脚本或脚本工具中可用。探索 ArcGIS 帮助文档,以获取通过Describe方法可用的属性的更多信息。

生成文件地理数据库和特征类

文件地理数据库不需要在运行脚本之前存在;相反,它们可以在使用CreateFileGDB工具执行脚本时生成,该工具也是 ArcPy 的一个方法。一旦创建了File地理数据库,就可以添加Feature Datasets

生成文件地理数据库非常简单。唯一的参数是将它放置在内的文件夹和地理数据库的名称:

import arcpy
folderPath = r"C:\Projects"
gdbName = "ArcPy.gdb"
arcpy.CreateFileGDB_management(folderPath, gdbName)

特征数据集的创建比较困难,因为有一个可选的空间参考参数,需要生成空间参考对象。虽然空间参考对象是可选的,但强烈建议使用。

有几种选项可以生成SpatialReference对象。其中之一是使用之前定义的return specialReference()函数;通过向函数传递一个特征类,创建一个空间参考对象。另一种方法是将投影文件的文件路径.prj作为可选的第三个参数传递。第三种方法是通过使用arcpy.SpatialReference方法并传递投影代码或投影字符串来生成空间参考对象:

spatialReference = arcpy.SpatialReference(2227)

无论如何生成,它随后都会传递给arcpy.CreateFeatureDataset方法,包括文件地理数据库的文件路径和特征数据集的名称:

spatialReference = arcpy.SpatialReference(2227)
fileGDB = r"{0}\{1}".format(folderPath,gdbName)
featureDataset = "Chapter12Results"
arcpy.CreateFeatureDataset_management(fileGDB, featureDataset,   spatialReference) 

生成特征类

现在已经创建了一个文件地理数据库和特征数据集,让我们在特征数据集中生成一个特征类。这是通过使用arcpy.CreateFeatureClass方法完成的。此方法有许多可选参数,包括用作模板的特征类和一个空间参考。对于这个例子,不需要使用空间参考参数,因为它将被写入特征数据集,这决定了使用哪个空间参考。模板参数将复制模板特征类的字段,但到目前为止,我们只创建形状字段:

featureClass = "BufferArea"
geometryType = "POLYGON"
featurePath = r"{0}\{1}".format(fileGDB,featureDataset)
arcpy.CreateFeatureclass_management(featurePath, featureClass,  geometryType) 

创建的特征类需要一些字段,这些字段将包含稍后填充的属性信息。字段有许多参数,取决于字段类型,包括长度、精度和别名等:

fieldName = "STOPID"
fieldAlias = "Bus Stop Identifier"
fieldType = "LONG"
fieldPrecision = 9
featureClassPath = r"{0}\{1}".format(featurePath,featureClass)
arcpy.AddField_management(featureClassPath, fieldName,fieldType, fieldPrecision,
"", "", fieldAlias)

让我们添加一个字段来存储公交车站分析产生的平均人口值:

fieldName2 = "AVEPOP"
fieldAlias2 = "Average Census Population"
fieldType2 = "FLOAT"
featureClassPath = r"{0}\{1}".format(featurePath,featureClass)
arcpy.AddField_management(featureClassPath, fieldName2, fieldType2, "", "", "", fieldAlias2)

文件地理数据库、特征数据集和特征类字段现在已经生成。让我们通过添加公交车站分析函数并将结果写入生成的特征类来将脚本扩展为脚本工具。创建一个填充特征类的脚本工具。

此脚本工具将借鉴第十章中概述的思想,高级几何对象方法,并将创建与缓冲公交车站相交的 Polygon 几何对象集合,以填充Shape 字段,以及与每个缓冲区相交的区块的公交车站 ID和平均人口。

打开脚本Chapter12_3.py并探索其内容。结合前面提到的代码片段和arcpy.GetParameterAsText的使用来从脚本工具获取数据,生成的数据将通过以下代码写入要素类:

arcpy.AddMessage("Beginning Analysis")
insertCursor = arcpy.da.InsertCursor(featureClassPath,['SHAPE@',fieldName, fieldName2])
arcpy.MakeFeatureLayer_management(censusBlocks2010,"census_lyr")

with arcpy.da.SearchCursor(busStops, ['SHAPE@', busStopField],sql) as cursor:
 for row in cursor:
 stop = row[0]
 stopID = row[1]
 busBuffer = stop.buffer(400)
 arcpy.SelectLayerByLocation_management("census_lyr","intersect",busBuffer,'','NEW_SELECTION')
 censusShapes = []
 censusPopList = []
 with arcpy.da.SearchCursor("census_lyr", ['SHAPE@',censusBlockPopField]) as ncursor:
 for nrow in ncursor:
 censusShapes.append(nrow[0])
 censusPopList.append(nrow[1])

 censusUnion = censusShapes[0]
 for block in censusShapes[1:]:
 censusUnion = censusUnion.union(block)

 censusPop = sum(censusPopList)/len(censusPopList)
 finalData = (censusUnion,stopID, censusPopulation)
 insertCursor.insertRow(finalData)
arcpy.AddMessage("Analysis Complete")

该脚本结合了本书中介绍的一些想法,使用户能够运行一个完整的流程,生成包含分析结果的要素类。通过仅添加感兴趣的字段,并用合并的Polygon对象填充它们,脚本消除了在运行空间分析时通常创建的大部分冗余,并生成可以在 ArcMap 中查看的结果数据集。

设置脚本工具参数

这是设置脚本工具参数时的样子:

设置脚本工具参数

参数列表很长,所以我使用了两张图片来展示它们。为每个参数选择正确的数据类型很重要,因为它将控制用于检索数据的对话框。

公交车站 ID 字段人口字段都是从它们各自的要素类中获取的。文件地理数据库名称是一个字符串,如果最初没有输入,代码将向输入字符串的末尾追加.gdb,以确保它可以正确生成。它不应已经存在;如果它存在,则不会生成(如果需要,可以在import语句之后将arcpy.env.overwriteOutput属性设置为True来更改此行为)。

一旦设置了参数,并且工具有了名称和描述,保存它然后打开工具。一旦填写完毕,它应该看起来像这样:

设置脚本工具参数

点击确定以运行工具。打开ArcMap并添加结果,以及来自第四章的旧金山多边形和Inbound71要素类,高级 ArcPy 脚本和泛化函数。经过一些制图符号化后,结果将类似于以下内容:

设置脚本工具参数

最终结果将按每个选定的公交车站有一行,包括平均人口和公交车站 ID 值。而不是使用电子表格作为输出,要素类将允许制作地图或进行进一步的空间分析。使用自定义脚本工具生成自定义数据使你在进行地理空间分析时处于驾驶员的位置,并使你的工具和 yourself 成为任何团队的有价值资产。

环境设置

ArcPy 模块允许通过 ArcPy 的env类控制全局设置,这些设置控制输入和输出过程。这些设置将对使用地理空间分析工具产生数据的精度产生影响。可以控制XYZM坐标的分辨率和公差,以及输出范围、栅格单元格大小、分析工作空间和许多其他设置。

要使用 ArcPy 访问环境设置,需要从arcpy导入env类:

>>> from arcpy import env

它也可以使用以下点表示法调用。设置工作空间消除了在后续调用workspace上的任何方法时传递文件路径的需要。以下是一个设置工作空间并调用ListDatasets()方法的示例,而不将文件路径作为参数传递:

>>> import arcpy
>>> arcpy.env.workspace = r"C:\Projects\SanFrancisco.gdb"
>>> arcpy.ListDatasets()
[u'SanFrancisco', u'Chapter3Results', u'Chapter4Results', u'Chapter5Results', u'Chapter7Results', u'Chapter11Results']

分辨率和公差设置

分辨率和公差设置控制由 ArcToolbox 中的工具或使用 ArcPy 运行脚本产生的任何数据的输出精度。这些设置(并且应该)为文件地理数据库或企业地理数据库中的要素数据集设置,但在内存中运行分析或使用形状文件时,或者如果地理空间分析需要比那些地理数据库使用的精度更高时,设置它们非常重要。

设置分辨率和公差需要了解您项目所需的精度。这些设置将限制对齐到线或找到与线相交的点的能力。线性单位需要反映所选的坐标系:

import arcpy
arcpy.env.MResolution = 0.0005
arcpy.env.MTolerance = 0.005
arcpy.env.ZResolution = "0.0025 Feet"
arcpy.env.ZTolerance = "0.001 Feet"
arcpy.env.XYResolution = "0.00025 Feet"
arcpy.env.XYTolerance = "0.0005 Feet"

其他重要的环境设置包括:

  • 范围设置,通过使用范围对象设置感兴趣区域的矩形,或使用当前坐标系中用空格分隔的坐标(XminYminXmaxYmax)的字符串来限制由分析产生的任何数据的范围。

  • 掩膜设置,将栅格分析限制为与要素类或作为字符串文件路径参数传递给设置的栅格相交的区域。

  • 单元格大小设置,控制使用栅格分析产生的数据的单元格大小。

抽时间探索强大的 ArcPy 环境设置,以减少编写代码所需的时间并确保高质量的数据生产。

摘要

本章和本书展示了 ArcPy 可以用于自动化地理空间分析的许多方法之一。通过应用这些课程,并通过在 ArcPy 的许多方法和属性上发挥创意,可以编写脚本并将重复且缓慢的地理空间过程制作成自定义工具,这将节省大量时间。

我希望您喜欢学习使用 ArcPy 和 Python 的基本脚本知识。我真心希望您甚至开始喜欢编程的想法,因为它强大且赋权。还有更多需要掌握的内容,但我认为您会发现,您编写的脚本越多,理解起来就越容易。

深入理解 ArcPy 的最佳资源是 ArcGIS 帮助文档,您可以通过 ArcCatalog 或 ArcMap 中的帮助菜单访问这些文档。文档也可在resources.arcgis.com/en/help/main/10.2/index.html找到。在 Google 中输入正确的问题也可能非常有帮助。编程论坛,如 Stack Exchange (gis.stackexchange.com/) 或 ESRI 的 GeoNet (geonet.esri.com/welcome),都是询问各种编程问题的宝贵资源。几乎每个问题都有答案(但永远不要害怕自己提问!)。

享受创造解决方案和工具的过程,并在你未来的所有地理空间编程挑战中好运!

posted @ 2025-10-01 11:27  绝不原创的飞龙  阅读(4)  评论(0)    收藏  举报