ARCGIS与python（3/3）

1185234520

7 在 在 Python  脚本中使用地图代数（Map Algebra ）
将 AML 程序移植到 Python中唯一令我犹豫的是它不能像在AML 中那样编写地图代
数语句来处理格网数据。但是，有一种方法可以很好地解决这个问题——简单修改一个
定义好的函数。Spatial Analyst 工具下的 MultiOutputMapAlgebra 工具允许使用完整的地
图代数赋值语句，包括嵌套结构（！ ） 。唯一的遗憾是这工具名字太长，所以我已经把它
嵌入一个自定义函数中，命名为“ma” ，当然函数中包括开始和结束“Spatial”扩展功
能。下面的例子说明了这种用法，虽然只有一点语句，但包括了在处理多工具中的两个
棘手的问题。
usingMapAlgebra.py
import arcgisscripting
gp=arcgisscripting.create(9.3)
def ma(expression):
gp.checkoutextension("spatial")
gp.MultiOutputMapAlgebra_sa(expression)
gp.checkinextension("spatial")
gp.workspace="c:\prog\hmbarea"
gp.overwriteoutput=1
if not gp.exists("streamsg"):
gp.FeatureToRaster_conversion("streams/arc","st-code","streamsg","60")
gp.CellSize=60
ma("slopeg=slope(elev,0.3048)")
ma("stream2=con(isnull(streamsg),0,streamsg)")
ma("histreams=( elev > 1000) and stream2 ")
另外一个方法是直接将方法赋值给变量，这样通过下面简单的一句话就可以实现前
面所讲的自定义函数了：
ma=gp.MultiOutputMapAlgebra_sa
但也要确保在初始化 ma 变量之前检查“Spatial”扩展功能，安全的编程习惯是你
在 Except: 后 添 加 gp.CheckInExtension(“spatial”) ， 这 是 为 了 避 免 在 你 尝 试
CheckOutExtension 时，其已经被检查过了，如下所示：
gp.checkoutextension("spatial")
ma=gp.MultiOutputMapAlgebra_sa
ma("slopeg=slope(elev,0.3048)")
ma("stream2=con(isnull(streamsg),0,streamsg)")
ma("histreams=( elev > 1000) and stream2 ")
gp.checkinextension("spatial")

8  数据管理和指针（Data Management and Cursors ）
8.1  数据管理（Data Management ）
处理表格数据、创建字段来储存这些数据或者其他数据操作时 GIS 中非常重要的工
作。 的那个你需要做一连串数据管理和分析时， 脚本是一个很好的选择。 在一些情况下，
我们需要处理数据字段时有一系列工具可以使用，比如， “Add Fields” 、 “Calculate values
for fields” 、 “delete fields”以及通过一个公共字段链接表格获取附加数据字段。其中有
些工具会利用不同的统计方法将输入数据输出为新的表格。其他情况下，我们需要处理
行数据，比如栅格数据中单个要素或栅格数据的属性值，接下来我们学习使用指针
（Cursor）逐个处理行数据。
首先， 我们看一下一些能够在处理表格数据 （主要是包括属性字段） 时用到的工具。
有些工具根据属性选择记录，或复制或删除选中记录。
Toolbox  Tool  做什么？  Output
Analysis
frequency  频率统计  Table
statistics  总结统计  Table
select
利用 where 子句选
择要素
新要素集
table_select
用 SQL 语句选择并
提取选中的属性
Table
Data Management
Alias:management
addField  添加字段  表格中字段
calculateField
通过表达式为字段
赋值
已存在字段的值
deleteField  删除字段
addXY
为点要素添加 X&Y
值
X&Y 值
selectLayer By
Attribute
对图层或表格通过
属性查询选择、 更新
或移除选择
copyFeatures  复制选中的要素  新要素集
deleteFeatures  删除选中的要素
addJoin
根据一个共同的字
段将一表格连接到
图层（或表格）
连接关系
removeJoin  移除已存在连接
Coverage
Alias:Arc
select(reselect)
根据逻辑表达式提
取 Coverage 地图要
素
输出 Coverage （包含
属性表）
additem
为信息表添加属性
字段
dropitem
从信息表中移除属
性字段

joinitem
根据相关属性融合
数据表
Table
addxy
为 point、 label、 node
表格添加 x、y 字段
字段添加到已存在
表格
  创建新脚本（addCalcField.py） ，添加字段“elevm”到“contours.shp” （surf_bld 文
