ArcGIS中Python脚本学习.docx
《ArcGIS中Python脚本学习.docx》由会员分享，可在线阅读，更多相关《ArcGIS中Python脚本学习.docx（49页珍藏版）》请在冰点文库上搜索。
ArcGIS中Python脚本学习
[初学入门]ArcGIS中Python脚本的使用
By：
飞天小猪
写在前面的话
一直想学习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脚本中，这些脚本可以生成脚本工具，像其他地理处理工具一样使用。
一、Python语言基础
安装PythonWin，在…\ArcGisDesktop9.3.iso\Desktop\PythonWin目录下可以找到PythonWin的安装程序，默认是不安装的，
。
同时会安装win32com以及允许任何脚本在基于Dispatch的地理处理过程中工作。
ArcGIS10中引入了全新的PythonWindow来增强内嵌的Python体验。
警告：
不要尝试更新随ArcGIS安装的Python到一个新的版本！
下面介绍Python的一些简单语法和规则。
1数学运算符
Python提供了多样化的通用数学运算符——多数编程语言的特征，以及许多通过import的modules提供的符号。
常用的有+,-,*,/,**（幂）,%（取模，即除后的余数）。
下面的表格显示了整型（Integer）和浮点型（Float）各种组合运算的结果，记住一条规则，只要参与运算的有浮点型，则结果为浮点型；全为整型时，结果才为整型。
输入表达式
结果
Notes
2+3
5
整型结果
2.+3
5.0
2.是浮点型，结果浮点型
2-3
-1
2*3
6
整型结果
2.*3
6.0
浮点型
5/2
2
整型
5./2
2.5
5%2
1
取模
Az=270
Newaz=az+180
Printnewaz%360
90
取模的用途之一——方位角加180后逆转方向
5**2
25
25**0.5
5.0
没有sqrt（）功能，除非添加math模块
2字符串操作
注：
使用Python帮助：
有超过30种内置方法来处理字符，请到SequenceTypes下的StringMethods寻找帮助！
字符串是一串字母，比如’SanFrancisco’，字符串下标从0开始。
学习字符串语法的最好方法是自己动手尝试，下标展示之：
输入
结果
Notes
print'zhulj'.capitalize（）
Zhulj
s.capitalize（）即将capitalize（）方法用于s
s='zhulj'
prints.capitalize（）
prints[0]
z
Strings可以像一个字母列表一样处理，第一个字母下标为0，某个字符段可以用1:
3来格式化：
从第1个的开头到第3个的开头，不包括下标为3的字母；s[-1]表示倒过来第一个，相当于s[len（s）-1]
s1=s[1]
prints1
h
prints[-2:
]
lj
prints[2:
3]
u
prints[2:
4]
ul
prints[2:
],s[:
5]
uljzhulj
s2=s.upper（）
prints2
ZHULJ
我们可以将字符串方法的结果赋给新的变量
s3=s+s2
prints3
zhuljZHULJ
字符串组合用“+”
prints*3
zhuljzhuljzhulj
字符串重复用“*”，后为重复次数
selstr='"elev">1000'
printselstr
"elev">1000
字符串可以使用单引号或双引号，跨行时用双引号。
othersel=”’elev’>1000”
printothersel
‘elev’>1000
prints.isupper（）
False
一些方法返回值为布尔型（True或False），一些返回索引值（下标值）
prints2.isupper（）
True
p='d:
/work/lu.shp'
printp.find（‘.’）
10
printp.find（‘/’）
2
plist=p.split（'/'）
printplist
['d:
','work','lu.shp']
你可以用split（）方法解析出不同的字符串片段，并创建一个列表（List），我们可以使用其中不同的元素
printplist[0]
d:
printplist[1]
work
p2='d:
\\work\\soil.shp'
printp2
d:
\work\soil.shp
反斜线“\”和某些字母一起有特殊用法，如\n为换行，“\”为转义字符，如“\\”则表示“\”
print'Jerry\'sKids'
Jerry'sKids
print'Jerry\'s\nKids'
Jerry's
Kids
p3=r'd:
\work\soil.shp'
printp3
d:
\work\soil.shp
字符串前加“r”则强制“\”代表其本身，而非转义字符，这对于文件路径的操作很方便
3模块的使用（Modules）
Python提供了一系列内置的方法（大量依赖于模块）用于通用编程。
Python安装时自带了大量Modules，最常用的有math,sys,random,array以及os.path。
当然还有好多Modules可以下载，比如数字处理（Numeric）——numpy，可在www.python.org或里搜索。
www.python.org/moin/NumericAndScientific页面中列举了一些。
使用Module前，必须import之。
通常我们会将一行import放在程序顶部，比如：
importarcgisscripting
当然，这不必成为你程序的第一行，但必须在使用它里面方法之前。
当要引用多个模块是，中间用逗号分隔，比如：
importarcgisscripting,sys,string,os,math
我们也可以自己为频繁使用的方法创建Module，下面，我们开始体验内置的Modules。
math和random模块
很多常用的数学计算功能都可以通过math找到，比如三角计算或对数计算，如果要使用复杂数字，就使用cmath模块。
和之前一样，通过以下表格来体现模块的使用：
输入
结果
Notes
importmath
printmath.log10（100）
2.0
以10为底的对数
printmath.log（100）
4.599
自然对数
printmath.pi
3.
π是一个静态常量，所以不需要括弧
pi=math.pi
printpi
3.
如果不想总是输入“math.pi”可以将其赋给一变量
pi
3.
不需要print即可查看变量值
printmath.sin（radians）
printmath.cos（radians）
printmath.tan（radians）
三角函数的计算是弧度制
degrad=pi/180
45*degrad
0.785398
度转化为弧度
sin=math.sin
sin（45*degrad）
sin（90*degrad）
0.6
1.0
即使功能函数（像sin）都可以赋给一个变量
math.e
2.71828
math.hypot（3,4）
5.0
此方法是求三角形的斜边
x1=520382;y1=4152373
x2=520475;y2=4152963
不同赋值语句间用“；”分隔
xr=x2-x1
yr=y2-y1
math.hypot（xr,yr）
math.sqrt（xr**2+yr**2）
（xr*xr+yr*yr）**0.5
597.956189
不同的方式，相同的结果
importrandom
random.random（）
0.8
random（）方法，每次结果都不同，值域为[0.0,1.0）
rnd=random.random
rnd（）
0.
mu=50
s=10
printrandom.gauss（mu,s）
46.4
4使用def构建函数
有点像Module，但更简单，函数是一个自己定义功能，用在之后的代码中，并且提供任何你想要使用的参数。
这个函数从此可像变量那样在程序中使用，结合例子更容易理解。
接下来的代码定义了一个将度转换为弧度的简单函数，同时也定义了一个弧度转换为度的函数，它们和Excel内置的函数类似。
importmath
defradians（angdeg）:
returnangdeg*math.pi/180
defdegrees（angrad）:
returnangrad*180/math.pi
printmath.sin（radians（45））
printdegrees（math.acos（0.5））
运行之，得到结果：
0.118760.0
5流程控制结构：
If，While，For
任何脚本或编程语言的一个重要特征就是执行一系列不同情形语句的能力。
你想要创建一系列山影栅格来代表夏天、冬天和春秋分。
山影（hillshade）工具需要有太阳高度角和方位角作为输入参数。
重要日期
太阳倾角
夏至（6月21日）
23.44
春秋分（3月21日，9月21日）
0
冬至（12月21日）
-23.44
接下来是一段相当简单的代码，通过太阳倾角（太阳光线正午垂直照射的纬度）获取太阳角和方位角以及纬度。
输入两个参数：
lat（研究区域的纬度，南半球为负）和decl（太阳倾角），由此得到sunangle和azimuth：
lat=30
decl=20
sunangle=90-lat+decl
azimuth=180
ifsunangle>90:
sunangle=180-sunangle
azimuth=0
printsunangle,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.path
if"d:
/testfolder/test.txt"）:
print"测试文件夹存在"
print"txt文件存在"
elif"d:
/testfolder"）:
print"测试文件夹存在"
print"测试文件夹存在，但txt文件不存在"
else:
print"两者都不存在"
可选探索示例
接下来的例子做的事情对GIS非常重要，但是实际上不用任何地理处理代码。
USGS7.5米分辨率DEM（数字高程模型）是文本文件（USGSDEM文件），投影为UTM，UTM北向和东向单位是米，但是高程单位可能是英尺（feet）或米（meters）。
因此在获取垂直或水平距离信息时会有问题，比如坡度可以通过垂直距离/水平距离获得。
如果你不在使用Z值之前设置为0.3048，将会出现错误结果。
但是不幸的是，你可能不知道DEM文本文件的垂直单位是英尺还是米。
这些信息保存在第539个字符里，“1”代表英尺，“2”代表米，所以可以通过读取这个文件判断。
下面的脚本演示了上述内容：
importfileinput
infile=r"c:
\prog\pendata\woodside.dem"
firstline=fileinput.input（infile）[0]
unitchar=firstline[539]
unit="（unknown:
nota7.5'DEM?
）"
ifunitchar=="1":
unit="feet"
ifunitchar=="2":
unit="meters"
print"\nElevationin"+""+unit
fileinput.close（）
输出结果：
Elevationinfeet
while（continued）
Ø运行下面的代码，说明了一种while循环：
x=1
whilex<10:
printx
x=x+1
屏幕依次输出1~9
Ø下面说明一下“==”（等于）的概念：
x=5
z=x==4
printz
输出False
x=5
z=x==5
printz
输出True
“==”是逻辑运算符之一，其他有“<”（小于）、“>”（大于）、“>=”（大于或等于）、“<=”（小于或等于）、“<>”（不等于）。
使用逻辑运算符计算得到结果为布尔型：
true
（1）和false（0）。
下面例子简单体会一下布尔型表达式：
x=1
whilex<10:
printx
x=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’结尾）
importos
ws="c:
/prog/hmbarea"
ilist=os.listdir（ws）#创建一个列表保存工作文件夹内的文件
fcs=[]#创建一个空列表，保存结尾为‘.shp’的文件
foriinilist:
ifi.endswith（".shp"）:
fcs.append（i）
forfcinfcs:
printfc（输出结果如右图所示）
Ø下面这个例子的循环较多次数。
变量mu是算术平均数，s是标准差——这两个是random.gauss（）用到的参数，可以尝试改变n的值查看结果！
importrandom
mu=50
s=10
z5=mu-s*1.96
z95=mu+s*1.96
count=0
n=100000
foriinrange（n）:
x=random.gauss（mu,s）
ifxz95:
count=count+1
printfloat（count）/n（每次运行的结果都不同，但都在0.05左右，即统计结果在5%左右）
6简单输入和输出
我们现在利用前面计算太阳角代码的片段，之前是直接指定参数值，现在我们有很多种方法提供输入参数，现在我们用sys方法。
尝试下面的代码，点击
（run运行）之后，在对话框中函数自变量（Arguments）一栏填入：
4023.44，如下图：
importsys
lat=float（sys.argv[1]）
decl=float（sys.argv[2]）#使用arguments（argv）方法从“Arguments”一栏中获取输入参数，并指定一个浮点型转换将字符型输入值传递给lat和decl
sunangle=90-lat+decl
azimuth=180
ifsunangle>90:
sunangle=180-sunangle
azimuth=0
print"正午太阳角="+str（sunangle）
print"方位角="+str（azimuth）（结果：
正午太阳角=73.44方位角=180）
二、ArcGIS&Python
1如何创建地理处理对象（geoprocessorobject）
所有geoprocessing的Python脚本都可以通过importarcgisscripting或者win32com去穿件geoprocessorobject。
下面的例子显示二者区别：
arcgisscriptingmodule需要在；通过win32com创建的geoprocessor可以在不同的Python版本上运行。
#9.3
importarcgisscripting
gp=arcgisscripting.create（9.3）
gp.workspace=”c:
/Tongass”
gp.clip_analysis（“standb4”,”clipcov”,”standb4_clip”,”POLY”.”1.25”）
#Dispatch
importwin32com.client
gp=
gp.workspace=”c:
/Tongass”
gp.clip_analysis（“standb4”,”clipcov”,”standb4_clip”,”POLY”.”1.25”）
✧理解ArcGIS中数据多样性格式对我们理解地理处理工具很有帮助。
比如，特征数据可能是①单个shp文件；②geodatabase（地理数据库，我们可能指定地理数据库为工作空间）；③多边形、弧或点要素的coverage数据。
当我们想遍历工作空间里的coverage文件时，应使用ListDatasets而不是ListFeatureClasses。
2获取地理处理帮助
如果你基本熟悉了Python的语法，便可以开始熟悉ArcGIS的Geoprocessing啦，获取这些方法帮助的途径有两个：
1ArcGIS帮助系统，GoToGeoprocessing/Automatingyourworkwithscripts/Scriptingobject:
ProperitesandMethods.
这里你会看到GeoprocessorObject,这个能让你很方便地获得所有对你有用的条目、设置和其他操作文档。
比如，你想得到特定类型的文件，就找到listfeatureclasses方法：
fcList=gp.ListFeatureClasses（“w*”,”polygon”）
注：
此方法的语法为：
object.ListFeatureClasses（{wildCard}AsString,{FeatureType}AsString）AsPythonList
{}=optionalwildcard为通配符，和星号（*）组合使用，如果没有通配符则返回工作目录下的所有featureclasses。
2GeoprocessorProgrammingModel（PDF），包含方法（左箭头表示）、属性（可读写的表示为杠铃形，只读的表示为部分杠铃形）
2.1举例：
如何使用GeoprocessorProgrammingModel中的Lists
Lists（列表）及其属性和方法在图表中用紫色标出，如下：
现在我们试着编写一段脚本列举出属性表中所有属性值（Fields）（以hmbarea栅格土地利用为例，文件存在c:
\prog\hmbarea下）
importarcgisscripting,sys
gp=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（此时输出结果如右图）
3使用地理处理工具——Toolboxes和Aliases
总所周知，地理处理工具在脚本中的使用和ArcToolbox中相同，但是需要提供工具名称来
升级会员 