Follow Me: ArcPy(1)

引子

做地理信息系统(GIS)这么多年了,说来惭愧,也许是因为做过的项目深度不够,用ESRI产品还真不多,以至于在扎堆的ESRI拥趸中显得格格不入。只好标榜自己是开源先锋。扛起了开源这面大旗,迈开大步,自觉脚下生风。只吹的我脊梁发冷,这才发觉当年VB for ArcObject 就没学好,也是自称是.Net爱好者而搪塞过去的。

说起这个系列的引子,那还是在一次会议上遇到了一位资深的,退休的,海洋气象局的老专家。由于在淡水湖物理模型和预报方面的杰出贡献,这位大爷同时被两所公立大学返聘,既授课又做研究。也许是和课堂里的年轻人接触多了,老头看上去意气风发,攀谈起来完全像一个正当年的叫兽:

“你知道吗,小笨?我也开始学习写程序啦!”。

”呦,您这院士级别的大神也亲自写程序啊?哪能劳您大驾啊,找个学生写呗。万一累着了,那是国家的损失啊,我们不答应!“,我还是很佩服自己谄媚的能力,说出这番话脸不红心不跳,用的还不是自己的母语。不由自鸣得意起来。。。

”嘿嘿,我在学用Python!”。

这句话到叫我不爽,顿时语塞。这条大蟒蛇近些年来越盘越大,已然有了佛罗里达缅甸巨蟒的阵势。特别在科学研究领域,俨然成为大家喜闻乐见的常用工具了!在GIS领域,龙头老大ESRI也是耳根发软,从ArcGIS 10.0起,Python正式登堂入室,成为了地图/数据处理自动化(map automation)的主推标准语言。记得刚推出不久,就有一个在挂羊头卖狗肉的研究所里,不务正业做GIS的中年妇女和我炫耀过:

“你知道吗?伦家现在用ArcMap根本不用菜单和鼠标噢!分分钟命令行Python搞定哦!”

作为一个每天几个小时盯着putty工作的人 VS 一个听到putty只能联系到HomeDepot或者Lowe’s的大妈,command line这种高级词汇,我是肿么也莫想到会出自她口。脸上笑意盈盈,心中暗想:

“这有什么新鲜,老子当年在大学里熬夜赶AutoCAD图以换取机械制图免考资格时,天天敲命令行,也没觉得自己是个金领啊。”

有意无意间,这加深了我对ArcPy的误解。特别是本人对Python就不感冒,至今无法明白为嘛大家能接受一种以缩进来标识作用域的程序语言。说实话,Python我是有在用的。特别是服务器端复杂一点的脚本。 Perl我是不愿学的,bash script只用来show off, 可惜我又不懂awk,无法吸引神仙妹妹的注意。。。

近来有一个项目,由于需要Network Tracing功能,终于说服自己还是用ArcGIS Server吧。可笑的是同事们都很惊讶,他们原以为我又要自己reinvent wheel了。可惜我没有(多么痛的领悟,期待好声音第3季),原因我能告诉他们吗:世界杯开始了耶!在这个Throwball被称做Football的国度,足球只能改名叫Soccer。你能想象ESPN免费现场直播重要比赛吗?World Cup is free for streaming on WatchESPN。。。

数据交给刚招的小弟去做了,模型原型也发给小弟们去完善了。 Web App prototype is ready。偶搞了一个OpenLayers+JQuery直接访问ArcGIS Server GPService through REST API。突然感觉无事可做,抄起电话打给一个大神聊天。这位仁兄当年任职于一个10年前业界大大出名的GIS公司:皇上,您还记的黄浦江畔的。。。

”Arcpy咋样啊?“
”没啥用。复杂一点的在服务器端还要写SOE”
“哦”,我刚一放松,他又说道;
“可是现在EDN license不再对大学免费了,贵啊!从10.0起就再没用过SOE了。“。
既然大学已经不免费了,我们自然也得不到实惠的价格。ESRI自从ArcGIS 10.0起,有了太多的吐槽点。闲着也是闲着,让我来看看能不能对ArcPy进行一番深度吐槽呢?

Arcpy前世今生