件夹内） ，并计算其值为“[elev]*0.3048” 。
代码如下：
import arcgisscripting
gp=arcgisscripting.create(9.3)
gp.workspace="c:/prog/surf_bld"
gp.toolbox="management"#设置 gp.toolbox 为“management” 。
try:
if not gp.listfields("contour.shp","elevm"):#判断此字段是否已经存在
gp.addField("contour.shp","elevm","float")
gp.calculateField("contour.shp","elevm","[elev]*0.30
48")
except:
print gp.GetMessages(2)
  创建新脚本（addXY.py）使用 management 里的 addxy 工具为 Marbes 文件夹下的
“samples.shp”添加 x&y 值。
代码如下：
import arcgisscripting
gp=arcgisscripting.create(9.3)
gp.overwriteoutput=1
gp.workspace="c:/prog/marbles"
gp.addxy_management("samples.shp")
8.2  指针（Cursors ）
指针给了一种访问数据值的通道，允许遍历数据表格中的所有记录。由于记录（行
数据）可能是一个矢量要素或一栅格值，所以，指针在处理数据是十分强大。
指针有三种类型：
  SearchCursor：读取一行中的值
  InsertCursor：插入新的行
  UpdateCursor：改变行中值以及删除行
  尝试下面的代码，使用了 SearchCursor 来显示 addCalcField.py 计算的结果：
import arcgisscripting
gp=arcgisscripting.create(9.3)
gp.workspace="c:/prog/surf_bld"
cur=gp.SearchCursor("contour.shp")
row=cur.Next()
while row:
print row.elev,row.elevm

row=cur.Next()
  尝试一下代码，介绍了读取和显示“stream.shp”内所有顶点，不要忘记 while 循
环内 Next()语句，否则会出现无限循环！也展示了几何要素的用法！
import arcgisscripting
gp=arcgisscripting.create(9.3)
rows=gp.SearchCursor("c:/prog/surf_bld/stream.shp")
row=rows.Next()
while row:
feat=row.shape
a=0
while a<feat.PartCount:
stArray=feat.GetPart(a)
stArray.Reset
pnt=stArray.Next()
while pnt:
print str(pnt.id)+";"+str(pnt.x)+";"+str(pnt.y)
pnt=stArray.Next()
a=a+1
row=rows.Next()
&#61656;  尝试下面代码， 展示了 InsertCursor 的用法， 同样的基本方法可用来读取外部数据，
比如文本文件。
import win32com.client, sys, math, string
gp = win32com.client.Dispatch("esriGeoprocessing.GPDispatch.1")
#import arcgisscripting
#gp=arcgisscripting.create(9.3) 也可以这样引用
gp.workspace = "c:/prog/Marbles"
try:
cur = gp.InsertCursor("mvalley_pts.shp")
feat = cur.NewRow()
feat.id = 12
feat.Name = "Sky High Lake Camp"
pnt = gp.CreateObject("Point")
pnt.x = 485339
pnt.y = 4600001
feat.shape = pnt
cur.InsertRow(feat)
del cur
except:
print gp.getmessages()
if cur: del cur
挑战 1 ：
下面的代码可以“打开” 、 “写入” 、 “关闭”一个文本文件，修改上面的脚本输出
节点值信息。(对“stream.shp”的操作)

import arcgisscripting,os,sys
gp=arcgisscripting.create(9.3)
wspath="c:/prog/surf_bld"
txtfile="strout.txt"
fpath=wspath+"/"+txtfile
if gp.exists(fpath):os.remove(fpath)
txtfile=open(fpath,"w")
gp.workspace=wspath
rows=gp.SearchCursor("stream.shp")
row=rows.Next()
while row:
feat=row.shape
a=0
while a<feat.partCount:
geomArray=feat.getpart(a)
fid=row.GetValue("ID")
geomArray.Reset
pnt=geomArray.Next()
while pnt:
txtfile.write(str(fid) + "; " + str(pnt.x) + "; " + str(pnt.y)+"\n")
pnt=geomArray.Next()
a=a+1
row=rows.next()
txtfile.close()
挑战 2 ：
写一个可以读取刚才输出的文本文件并且创建一个新的 shp 文件的脚本。
如果你使用的是 ArcGIS10， 这里有最新的帮助， 如果使用的是 ArcGIS9.3.1 或 9.3，
这里有相应帮助。
import arcgisscripting,sys,fileinput
gp=arcgisscripting.create(9.3)
inFile="c:/prog/surf_bld/strout.txt"
gp.workspace="c:/prog/surf_bld"
fcName="newsamp.shp"
gp.overwriteoutput=1
gp.CreateFeatureClass_management(gp.workspace,fcName,"point")
try:
cur=gp.InsertCursor(fcName)
for line in fileinput.input(inFile):
values=line.split(";")
id=int(values[0])
x=float(values[1])
y=float(values[2])
feat=cur.NewRow()

pnt=gp.CreateObject("Point")
pnt.x=x
pnt.y=y
feat.shape=pnt
feat.id=id
cur.InsertRow(feat)
print id,x,y
del cur
fileinput.close()
except:
print gp.getmessages()
if cur:
del cur
接下来该做什么呢？
我们已经做了很多， 但是大多数都是创建处理特定事情的部分， 比如遍历一系列
数据。但是到目前为止我们并没有花太多时间在这些循环上。
另一间我们没有花太多时间处理的事情是， 让脚本工具更好地工作。 输入和输出
是这个处理过程中有挑战性的部分， 所以在这方面下手可能会更适合你。 我已经发
现这是最令人困惑的部分， 所以我将能发现的尽可能多的输入/输出情况做了总结，
请看附录 1。
最后一个建议是在模型中使用你的脚本工具，这是一个非常值得和有用的工作。
附录 1 ：地理处理脚本中输入& 输出方法指南
1、 输入数据集（Datasets）
&#61548;  脚本：使用 GetParameterAsText 或 sys.argv
&#61548;  工具参数： （要素类、栅格数据或图层等）必需的输入
2、 输出数据集（Datasets） （如，Spatial Statistics/测量地理分布里的平均中心）
&#61548;  脚本：使用 GetParameterAsText 或 sys.argv
&#61548;  工具参数： （要素类、栅格数据等）必需的输入
3、 输出值和字段（如，Spatial Statistics/Utilities 下的计算面积）
&#61548;  脚本：输入和输出的数据集、复制输入到输出
&#61548;  工具参数：输入要素类、输出要素类
4、 输出单个值（如，Spatial Statistics/Analyzing Patterns 下的空间自相关）
&#61548;  脚本：SetParameterAsText
&#61548;  工具参数： （Long、Double 等）获取的输出
5、 输入和输出表格或要素类 （如， Management/Fields 下的 Transpose Time fields）
&#61548;  脚本：GetParameterAsText
&#61548;  工具参数：输入—表格视图，必需的输入；输出—表格，必需的输出
6、 修改输入到输出， 没有找到简单的例子， 除了Mosaic （Samples/Conversion/Raster）
&#61548;  脚本：输出用 GetParameterAsText ；类似的 SetParameterAsText

&#61548;  工具参数：输出是必须的输入，然后作为获取的输出，设置从原始处获取
7、 批处理（如，Management/Projections 和 Transformations/Feature 下的 Batch
Project）
&#61548;  脚本：输入和输出用 GetParameterAsText，获取的输出不需要
&#61548;  工具参数：输入地理数据集，必需的输入，允许多值；输出工作空间，工
作空间或要素数据集，必需的输入；输出空间的单独的获取的输出也是一
样
8、 如果成功完成则返回True （如， Management/Raster下的Batch Calculate Statistics）
&#61548;  脚本：nothing
&#61548;  工具参数：布尔型参数，获取的输出默认为 True
什么时候用 Derived（获取的，衍生的）？如果工具更新输入到工具里，或者
如果你需要一个输出并用在一个模型里，但是你不想让输出参数在脚本对话
框中出现。
附录 2 ：其他
1、 工作空间和临时工作空间设置
支持“临时工作空间”环境设置的工具可将指定的位置用作输出数据集的默认
工作空间。 “临时工作空间”专门用于存放不愿保留的输出数据。
“临时工作空间”环境的主要用途是供“模型构建器”使用。 “模型构建器”
需要使用一个工作空间来写入中间数据集（模型运行后便不再使用） 。尽管它
主要服务于“模型构建器” ，但有时可能也需要为各工具对话框设置临时工作
空间。
使用工具对话框时， 输出数据集名称会按照当前工作空间环境和临时工作空间
环境设置自动生成。
2、 Python 和 PythonWin 的版本
不要急于升级你的 Python 或 PythonWin 版本。最好是使用随 ArcGIS 安装的版
本，而且 PythonWin 不是自动安装的，在前文中有介绍。

panda小正太

学习学习！

wy1018651314