[初学入门]arcgis中python脚本学习.pdf

《[初学入门]arcgis中python脚本学习.pdf》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

[初学入门]ArcGIS中Python脚本的使用By：飞天小猪目录写在前面的话....................................................................................................................................2前言...................................................................................................................................................2一、PYTHON语言基础.....................................................................................................................31数学运算符.......................................................................................................................................32字符串操作.......................................................................................................................................43模块的使用（MODULES）.................................................................................................................54使用DEF构建函数............................................................................................................................65流程控制结构：IF，WHILE，FOR......................................................................................................76简单输入和输出...............................................................................................................................9二、ARCGIS&PYTHON.....................................................................................................................101如何创建地理处理对象（GEOPROCESSOROBJECT）..........................................................................102获取地理处理帮助.........................................................................................................................112.1举例：如何使用GeoprocessorProgrammingModel中的Lists............................................113使用地理处理工具——TOOLBOXES和ALIASES.................................................................................124在建模中使用脚本（SCRIPTSINMODELBUILDER）............................................................................135在PYTHONWIN里调试地理处理脚本.............................................................................................195.1调试选择和消息.....................................................................................................................205.2PythonWin的调试工具............................................................................................................215.3地理处理工具举例..................................................................................................................226使用描述（DESCRIBE）和存在（EXISTS）获取数据信息...............................................................226.1描述..........................................................................................................................................236.2存在（Exists）.........................................................................................................................246.3在循环中使用描述和存在......................................................................................................247在PYTHON脚本中使用地图代数（MAPALGEBRA）........................................................................278数据管理和指针（DATAMANAGEMENTANDCURSORS）.....................................................................288.1数据管理（DataManagement）...........................................................................................288.2指针（Cursors）......................................................................................................................29附录1：地理处理脚本中输入&输出方法指南..............................................................................32附录2：其他...................................................................................................................................331/33 写在前面的话一直想学习ArcGIS中的Python脚本，大四下半学期终于有了时间，可是想找到这么一本好的教材不容易。茫茫互联网，终于找到了旧金山州立大学JerryDavis教授的个人主页，对其中《GeoprocessingScriptsWithPython》如获至宝，独乐乐不如众乐乐，现在将其教程翻译并结合自己的学习情况给出总结。希望能够给更多想学习Python的同学一个参考。另外，在我刚开始接触Python时，是看了台湾辅仁大学一位老师的视频课件，在此致谢。我想从两个大部分总结：一、Python语言基础；二、ArcGIS&Python。其中第一部分参考了《Python精要参考（第二版）》、《Python编程金典（读书笔记）》等书籍文献。对于多数读者来说，可能或多或少有一些编程基础，所以理解起来应该不成问题。文中多数数据来自JerryDavis教授的主页，放在“C:prog”目录下，为了直观，我将运算结果一并编辑，方便参考。值得一提的是ArcGIS的在线帮助文档，一个实时更新的GIS宝库，很多专业性知识都可以找到答案，点击链接ArcGIS10中文帮助、ArcGIS9.3.1或9.3英文帮助。获取更过脚本例子来学习：ESRI的地理处理模型和脚本工具库是个不错的选择。由于我也是初次接触，翻译或者心得难免有纰漏之处，希望同仁们可以多多交流！前言在GIS建模或GIS数据管理中，你可能经常需要处理一系列步骤才可以完成的工作；你可能有一个工作目录下的数据需要重投影、裁剪到研究区域，或者用某种方法组合成期望的结果；我们也经常需要根据不同情形用不同方法处理数据，因此我们需要作出选择，而高质量的决策需要考虑很多低水平的决策，这可以通过脚本程序模型辅助完成。脚本编程的主要目的是使枯燥的处理数据工作自动化，通过逻辑来指挥处理过程。我想自动化和逻辑是关键，它们区别于我们多数使用计算机时的交互活动。我们发E-mail，写文章或者设计地图，都需要和计算机交互，而处理一系列数据，我们需要自动化和利用逻辑来指导自动化。在地理处理脚本逻辑中，我们需要在允许我们做的事情中作出决定，比如，处理栅格数据不同于矢量数据，或为没投影的数据设置投影，或处理仅在特定时间搜集的数据集。对于重要的GIS工作来说，脚本以及其他形式的程序是必需的，而非可有可无。在接下来的联系中，我们会探索Python的使用以及创建脚本来使用ArcGIS里众多的地理处理工具。所有你能在ArcToolbox或Model中使用的工具都能够用在Python脚本中，这些脚本可以生成脚本工具，像其他地理处理工具一样使用。2/33 一、Python语言基础安装PythonWin，在…ArcGisDesktop9.3.isoDesktopPythonWin目录下可以找到PythonWin的安装程序，默认是不安装的，。同时会安装win32com以及允许任何脚本在基于Dispatch的地理处理过程中工作。ArcGIS10中引入了全新的PythonWindow来增强内嵌的Python体验。警告：不要尝试更新随ArcGIS安装的Python到一个新的版本！下面介绍Python的一些简单语法和规则。1数学运算符Python提供了多样化的通用数学运算符——多数编程语言的特征，以及许多通过import的modules提供的符号。常用的有+,-,*,/,**(幂),%(取模，即除后的余数)。下面的表格显示了整型（Integer）和浮点型（Float）各种组合运算的结果，记住一条规则，只要参与运算的有浮点型，则结果为浮点型；全为整型时，结果才为整型。3/33 输入表达式结果Notes2+35整型结果2.+35.02.是浮点型，结果浮点型2-3-12*36整型结果2.*36.0浮点型5/22整型5./22.55%21取模Az=27090取模的用途之一——方位角加180后逆转方向Newaz=az+180Printnewaz%3605**22525**0.55.0没有sqrt()功能，除非添加math模块2字符串操作注：使用Python帮助：有超过30种内置方法来处理字符，请到SequenceTypes下的StringMethods寻找帮助！字符串是一串字母，比如’SanFrancisco’，字符串下标从0开始。学习字符串语法的最好方法是自己动手尝试，下标展示之：输入结果Notesprint'zhulj'.capitalize()s='zhulj'Zhuljs.capitalize()即将capitalize()方法用于sprints.capitalize()prints[0]zs1=s[1]hStrings可以像一个字母列表一样处理，第一个字prints1母下标为0，某个字符段可以用1:3来格式化：从prints[-2:]lj第1个的开头到第3个的开头，不包括下标为3prints[2:3]u的字母；s[-1]表示倒过来第一个，相当于s[len(s)-1]prints[2:4]ulprints[2:],s[:5]uljzhuljs2=s.upper()ZHULJ我们可以将字符串方法的结果赋给新的变量prints2s3=s+s2zhuljZHULJ字符串组合用“+”prints3prints*3zhuljzhuljzhulj字符串重复用“*”，后为重复次数selstr='"elev">1000'"elev">1000printselstr字符串可以使用单引号或双引号，跨行时用双引othersel=”’elev’>1000”号。‘elev’>1000printotherselprints.isupper()False一些方法返回值为布尔型（True或False），一些prints2.isupper()True返回索引值（下标值）4/33 p='d:/work/lu.shp'10printp.find(‘.’)printp.find(‘/’)2plist=p.split('/')['d:','work',你可以用split()方法解析出不同的字符串片段，并printplist'lu.shp']创建一个列表（List），我们可以使用其中不同的printplist[0]d:元素printplist[1]workp2='d:\work\soil.shp'd:worksoil.shpprintp2反斜线“”和某些字母一起有特殊用法，如 print'Jerry'sKids'Jerry'sKids为换行，“”为转义字符，如“\”则表示“”Jerry'sprint'Jerry's Kids'Kidsp3=r'd:worksoil.shp'字符串前加“r”则强制“”代表其本身，而非d:worksoil.shpprintp3转义字符，这对于文件路径的操作很方便3模块的使用（Modules）Python提供了一系列内置的方法（大量依赖于模块）用于通用编程。Python安装时自带了大量Modules，最常用的有math,sys,random,array以及os.path。当然还有好多Modules可以下载，比如数字处理（Numeric）——numpy，可在www.python.org或www.google.com里搜索。www.python.org/moin/NumericAndScientific页面中列举了一些。使用Module前，必须import之。通常我们会将一行import放在程序顶部，比如：importarcgisscripting当然，这不必成为你程序的第一行，但必须在使用它里面方法之前。当要引用多个模块是，中间用逗号分隔，比如：importarcgisscripting,sys,string,os,math我们也可以自己为频繁使用的方法创建Module，下面，我们开始体验内置的Modules。math和random模块很多常用的数学计算功能都可以通过math找到，比如三角计算或对数计算，如果要使用复杂数字，就使用cmath模块。和之前一样，通过以下表格来体现模块的使用：输入结果Notesimportmath2.0以10为底的对数printmath.log10(100)printmath.log(100)4.60517018599自然对数printmath.pi3.14159265359π是一个静态常量，所以不需要括弧5/33 pi=math.pi如果不想总是输入“math.pi”可以将3.14159265359printpi其赋给一变量pi3.1415926535897931不需要print即可查看变量值printmath.sin(radians)printmath.cos(radians)三角函数的计算是弧度制printmath.tan(radians)degrad=pi/1800.78539816339744828度转化为弧度45*degradsin=math.sin0.70710678118654746即使功能函数（像sin）都可以赋给一sin(45*degrad)1.0个变量sin(90*degrad)math.e2.7182818284590451math.hypot(3,4)5.0此方法是求三角形的斜边x1=520382;y1=4152373不同赋值语句间用“；”分隔x2=520475;y2=4152963xr=x2-x1yr=y2-y1math.hypot(xr,yr)597.28468923956189不同的方式，相同的结果math.sqrt(xr**2+yr**2)(xr*xr+yr*yr)**0.5importrandom0.27281588185756478random.random()rnd=random.random0.4456393835072503rnd()random()方法，每次结果都不同，值mu=50域为[0.0,1.0)s=1046.528202194printrandom.gauss(mu,s)4使用def构建函数有点像Module，但更简单，函数是一个自己定义功能，用在之后的代码中，并且提供任何你想要使用的参数。这个函数从此可像变量那样在程序中使用，结合例子更容易理解。接下来的代码定义了一个将度转换为弧度的简单函数，同时也定义了一个弧度转换为度的函数，它们和Excel内置的函数类似。importmathdefradians(angdeg):returnangdeg*math.pi/180defdegrees(angrad):returnangrad*180/math.piprintmath.sin(radians(45))printdegrees(math.acos(0.5))运行之，得到结果：0.70710678118760.06/33 5流程控制结构：If，While，For任何脚本或编程语言的一个重要特征就是执行一系列不同情形语句的能力。你想要创建一系列山影栅格来代表夏天、冬天和春秋分。山影（hillshade）工具需要有太阳高度角和方位角作为输入参数。重要日期太阳倾角夏至（6月21日）23.44春秋分（3月21日，9月21日）0冬至（12月21日）-23.44接下来是一段相当简单的代码，通过太阳倾角（太阳光线正午垂直照射的纬度）获取太阳角和方位角以及纬度。输入两个参数：lat（研究区域的纬度，南半球为负）和decl（太阳倾角），由此得到sunangle和azimuth：lat=30decl=20sunangle=90-lat+declazimuth=180ifsunangle>90:sunangle=180-sunangleazimuth=0printsunangle,azimuth上面的例子中lat和decl强制赋了值。有三种流程控制操作：if仅在一个特定情形下才执行语句；while当一种情形存在下，持续执行语句for遍历一系列值这些语法和def有些相似：初始语句后加顿号、需要执行的语句块有缩进。这三个结构的一些重要的公共特征：①if、while、for语句均以冒号结尾，接下来是缩进的代码块，用于if、while、for定义的情形。在脚本编写窗口，你会发现，你在一行末尾打上冒号后，下一行自动缩进，在接下来的一行按下退格键取消缩进。②如果你只需做一件事情，你可以在冒号后面同一行添加简短的语句，比如：ifx>0:print‘x比0大’print‘下一行不要缩进了。。。’if（continued）接下来，我们会探索一下另一个方便的模块：os.path:开始之前，在d:/下创建一个“testfolder”文件夹，然后新建一个“test.txt”文件；尝试以下代码段，确保print语句前有缩进。importos.pathifos.path.exists("d:/testfolder/test.txt"):print"测试文件夹存在"print"txt文件存在"7/33 elifos.path.exists("d:/testfolder"):print"测试文件夹存在"print"测试文件夹存在，但txt文件不存在"else:print"两者都不存在"可选探索示例接下来的例子做的事情对GIS非常重要，但是实际上不用任何地理处理代码。USGS7.5米分辨率DEM（数字高程模型）是文本文件（USGSDEM文件），投影为UTM，UTM北向和东向单位是米，但是高程单位可能是英尺（feet）或米（meters）。因此在获取垂直或水平距离信息时会有问题，比如坡度可以通过垂直距离/水平距离获得。如果你不在使用Z值之前设置为0.3048，将会出现错误结果。但是不幸的是，你可能不知道DEM文本文件的垂直单位是英尺还是米。这些信息保存在第539个字符里，“1”代表英尺，“2”代表米，所以可以通过读取这个文件判断。下面的脚本演示了上述内容：importfileinputinfile=r"c:progpendatawoodside.dem"firstline=fileinput.input(infile)[0]unitchar=firstline[539]unit="(unknown:nota7.5'DEM?)"ifunitchar=="1":unit="feet"ifunitchar=="2":unit="meters"print" Elevationin"+""+unitfileinput.close()输出结果：Elevationinfeetwhile（continued）运行下面的代码，说明了一种while循环：x=1whilex<10:printxx=x+1屏幕依次输出1~9下面说明一下“==”（等于）的概念：x=5x=5z=x==4z=x==5printzprintz输出False输出True“==”是逻辑运算符之一，其他有“<”（小于）、“>”（大于）、“>=”（大于或等于）、“<=”（小于或等于）、“<>”（不等于）。使用逻辑运算符计算得到结果为布尔型：true（1）和false（0）。下面例子简单体会一下布尔型表达式：x=1whilex<10:8/33 printxx=x+1表达式“x<10”结果是true或false，所以这样允许我们在计算完一种情况时运行一系列代码。许多情况下我们需要使用条件代码。while循环的一个优点是允许我们跳过整个部分，如果条件不满足初始情况。由于while提供一种容易结束循环的方法，我们甚至用它代替if语句。当循环一个数据集时（GIS中很常用的工作）while循环很有用。在后面地理处理中我们会接触一些例子。for尝试下面代码，演示了for循环：forxin[1,2,3,4]:（注：[1,2,3,4]用range(1,5)代替是一样的）msg="helloworld"printstr(x)+""+msg（注：当我们希望一个数字x和字符串连接时，必须先对数字进行格式转换即str(x)，否则系统报错）下面的代码创建并输出指定文件夹内shp文件名列表（每个都以‘.shp’结尾）importosws="c:/prog/hmbarea"ilist=os.listdir(ws)#创建一个列表保存工作文件夹内的文件fcs=[]#创建一个空列表，保存结尾为‘.shp’的文件foriinilist:ifi.endswith(".shp"):fcs.append(i)forfcinfcs:printfc（输出结果如右图所示）下面这个例子的循环较多次数。变量mu是算术平均数，s是标准差——这两个是random.gauss()用到的参数，可以尝试改变n的值查看结果！importrandommu=50s=10z5=mu-s*1.96z95=mu+s*1.96count=0n=100000foriinrange(n):x=random.gauss(mu,s)ifxz95:count=count+1printfloat(count)/n（每次运行的结果都不同，但都在0.05左右，即统计结果在5%左右）6简单输入和输出我们现在利用前面计算太阳角代码的片段，之前是直接指定参数值，现在我们有很多种9/33 方法提供输入参数，现在我们用sys方法。尝试下面的代码，点击（run运行）之后，在对话框中函数自变量（Arguments）一栏填入：4023.44，如下图：importsyslat=float(sys.argv[1])decl=float(sys.argv[2])#使用arguments（argv）方法从“Arguments”一栏中获取输入参数，并指定一个浮点型转换将字符型输入值传递给lat和declsunangle=90-lat+declazimuth=180ifsunangle>90:sunangle=180-sunangleazimuth=0print"正午太阳角="+str(sunangle)print"方位角="+str(azimuth)（结果：正午太阳角=73.44方位角=180）二、ArcGIS&Python1如何创建地理处理对象（geoprocessorobject）所有geoprocessing的Python脚本都可以通过importarcgisscripting或者win32com去穿件geoprocessorobject。下面的例子显示二者区别：arcgisscriptingmodule需要在Python2.5.1版本上创建并且需要此版本创建geoprocessor；通过win32com创建的geoprocessor可以在不同的Python版本上运行。#9.3importarcgisscriptinggp=arcgisscripting.create(9.3)gp.workspace=”c:/Tongass”gp.clip_analysis(“standb4”,”clipcov”,”standb4_clip”,”POLY”.”1.25”)#Dispatchimportwin32com.clientgp=win32com.client.Dispatch(“esriGeoprocessing.GpDispatch.1”)gp.workspace=”c:/Tongass”gp.clip_analysis(“standb4”,”clipcov”,”standb4_clip”,”POLY”.”1.25”)理解ArcGIS中数据多样性格式对我们理解地理处理工具很有帮助。比如，特征数据可能是①单个shp文件；②geodatabase（地理数据库，我们可能指定地理数据库为10/33 工作空间）；③多边形、弧或点要素的coverage数据。当我们想遍历工作空间里的coverage文件时，应使用ListDatasets而不是ListFeatureClasses。2获取地理处理帮助如果你基本熟悉了Python的语法，便可以开始熟悉ArcGIS的Geoprocessing啦，获取这些方法帮助的途径有两个：①ArcGIS帮助系统，GoToGeoprocessing/Automatingyourworkwithscripts/Scriptingobject:ProperitesandMethods.这里你会看到GeoprocessorObject,这个能让你很方便地获得所有对你有用的条目、设置和其他操作文档。比如，你想得到特定类型的文件，就找到listfeatureclasses方法：fcList=gp.ListFeatureClasses(“w*”,”polygon”)注：此方法的语法为：object.ListFeatureClasses({wildCard}AsString,{FeatureType}AsString)AsPythonList{}=optionalwildcard为通配符，和星号（*）组合使用，如果没有通配符则返回工作目录下的所有featureclasses。②GeoprocessorProgrammingModel（PDF），包含方法（左箭头表示）、属性（可读写的表示为杠铃形，只读的表示为部分杠铃形）2.1举例：如何使用GeoprocessorProgrammingModel中的ListsLists（列表）及其属性和方法在图表中用紫色标出，如下：现在我们试着编写一段脚本列举出属性表中所有属性值（Fields）（以hmbarea栅格土地利用为例，文件存在c:proghmbarea下）importarcgisscripting,sysgp=arcgisscripting.create(9.3)gp.workspace="c:/prog/hmbarea"fieldList=gp.ListFields("landuse","*","all")dsc=gp.describe("landuse")print"landuse"+""+dsc.DataType+":"forfieldinfieldList:...printfield.Name,field.Type(此时输出结果如右图)11/33 3使用地理处理工具——Toolboxes和Aliases总所周知，地理处理工具在脚本中的使用和ArcToolbox中相同，但是需要提供工具名称来使用它们。但是记住一个名称可能有好几个工具，比如，裁剪工具（Clip）在Coverage->Analysis->Extract工具集里，另一个是在DataManagementTools下的Raster工具集下。区别是每个工具适用不同的数据类型，但是我们如何在脚本中区分它们呢？在ArcToolbox中，我们可以随意选择要使用的工具，但在脚本里就必须在某些方面加以区分。因此我们使用Aliases（别名）——已经成为工具名称的一部分。每一个工具都有自己的别名，我们可以通过右键->属性来查看：AliaseToolbox“conversion”Conversion“3d”3DAnalyst“geocoding”Geocoding“analysis”Analysis“ga”GeostatisticalAnalyst“arc”Coverage“lr”LinearReferencing“management”DataManagement“sa”SpatialAnalyst“cartography”Cartography“stats”SpatialStatistics现在我们用一段脚本来解释：importarcgisscripting,sysgp=arcgisscripting.create(9.3)gp.Workspace=”c:/prog/hmbarea”gp.overwriteoutput=1#OverWriteOutput:Boolean,为1表示允许覆盖已存在文件#将streams/arc转换为shp文件gp.copyfeatures_management("streams/arc","streams.shp")#利用转换后的shp文件，做200米的缓冲gp.buffer_analysis("streams.shp","stbuff200.shp",200)#用做过缓冲的shp裁剪gp.Clip_analysis("geol.shp","stbuff200.shp","geolstr200.shp")注：上面脚本用“#”注释的中文内容不要出现在脚本文件中，否则会出现错误结果截图：12/33 如果你一次使用一个工具集中的若干工具，可以通过环境设置省下一些文字：gp.Toolbox="Analysis"gp.Buffer("streams.shp","stbuff200.shp",200)gp.Clip("geol.shp","stbuff200.shp","geolstr200.shp")4在建模中使用脚本（ScriptsinModelBuilder）首先，需要记住的很重要的一点是，ArcToolbox里相当数量的工具实际上都是脚本。脚本都有一个图标。比如，空间统计分析工具（SpatialStaisticstools）几乎都是Python脚本（一些是模型中使用了脚本）。实际上，我们可以复制并编辑这些脚本作为其他用途。为了在ModelBuilder中使用脚本或将脚本当做ArcToolbox中工具使用，我们需要考虑如何给它输入值以及让其设置输出值。仍然以太阳角计算代码为例，我们给其加上两句引用，四句输入输出语句，就可以用作Modelbuilder中的一个步骤了。编辑下面的脚本，不过不要再PythonWin中运行之，因为getparameterastext和setparameterastext只能用在脚本工具（ArcToolbox或Modelbuilder）中。importarcgisscriptinggp=arcgisscripting.create(9.3)13/33 lat=float(gp.getparameterastext(0))decl=float(gp.getparameterastext(1))sunangle=90-lat+declazimuth=180ifsunangle>90:sunangle=180-sunangleazimuth=0gp.setparameterastext(2,str(sunangle))gp.setparameterastext(3,str(azimuth))这段代码中的“新面孔”：getparameterastext：是Python传递给geoprocessor（我们称之为gp）的一个方法，允许工具提供输入参数。每次你运行这个工具时，都会看到一个对话框，提示输入参数，这个方法允许你在接下来的程序中使用。索引0和1指第一个和第二个参数。setparameterastext：和getparameterastext相反，它传递第二个条目（比如str(sunangle)的值）给指定的输出参数。前两个参数索引为0和1，所以接下来输出参数的索引为2和3。然后，我们需要将脚本加进工具（Makingascriptintoatool），那样才能在ArcToolbox或ModelBuilder或Commandline中使用。如前面提到的那样，这个脚本只能用于工具，包括输入/输出方法是PythonWin不能处理的，但这些是多数工具必需的。在ArcCatalog里，定位到Python脚本保存的文件夹，最好和数据在一个盘里，右键->新建toolbox。当然也可以使用之前创建的toolbox。在ArcToolbox里，右键toolbox，选添加->scripts，填写如下图：14/33 注意：脚本文件是一个脚本工具的参考！非常重要的一点，你使用脚本创建了一个脚本工具，但是脚本本身并没有和工具一起保存，脚本工具作为toolbox的一部分保存在“*.tbx”文件中。你也可以右键脚本工具点击“编辑”，出现PythonWin或其他任何IDE窗口，这看起来好像是脚本存在于工具中。所以，记得移动时将脚本工具文件和脚本本身一起拷贝。“下一步”后是参数配置页面，如下图设置各参数如表格所示：DisplayNameDatatypeDirectionDefaultLatitudeDoubleinput0DeclinationDoubleinput0SunAngleDoubleoutput35AzimuthDoubleoutput30015/33 现在可以运行此toolbox了，对话框提示你输入latitude和declination。工具是否正确运行呢？如果是的话，你应该可以看到“Completed”，你可能会看到有一黑色窗体一闪而过，不用担心。那么，它是干什么的呢？还记得结果是输出两个数字参数，那么，这些数字哪去了呢？很好的问题，这仅能说明你能创建一个工具，但是不能想ArcToolbox那样运行。比如输出一种数据，栅格或特征数据（.shp）之类的。但是它能在作为Modelbuilder工具正常运行，下面我们将使用它的输出参数创建一个hillshade。在你的toolbox中新建一个model，将刚才创建的脚本工具（scripttool）拖进来。双击工具，输入参数，初始化之。打开“SunAngle”和“Azimuth”发现它们还是默认值，说明此脚本工具还没有运行。右键单击工具，选择Run，然后发现两个输出参数已经改变！16/33 需要注意的是：latitude范围是-90~90，declination范围是-23.44~23.44。尝试输入latitude-70，declination23.5，你会得到什么？为什么？确保你已经得到值域内的太阳角和方位角，它们将是构建hillshade的输入参数。首先，添加hillshade工具，双击指定一个elevation栅格数据（这里我选择了marbles文件夹下的elevation），用下拉条指定azimuth和altitude值为azimuth和sunangle。运行，然后右键单击输出文件，选“AddtoDisplay”在ArcMap里查看结果。探索1：如果你想通过日期获取太阳倾角，该怎么做呢？尝试下面代码，保存为getnoonsunfromdata.py，现在有五个参数，如下表进行正确设置：DisplayNameDatatypeTypeDirectionDefaultMonthLongRequiredinput1DateLongRequiredinput1LatitudeDoubleRequiredinput0SunAngleDoubleDerivedoutput35AzimuthDoubleDerivedoutput300importwin32com.client,sys,mathgp=win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")#构建两个函数，首先判断输入月/日的合法性，然后获取是一年中的第几天defjdate(im,id):#通过年月日起返回一年中的第几天lastD=[31,28,31,30,31,30,31,31,30,31,30,31]jd=[0,31,59,90,120,151,181,212,243,273,304,334]if((im>0)and(im<13)):if((id>0)and(id<=lastD[im-1])):returnjd[im-1]+idelse:print"Datenotinmonth"return0else:print"Monthshouldbebetween1and12inclusive"return0#计算倾角的函数defdeclin(day):17/33 rad=math.pi/180xlong=279.164+0.985647*dayxanom=356.381+0.985600*dayg=xanom*radelong=xlong+1.915*math.sin(g)+0.02*math.sin(2*g)elong=elong*radoblecl=23.44*radreturnmath.asin(math.sin(oblecl)*math.sin(elong))/rad#主程序month=int(sys.argv[1])date=int(sys.argv[2])lat=float(sys.argv[3])decl=declin(jdate(month,date))sunangle=90-lat+declazimuth=180ifsunangle>90:sunangle=180-sunangleazimuth=0gp.setparameterastext(3,str(sunangle))gp.setparameterastext(4,str(azimuth))我们在思考远一点，如何根据一天内的任何时间而不是正午时间，获取太阳角和方位角呢？感兴趣的可以从本文开头的链接中下砸相应代码学习。探索2：如何在PythonWin里运行这个脚本？首先我们得明确几个点：①我们将把hillshade作为脚本的一部分使用，并为其提供输入参数：一个高程栅格（elevationraster）；②GetParameterastext仅在用作工具时起作用。下表为脚本工具的参数设置：DisplaynameDatatypeTypeDirectionDefaultMonthLongRequiredInputDateLongRequiredInputLatitudeDoubleRequiredInputWorkspaceWorkspaceRequiredInputElevationRasterDatasetRequiredInputHillshadeStringRequiredInput代码如下：importwin32com.client,sys,mathgp=win32com.client.Dispatch("esriGeoprocessing.GPDispatch.1")#函数定义，同上，下面仅给出函数名称：defjdate(im,id):defdeclin(day):#主程序，使用sys.argv[]代替getparameterastext()month=int(sys.argv[1])18/33 date=int(sys.argv[2])lat=float(sys.argv[3])gp.Workspace=sys.argv[4]#输入时注意，路径应为反斜杠“”elev=sys.argv[5]hillsh=sys.argv[6]#给输出hillshade指定文件名decl=declin(jdate(month,date))sunangle=90-lat+declazimuth=180ifsunangle>90:sunangle=180-sunangleazimuth=0gp.addmessage("abouttotry")try:gp.OverwriteOutput=1gp.addmessage("afteroverwriteoutputssetting,"+gp.workspace+"/"+hillsh)gp.CheckOutExtension("spatial")gp.addmessage(gp.workspace+"/"+hillsh)gp.HillShade_sa(elev,gp.workspace+"/"+str(hillsh),azimuth,sunangle)gp.addmessage("donewithhillshade")gp.CheckInExtension("spatial")except:##printgp.getmessages()gp.addmessage(gp.getmessages())gp.CheckInExtension("spatial")阅读代码发现，第一个输入参数不适用的getparameterastext(0)，而是sys.argv[1]，这是因为getparameterastext()方法只在工具中起作用，而sys.argv[]有同样的效果，索引从1而不是0开始，当然，这需要首先引用sys模块。这里我们直接指定输出文件在输入数据文件夹内，省去了setparameterastext()方法，当然这个方法在PythonWin中也无法运行。尝试输入参数如下图，得到右下结果：5在PythonWin里调试地理处理脚本既然我们已经认真地学会了从Python中创建并运行地理统计工具，那么现在需要考虑如何调试我们的程序了。我们需要经常在Python和添加的地理处理系统引用之间调试程序。当一个地理处理工具运行失败后，我们需要从Pythonwin中得到一个丰富的消息，而不是“未知错误”。19/33 5.1调试选择和消息Python优于AML的优点之一是它有更好的调试方法，调试程序有很多选择，这里不能一一列举。打印语句（Printstatements）一开始就养成良好的调试方法是：将变量的当前值或脚本的处理过程打印在屏幕上。比如，对之前的一段脚本加以修改：importarcgisscripting,sysgp=arcgisscripting.create(9.3)gp.Workspace="c:/prog/hmbarea"gp.overwriteoutput=1gp.copyfeatures_management("streams/arc","streams.shp")gp.Toolbox="Analysis"gp.buffer("streams.shp","stbuff200.shp",200)print"FinishedBuffer.NowCliping..."gp.Clip("geol.shp","stbuff200.shp","geolstr200.shp")可以看出成功运行脚本！然而当在工具中运行时，print语句不会产生错误，但也不会输出任何东西，因此，我们用gp.addmessage("FinishedBuffer.NowCliping...")代替print语句。那么，如果想无论在工具中或Pythonwin中都可以显示消息，就可以这两句都写上。我喜欢的做法是定义一个‘sendmsg’函数来输出消息：defsendmsg(msg):printmsggp.addmessage(msg)….sendmsg("FinishedBuffer.NowCliping...")获取工具消息和（try:…except:）上面的方法很有用，但当我们运行我们并不了解很多信息的地理处理工具时就显得无能为力了。如果在Pythonwin中运行Buffer时出现错误，比如输入文件不存在等，只能看到“未知错误”的提示，这就引出了GetMessage上的方便！现在我们要体验两种调试技巧：①GetMessages地理处理方法；②Python语言的错误处理程序：try…except。①现在我们添加GetMessages查看错误信息：还是上面的代码，把“streams.shp”修改为“stream.shp”，查看错误信息：gp.Toolbox="Analysis"gp.buffer("stream.shp","stbuff200.shp",200)gp.Clip("geol.shp","stbuff200.shp","geolstr200.shp")运行之，查看错误信息，然后修改代码如下：try:gp.Toolbox="Analysis"20/33 gp.buffer("stream.shp","stbuff200.shp",200)gp.Clip("geol.shp","stbuff200.shp","geolstr200.shp")except:printgp.GetMessages()GetMessages():所以信息运行可以看到以下错误提示：GetMessages(0):只显示消息顺便说一下GetMessages的集中形式和代表的含GetMessages(1):只显示警告义GetMessages(2):只显示错误②运行下面简单的代码，并改变x的值为非零，查看结果：x=0.y=15.try:printy/xexcept:print"错误！0不能是被除数！"5.2PythonWin的调试工具PythonWin提供了一些工具：①单步执行代码②插入断点③观测变量，或其他事情。这些都可以在PythonWin的帮助系统里找到相关示例教程。这里我们演示单步执行代码的过程。以前面用过的太阳角计算脚本为例，点击运行之后，在Debugging下拉框里看到调试选项。选择“Step-throughinthedebugger”。然后，一系列调试工具出现；一个黄色三角出现在第一行代码前。按下（Stepover）按钮进入下一行，继续执行；点击（Watch）按钮，出现，在Expression下里输入lat等可以查看当前变量值；当你遍历完程序，点击关闭调试；尝试添加一个断点，当定位到某一行时，在点，然后运行此脚本在“Runinthedebugger”下，发现没有任何影响，但是修改第一行为“lat=10”，运21/33 行发现在断点处停止，，点击按钮继续程序。得到结果“Noonsunangle=76.56Azimuth=0”。5.3地理处理工具举例转换脚本（conversionscript）既然我们已经会在帮助系统里寻找答案，也掌握了几种程序调试方法，现在让我们建立一个脚本。①添加引用importarcgisscripting②创建“sendmsg”函数defsendmsg(msg):printmsggp.addmessage(msg)③gp.OverwriteOutput=1④使用DataManagementtoolbox中的Workspace工具集中的CreateFolder工具在“c:/prog”中创建新文件夹“woodside”gp.CreateFolder_management("c:prog","woodside")⑤设置gp的工作文件夹为woodsidegp.workspace="c:progwoodside"⑥使用Conversion里的ToRaster工具集中的DEMtoRaster工具将c:/prog/pendata下的woodside.dem转换为woodelev栅格数据，不要提供任何可选参数值gp.DEMtoRaster_conversion('c:progpendatawoodside.dem','woodelev')⑦在第一行前添加try:，最后一行后添加except:，最后加上，sendmsg(gp.GetMessages()).要素转换为栅格（featuretorasterconversion）相信有了前面的联系，本例很简单了，将d:/prog/pendata/landuse(polygon)根据“LU-CODE”字段转换为“lugrid”栅格数据，分辨率为30。try:importarcgisscriptingdefsendmsg(msg):printmsggp.addmessage(msg)gp.OverwriteOutput=1gp.workspace="c:progpendata"gp.FeatureToRaster_conversion('landuse/polygon','LU-CODE','lugrid','30')except:sendmsg(gp.messages())6使用描述（Describe）和存在（Exists）获取数据信息有很多情况下我们需要使用某个GIS数据的特征去处理其他数据。比如，在栅格运算中，22/33 我们可能想用一个数据的边界去界定另个数据集，很像裁剪操作，或者检查一个数据的拓扑错误。6.1描述gp.Describe方法适用于生成一个数据集的若干属性。GeoprocessorProgrammingModel中的绿色区域部分，有按等级划分的很多属性：所有数据都有一个数据类型（DataType）和目录路径（CatalogPath）。FeatureClass属性包括所有表个属性和数据集属性；栅格数据属性包括栅格波段属性和数据集属性，所以数据集属性在一定程度上是共有的（FeatureClass、Rasters、Coverages）。利用描述系统用一个数据边界裁剪另一个栅格数据，将裁剪后的土地利用图存在woodside文件夹下。importarcgisscriptinggp=arcgisscripting.create(9.3)gp.workspace="c:progwoodside"dscRas=gp.Describe("woodelev")#printdscRas.Extentenvel=str(dscRas.Extent.XMin)+''+str(dscRas.Extent.YMin)+''+str(dscRas.Extent.XMax)+''+str(dscRas.Extent.YMax)23/33 #printenvel因为从9.3版本之后，Extent是一个Object，所以直接使用dscRas.Extent在clip语句中会出现错误，而clip语句的第二个参数Rectangle是一个Envelope类型数据，所以这里我使用了envel接收Extent的四个值。try:gp.clip_management("c:progpendatalugrid",envel,"luwood")except:printgp.GetMessages()编写一个脚本输出“c:/prog/hmbarea/tmcomp.bil”的波段个数，命名为countBands.py。在栅格数据集属性中寻找这个方法：importarcgisscripting,sysgp=arcgisscripting.create(9.3)gp.workspace="c:/prog/hmbarea"dta="tmcomp.bil"dsc=gp.Describe(dta)printdta+""+dsc.DataType+":"+str(dsc.BandCount)+"bands."运行结果：“tmcomp.bilRasterDataset:3bands.”6.2存在（Exists）关于一个数据集的一个更基本的信息是其是否存在。测试一个数据是否存在也可帮助我们避免错误。一个非常有用的技术是通过gp.Exists判断一个特定的数据是否存在。gp.Exists能用在各种类型数据上，比如在使用数据之前添加一个判断是否存在的语句：ifgp.Exists(“streams.shp”)gp.Buffer(“streams.shp”,”stbuff200.shp”,200)在下面的代码中我们将使用这种方法判断数据是否存在，其可以成为overwriteoutput设置的替代语句。判断一个字段属性是否存在。上面的方法无法测试一个数据的字段是否存在，但线面的listFields枚举是个技巧：(listFields返回一个Python属性对象列表)ifnotgp.listFields(“streams.shp”,”st-code”):我们将在后面使用这句判断字段是否存在。6.3在循环中使用描述和存在这会让你感到脚本在处理多种数据方面的强大作用。我们将裁剪一整系列的栅格数据，并存在新文件夹内。在这个处理过程中，我们也会在创建输出数据之前用“存在工具”中一个很方面的方法检查其是否已经存在。在本例中，我们使用了ArcInfoWorkspace，并把格网栅格存在这里，当然你也可以将他们存在Geodatabase中。代码如下：24/33 importarcgisscriptinggp=arcgisscripting.create(9.3)#gp.checkextension("spatial")gp.workspace="c:/prog/hmbarea"ifnotgp.Exists(gp.workspace+"/cities_polygon.shp"):gp.FeatureclassToShapefile_conversion("c:/prog/pendata/cities/polygon",gp.workspace)#因为没有找到citypoly.shp和city栅格数据，所以我使用pendata文件夹下的cities里的polygon转换的shp文件取代。else:printgp.workspace+"/cities_polygon.shpexists"ifnotgp.Exists(gp.workspace+"/cities"):gp.PolygonToRaster_conversion("cities_polygon.shp","CITY_CODE","cities","","","60")else:printgp.workspace+"/citiesexists"ifnotgp.Exists(gp.workspace+"hmbcity"):gp.createArcInfoWorkspace(gp.workspace,"hmbcity")else:printgp.workspace+"/hmbcityexists"try:gp.mask="cities"gp.checkoutextension("spatial")#检查有spatial的许可raslist=gp.listrasters()forrasinraslist:outputname=gp.workspace+"/hmbcity/"+rasprintoutputnameifnotgp.exists(outputname):try:gp.ExtractByMask_sa(ras,gp.mask,outputname)#用掩膜裁剪except:gp.getmessages()except:gp.getmessages()探索1：selSanMateoSel.py用Select_Analysis创建SanMateo（city_code=26）城市的shp文件。创建一个小的工作空间命名为“Sanmateo”，然后用SanMateo.shp裁剪pendata文件夹下文件名以‘g’开头（“g*”）的coverage文件的polygon，存入此空间。建议：你需要使用gp.listdatasets(“g*”)代替gp.listfeatureclasses，因为你需要寻找coverage，在你处理列表里成员时，可以将成员名称赋给变量f，然后用gp.exists(f+”/polygon”)查找多边形要素集。你也会发现gp.CreateFolder_management很有用，你也可以裁剪所有数据时设置通配符为“*”，或者不适用通配符。代码如下：25/33 importarcgisscriptinggp=arcgisscripting.create(9.3)gp.overwriteoutput=1gp.workspace="c:/prog/pendata"workspaceLoc="c:/prog"try:gp.Select_Analysis("cities/polygon","San_Mateo.shp",'"CITY-CODE"=26')wsName="San_Mateo"newWS=workspaceLoc+"/"+wsNameifnotgp.exists(newWS):gp.CreateFolder_management(workspaceLoc,wsName)fList=gp.listdatasets("g*")forfinfList:ifgp.exists(f+"/polygon"):gp.Clip_analysis(f+"/polygon","San_Mateo.shp",newWS+"/"+f)except:printgp.getmessages()探索2：multiFeatureRas.py创建两个Python列表：cov=[“geology”,”landuse”,”publands”,”flood”,”streams”,”roads”,”vegetation”]fld=[“TYPE-ID”,”LU-CODE”,”PUB-CODE”,”FLOOD-CODE”,”ST-CODE”,”ROAD-CODE”,”VEG-CODE”]将每个coverage根据对应的字段转换为栅格，分辨率设置为60。代码如下：importarcgisscriptinggp=arcgisscripting.create(9.3)gp.workspace="c:/prog/pendata"cov=["geology","landuse","publands","flood","streams","roads","vegetation"]fld=["TYPE-ID","LU-CODE","PUB-CODE","FLOOD-CODE","ST-CODE","ROAD-CODE","VEG-CODE"]gp.overwriteoutput=1ifnotgp.Exists("multiFeatRas"):#新建一个文件夹保存处理后数据gp.CreateFolder_management(gp.workspace,"multiFeatRas")foriinrange(len(cov)):#遍历coverage数据ifgp.Exists(cov[i]+"/polygon"):covTop=cov[i]+"/polygon"elifgp.Exists(cov[i]+"/arc"):covTop=cov[i]+"/arc"gp.FeatureToRaster_conversion(covTop,fld[i],"multiFeatRas/"+cov[i]+"g","60")26/33 7在Python脚本中使用地图代数（MapAlgebra）将AML程序移植到Python中唯一令我犹豫的是它不能像在AML中那样编写地图代数语句来处理格网数据。但是，有一种方法可以很好地解决这个问题——简单修改一个定义好的函数。SpatialAnalyst工具下的MultiOutputMapAlgebra工具允许使用完整的地图代数赋值语句，包括嵌套结构（！）。唯一的遗憾是这工具名字太长，所以我已经把它嵌入一个自定义函数中，命名为“ma”，当然函数中包括开始和结束“Spatial”扩展功能。下面的例子说明了这种用法，虽然只有一点语句，但包括了在处理多工具中的两个棘手的问题。usingMapAlgebra.pyimportarcgisscriptinggp=arcgisscripting.create(9.3)defma(expression):gp.checkoutextension("spatial")gp.MultiOutputMapAlgebra_sa(expression)gp.checkinextension("spatial")gp.workspace="c:proghmbarea"gp.overwriteoutput=1ifnotgp.exists("streamsg"):gp.FeatureToRaster_conversion("streams/arc","st-code","streamsg","60")gp.CellSize=60ma("slopeg=slope(elev,0.3048)")ma("stream2=con(isnull(streamsg),0,streamsg)")ma("histreams=(elev>1000)andstream2")另外一个方法是直接将方法赋值给变量，这样通过下面简单的一句话就可以实现前面所讲的自定义函数了：ma=gp.MultiOutputMapAlgebra_sa但也要确保在初始化ma变量之前检查“Spatial”扩展功能，安全的编程习惯是你在Except:后添加gp.CheckInExtension(“spatial”)，这是为了避免在你尝试CheckOutExtension时，其已经被检查过了，如下所示：gp.checkoutextension("spatial")ma=gp.MultiOutputMapAlgebra_sama("slopeg=slope(elev,0.3048)")ma("stream2=con(isnull(streamsg),0,streamsg)")ma("histreams=(elev>1000)andstream2")gp.checkinextension("spatial")27/33 8数据管理和指针（DataManagementandCursors）8.1数据管理（DataManagement）处理表格数据、创建字段来储存这些数据或者其他数据操作时GIS中非常重要的工作。的那个你需要做一连串数据管理和分析时，脚本是一个很好的选择。在一些情况下，我们需要处理数据字段时有一系列工具可以使用，比如，“AddFields”、“Calculatevaluesforfields”、“deletefields”以及通过一个公共字段链接表格获取附加数据字段。其中有些工具会利用不同的统计方法将输入数据输出为新的表格。其他情况下，我们需要处理行数据，比如栅格数据中单个要素或栅格数据的属性值，接下来我们学习使用指针（Cursor）逐个处理行数据。首先，我们看一下一些能够在处理表格数据（主要是包括属性字段）时用到的工具。有些工具根据属性选择记录，或复制或删除选中记录。ToolboxTool做什么？Outputfrequency频率统计Tablestatistics总结统计Table利用where子句选Analysisselect新要素集择要素用SQL语句选择并table_selectTable提取选中的属性addField添加字段表格中字段通过表达式为字段calculateField已存在字段的值赋值deleteField删除字段为点要素添加X&YaddXYX&Y值值对图层或表格通过DataManagementselectLayerBy属性查询选择、更新Alias:managementAttribute或移除选择copyFeatures复制选中的要素新要素集deleteFeatures删除选中的要素根据一个共同的字addJoin段将一表格连接到连接关系图层（或表格）removeJoin移除已存在连接根据逻辑表达式提输出Coverage（包含select(reselect)取Coverage地图要属性表）素Coverage为信息表添加属性Alias:Arcadditem字段从信息表中移除属dropitem性字段28/33 根据相关属性融合joinitemTable数据表为point、label、node字段添加到已存在addxy表格添加x、y字段表格创建新脚本（addCalcField.py），添加字段“elevm”到“contours.shp”（surf_bld文件夹内），并计算其值为“[elev]*0.3048”。代码如下：importarcgisscriptinggp=arcgisscripting.create(9.3)gp.workspace="c:/prog/surf_bld"gp.toolbox="management"#设置gp.toolbox为“management”。try:ifnotgp.listfields("contour.shp","elevm"):#判断此字段是否已经存在gp.addField("contour.shp","elevm","float")gp.calculateField("contour.shp","elevm","[elev]*0.3048")except:printgp.GetMessages(2)创建新脚本（addXY.py）使用management里的addxy工具为Marbes文件夹下的“samples.shp”添加x&y值。代码如下：importarcgisscriptinggp=arcgisscripting.create(9.3)gp.overwriteoutput=1gp.workspace="c:/prog/marbles"gp.addxy_management("samples.shp")8.2指针（Cursors）指针给了一种访问数据值的通道，允许遍历数据表格中的所有记录。由于记录（行数据）可能是一个矢量要素或一栅格值，所以，指针在处理数据是十分强大。指针有三种类型：SearchCursor：读取一行中的值InsertCursor：插入新的行UpdateCursor：改变行中值以及删除行尝试下面的代码，使用了SearchCursor来显示addCalcField.py计算的结果：importarcgisscriptinggp=arcgisscripting.create(9.3)gp.workspace="c:/prog/surf_bld"cur=gp.SearchCursor("contour.shp")row=cur.Next()whilerow:printrow.elev,row.elevm29/33 row=cur.Next()尝试一下代码，介绍了读取和显示“stream.shp”内所有顶点，不要忘记while循环内Next()语句，否则会出现无限循环！也展示了几何要素的用法！importarcgisscriptinggp=arcgisscripting.create(9.3)rows=gp.SearchCursor("c:/prog/surf_bld/stream.shp")row=rows.Next()whilerow:feat=row.shapea=0whilea