Arcpy是ESRI在 2010年,随ArcGIS 10 (好吧,实际上就是9.4,命名为10.0明显高大上了)推出的Python第三方库。其实在ArcGIS 9.2中,就已经推出了arcgisscripting。只不过那个时候ESRI还比较暧昧,没有公开偏袒Python。文档中只是说提供对Python, VBScript, JScript和Perl的支持。随着微软停止了对VB,VBScript, JScript的支持,Perl让人感觉廉颇老矣,对Python的全面支持也是ESRI不得不做出的决定。事实证明,ESRI的宝押对了。让我郁闷的是:从10起,我们再也看不到那些华丽的UML
VB Programming Model了,偶们再也不能向师妹炫耀Basic programming的不简单了,Say Goodbye to VB for ArcObject。。。(假装感伤而已,本人除了在课堂上,莫有写过一行VB for ArcObject)说到ArcObject Object Model Diagram,甚是让人怀念啊。每当长夜漫漫,无心睡眠;每当思念师妹,心烦意乱;每当考试前夜,辗转反侧,只要拿出一张ESRI Object Model Diagram凝视端详一番,保证你立即昏昏欲睡。如果一张无效,OMD一共30张,总有一张适合你。

但是,Arcpy真的可以取代过去VB/VBA在我们心中的地位吗?小伙伴们能期待Arcpy for ArcObject 吗?Arcpy和ArcObject是神马关系尼?好消息是:ArcObject的中心地位从没有动摇过,而VB for ArcObject的真正继承者却是.Net(如果你认为VB 和 VB.Net是一样的,那我就只有暗笑了)。Arcpy对ArcObject来说,至多只能算是个外围女吧,因为她取代的只是VB的小姐妹:VBA在ArcMap里的Geoprocessing交互功能。从ArcGIS 9.2起,ArcObjects引入了一个新组件Geoprocessing
Object。这个新组件主要是用来自动化ArcMap中的Geoprocessing 工具箱的,包括:

Analysis toolbox

Cartography toolbox

Conversion toolbox

Data Management toolbox

Editing toolbox

Geocoding toolbox

Linear Referencing toolbox

Multidimension toolbox

Spatial Statistics toolbox

(抱歉啊,从小笨第一天用ArcGIS起,就没用过中文版的。)

在Arcpy中,所有的geoprocessing功能都是作为直接定义在Arcpy模块下的函数实现的。换句话说,在Arcpy中进行任何geoprocessing工作,只需要直接调用arcpy下直接调用函数,例如:

arcpy.env.workspace="C:\Users\micropentium6\Documents\ArcGIS\GIS Data"

arcpy.MinimumBoundingGeometry_management("Marinas.shp","Marinas_CH.shp","CONVEX_HULL","ALL")

在上面的例子中,生成点图层Marinas的凸包,只需要调用函数MinimumBoundingGeometry_management,并传递合适的参数,比如”CONVEX_HULL”。细心的读者会发现,这个函数中的参数设定和工具箱–〉数据管理–〉Minimum Bounding Geometry对话框中的设定是一一对应的。真的是自动化啊!如果你想了解更多的数据处理函数,请猛击这里。细心的读者也许已经注意到了,这些函数的命名似乎是有规律的:“工具英文名+下划线+工具箱别名“就是Python函数名。

工具名_工具箱别名

如果你不知道工具箱的别名,请右键工具箱-〉属性。

你一定开始抱怨Arcpy太没有技术含量了?她还能做什么其他的活儿吗?根据ESRI博客:

“ArcPy is a site-package that builds on (and is a successor to) the successful arcgisscripting module. Its goal is to create the cornerstone for a useful and productive way to perform geographic data analysis, data conversion, data management, and map automation
with Python.”

Arcpy是ESRI基于arcgiscripting 模块开发的Python第三方库。在传承了arcgiscripting的优秀品质的基础上,Arcpy可以高效的承担地理数据分析,数据转换,数据管理和地图自动化的工作。

Arcpy除了进行数据处理自动化外,在首作10.0中还提供了三个模块:Mapping(制图,arcpy.mapping), Spatial Analyst(空间分析,arcpy.sa),和 Geostatistical Analyst(地理统计分析,arcpy.ga)。10.1中加入了Data Access(数据访问模块,arcpy.da)和Time(时间模块,arcpy.time)。除此之外,Arcpy还定义了其他一些函数和类。看上去是不是有点儿乱啊?我们将在后续章节一一介绍这些模块,类,函数的使用。

抱怨技术含量低的读者们注意了,小笨要稍微展开一下Arcpy的底层实现,以帮助读者理解为甚Arcpy会有如此混乱的组织结构(以下内容可能会引起观众不安,请相关人士在家长带领下有选择的阅读 )。

ESR产品二次开发的基础是ArcObject。Arcobject是基于微软COM技术实现的。在ArcObject中,Geoprocessing(数据处理)组件实现的是GpDispatch组件对象(coclass)。调用数据处理函数的参数只能是字符串和调用数据处理对象返回的结果。因此,它被ESRI歧视性的称为粗粒度对象(coarse-grained object)。好吧,貌似对它不公平,那我们就叫高级对象吧。Arcpy数据处理函数作为ArcObject数据处理组件的高级实现,使高级中的高级,下里巴的低级活儿它是不干的。所以,您也不要期待Arcpy可以象VBA老婆一样上的厅堂,下的厨房啦。偶们是女神,只干”粗粒度“的高级活儿!等会儿,那Arcpy中的那些模块有是神马呢?小笨不仅笨,而且懒。不可能真的深入Arcpy实现去一探究竟,只是抄袭了Arcpy
mapping模块的帮助文件加上自己的想象YY一番:

Arcpy.mapping was built for the professional GIS analyst (as well as for developers). Traditionally, the scenarios listed above had to be done
using ArcObjects and it often proved to be a very difficult programming environment to learn for the average GIS professional. Arcpy.mapping is a courser-grained object model, meaning that the functions are designed in a way that a single arcpy.mapping function
can replace many lines of ArcObjects code.

上面大致的意思是说Arcpy.mapping是用来简化以前需要ArcObject来实现的一些功能。Arcpy.mapping一行的函数调用就可以解决以前使用ArcObject很多行才能实现的功能。以下的代码绘制mxd的图存到PDF中,真的只有两行哦!

mxd = arcpy.mapping.MapDocument("C:/Project/Watersheds.mxd")
arcpy.mapping.ExportToPDF(mxd, "C:/Project/Output/Watersheds.pdf")

不知道ArcObject的大神们看到这样的文档会不会吐血。我估计不会。因为他们实在是习惯了ESRI的一贯风格:你千辛万苦对ArcMap做的扩展,刚自鸣得意了一阵,就惊讶的发现一样的功能已经在下一个ArcMap的版本中实现了!不过,ArcObject developer们也不必过分担心,并不是所有的功能都可以用Arcpy实现的。

“Arcpy.mapping is not a replacement for ArcObjects but rather an alternative for the different scenarios it supports. ArcObjects is still necessary for finer-grain development and application customization, whereas arcpy.mapping is intended for automating
the contents of existing map documents and layer file”

ESRI向您保证:Arcpy只用来做自动化。细粒度的开发仍然需要直接调用ArcObject。

那说到底,Arcpy能直接调用ArcObject吗?俩个痣:不行!Arcpy公开的文档没有提供直接访问ArcObject的方法,但是,我们可以相信Arcpy的底层实现一定是调用ArcObject的。不信的话,小笨不才,愿意带你去看看代码。小笨座机上的Arcpy是安装在这个地方的:

C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy

在这个文件夹下,我们可以看到一个Python 库的典型结构。如果您不了解“典型结构”,请戳这里。注意到文件夹”arcobjects”没,see?神马,这不足以证明Arcpy调用ArcObject? 好吧,我们继续。。。其实,代码阅读其实就像艺术品鉴赏一样,各人的品味情趣不同,欣赏的要素,体会的意境便大相径庭。卖够了关子,该开始干活儿了。就让我们从最常用的MapDocument开始吧。

MapDocument的类定义在_mapping.py中。我们注意到文件以下划线开头,表明文件中的代码为库的内部实现。

class MapDocument(mixins.MapDocumentMethods,_BaseArcObject):

我们知道Python也是全面支持面向对象编程的。可是小P不学好啊,不学Java,不学C#,偏要学C++,搞多重继承。 MapDoucment的基类有两个: MapDocumentMethods和_BaseArcObject。后者的名字是不是再次证明了我的判断呢?神马,读者你真是不见黄河心不死啊。我们继续往下看。MapDocumentMethods是arcobjects\mixins.py中一个类 。mixins这个名字很有品味的说。不太了解Python,但有Java或者C#经历的同学们可以把Python中的mixins看作是Interface。由于上面的派生类MapDocument没有显式的定义构造函数,根据Python的MRO规则 (Method
Resolution Order, Python 2.3前用深度优先算法, 2.3后,使用C3算法,可以使用class.mro()的静态函数来查看Python MRO的结果),相应的基类(class MapDocumentMethods)构造函数被调用。

class MapDocumentMethods(object):
    def __init__(self, mxd):
        """MapDocument(mxd_path)

           Provides a reference to a map document ( .mxd ) stored on disk or to
           the map document that is currently loaded within the ArcMap
           application (using the CURRENT keyword)

             mxd_path(String):
           A string that includes the full path and file name of an existing map
           document ( .mxd ) or a string that contains the keyword CURRENT."""
        assert (os.path.isfile(mxd) or (mxd.lower() == "current")), gp.getIDMessage(89004, "Invalid MXD filename")
        super(MapDocumentMethods, self).__init__(mxd)

MapDocumentMethods的构造函数貌似直接调用了“基类”的__init__函数,如果super返回基类对象的话。然而,在多重继承这个“混乱”的世界里,如果我告诉你这里的super调用了MapDocument第二个基类_BaseArcObject的构造函数,你会觉得我在忽悠你吗?_BaseArcObject是arcobjects文件夹中_base.py中定义的一个内部类。它的构造函数十分强大,可以接受无限任何参数。为神马会是这样呢?那是为了配合前面说过的“忽悠人”的python多重继承MRO!这一切的乱象都是由super这个家伙导致的。要想了解super为毛super,Python官网推荐此篇文章:请戳这里

class _BaseArcObject(object):
    _arc_object = None
    def __init__(self, *args, **kwargs):
        """Wrapper for ArcGIS scripting Arc Objects --
           Create a new object instance if no reference is passed in."""
        super(_BaseArcObject, self).__init__()
        self._arc_object = gp._gp.CreateObject(self.__class__.__name__,
                    *((arg._arc_object if hasattr(arg, '_arc_object') else arg)
                        for arg in args))
        for attr, value in kwargs.iteritems():
            setattr(self._arc_object, attr, value)
        self._go()

我们可以看到,它偷偷的给自己定义了一个_arc_object的静态成员变量。更搞笑的是,这个CreateObject是定义在geoprocessing名空间的。说它搞笑是因为你会觉得geoprocessing会consume BaseArcObject,但现实(实现)却相反。这个叫_gp的是这样定义的:

class Geoprocessor(object):
    """Represents a geoprocessing object in ArcGIS"""
    def __init__(self):
        """Geoprocessor()"""
        self._gp = arcgisscripting.create(10.0)

哦,到这里,我想各位和小笨过不去的看官们算是了然了吧?说到底,Arcpy和ArcObject的交集是在这个arcgisscripting,也就是前面介绍的Arcpy的前辈。值得注意的是,这个arcgisscripting库是不在Arcpy文件夹里的。在小笨的座机上,它住在:C:\Program Files (x86)\ArcGIS\Desktop10.2\bin,也就是arcmap.exe他们家。而且它不是一个文本文件,它的全名是:arcgisscripting.pyd。其实它就是一个DLL。我想我们的代码追踪就到此为止吧。还没有看懂的朋友们,小笨再在为你们梳理一下:

Arcpy中的各个模块,无论早期的geoprocssing还是新加入的mapping,内部统统都是通过调用arcgisscripting中提供的ArcObject对象来实现各种功能的。arcgisscripting.pyd中一定实现了由COM到Python对象的转换(这种转换也许不能完全适合Arcpy,因此在Arcpy中小笨也找到了arcobjectconversion.py这样的文件)。但无论如何,Arcpy的底层实现是完全依托在arcgiscripting上的,并且由于历史的原因,所有的arcpy模块都会依赖geoprocessing模块中的部分函数来实现对arcgisscripting的访问(小笨觉得这种依赖完全可以被重写)。既然我们谈到了COM访问,下面有必要讨论一下,Python访问COM的性能问题。

我们都知道ESRI是微软的脑残粉,COM的忠实追随者 (SDE and Geodatabase for Uinx/Linux除外,那是因为Oracle)。如果不使用VC++,与COM交互是要付出代价的,需要interoperability。从VB/VBA迁移到Python的性能损失会有多大呢?没能找到最新的数据,又不愿意自己测试然后让读者骂。只找到这个:

CSDN Image Caching Service Is Down, Shame on you, CSDN!

貌似损失不小哦。注意测试环境是ArcGIS 9.3,四年过去了,不知是否有改进。不过,这个貌似没法记在ESRI头上吧?MS才是元凶!看着用户凶凶的表情,ESRI只能怨怨的说:”皇上,臣妾做不到啊!“

小笨猜测,ESRI会逐渐把ArcObject中和COM松耦合的,对性能要求不高的组件逐渐迁移到使用ArcPy访问,而其他迁移性能损失比较大的仍然使用传统方式访问。这个过程会持续几个版本。我们有理由看到Arcpy的模块会越来越多。但也不排除Arcpy在将来会被放弃,看看Web ADF for .net and java就明白了。其实ESRI也有苦衷啊,自己说到底只是个应用开发商,对上游的平台,语言,协议没有话语权。所以只好跟随大趋势变来变去。

在了解了Arcpy是什么,能干什么之后,从下一章开始,小笨会带领大家逐渐熟悉Arcpy的方访面面。So,follow me!

转载自:https://blog.csdn.net/micropentium6/article/details/37686557

You may also like...