ARCGIS与python（2/3）

1185234520

二、ArcGIS&Python
1  如何创建地理处理对象（geoprocessor object ）
所有 geoprocessing 的 Python 脚本都可以通过 import arcgisscripting 或者 win32com 去穿
件 geoprocessor object。下面的例子显示二者区别：arcgisscripting module 需要在
Python2.5.1 版本上创建并且需要此版本创建 geoprocessor；通过 win32com 创建的
geoprocessor 可以在不同的 Python 版本上运行。
#9.3
import arcgisscripting
gp=arcgisscripting.create(9.3)
gp.workspace=”c:/Tongass”
gp.clip_analysis(“standb4”,”clipcov”,”standb4_clip”,”POLY”.”1.25”)
#Dispatch
import win32com.client
gp=win32com.client.Dispatch(“esriGeoprocessing.GpDispatch.1”)
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 啦，
获取这些方法帮助的途径有两个：
① ArcGIS 帮助系统，Go To Geoprocessing/Automating your work with scripts/Scripting
object:Properites and Methods.
这里你会看到 Geoprocessor Object,这个能让你很方便地获得所有对你有用的条目、
设置和其他操作文档。比如，你想得到特定类型的文件，就找到 listfeatureclasses 方
法：fcList=gp.ListFeatureClasses(“w*”,”polygon”)
注： 此方法的语法为： object.ListFeatureClasses({wildCard} As String, {FeatureType}
As String) As Python List
{} = optional wildcard 为通配符，和星号（*）组合使用，如果没有通配符则返回
工作目录下的所有 feature classes。
② Geoprocessor Programming Model（PDF） ，包含方法（左箭头表示） 、属性（可读写
的表示为杠铃形，只读的表示为部分杠铃形）
2.1  举例：如何使用 Geoprocessor Programming Model  中的 Lists
Lists（列表）及其属性和方法在图表中用紫色标出，如下：
现在我们试着编写一段脚本列举出属性表
中所有属性值（Fields） （以 hmbarea 栅格土地利
用为例，文件存在 c:\prog\hmbarea 下）
import arcgisscripting, sys
gp = arcgisscripting.create(9.3)
gp.workspace = " c:/prog/hmbarea"
fieldList = gp.ListFields("landuse", "*",
"all")
dsc=gp.describe("landuse")
print "landuse"+" "+dsc.DataType+":"
for field in fieldList:
... print field.Name, field.Type (此时输出结果如右图)

3  使用地理处理工具 ——Toolboxes 和 和 Aliases
总所周知，地理处理工具在脚本中的使用和 ArcToolbox 中相同，但是需要提供工具
名称来使用它们。但是记住一个名称可能有好几个工具，比如，裁剪工具（Clip）在
Coverage->Analysis->Extract 工具集里， 另一个是在 Data Management Tools 下的 Raster 工
具集下。区别是每个工具适用不同的数据类型，但是我们如何在脚本中区分它们呢？
在 ArcToolbox 中，我们可以随意选择要使用的工具，但在脚本里就必须在某些方面
加以区分。因此我们使用 Aliases（别名）——已经成为工具名称的一部分。
每一个工具都有自己的别名，我们可以通过右键->属性来查看：
Aliase  Toolbox
“conversion”  Conversion
“3d”  3D Analyst  “geocoding”  Geocoding
“analysis”  Analysis  “ga”  Geostatistical Analyst
“arc”  Coverage  “lr”  Linear Referencing
“management” Data Management  “sa”  Spatial Analyst
“cartography”  Cartography  “stats”  Spatial Statistics
现在我们用一段脚本来解释：
import arcgisscripting,sys
gp=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")
注：上面脚本用“#”注释的中文内容不要出现在脚本文件中，否则会出现错误
结果截图：

如果你一次使用一个工具集中的若干工具，可以通过环境设置省下一些文字：
gp.Toolbox = "Analysis"
gp.Buffer("streams.shp", "stbuff200.shp", 200)
gp.Clip("geol.shp", "stbuff200.shp", "geolstr200.shp")
4  在建模中使用脚本（Scripts in ModelBuilder ）
首先，需要记住的很重要的一点是，ArcToolbox 里相当数量的工具实际上都是
脚本。 脚本都有一个 图标。 比如， 空间统计分析工具 （Spatial Staistics tools）
几乎都是 Python 脚本（一些是模型中使用了脚本
） 。实际上，我们可以复制并编辑这些脚本作为其他用
途。
为了在 ModelBuilder 中使用脚本或将脚本当做 ArcToolbox 中工具使用，我们需
要考虑如何给它输入值以及让其设置输出值。仍然以太阳角计算代码为例，我们给
其加上两句引用，四句输入输出语句，就可以用作 Modelbuilder 中的一个步骤了。
编辑下面的脚本，不过不要再 PythonWin 中运行之，因为 getparameterastext
和 setparameterastext 只能用在脚本工具（ArcToolbox 或 Modelbuilder）中。
import arcgisscripting
gp=arcgisscripting.create(9.3)

lat=float(gp.getparameterastext(0))
decl=float(gp.getparameterastext(1))
sunangle=90-lat+decl
azimuth=180
if sunangle>90:
sunangle=180-sunangle
azimuth=0
gp.setparameterastext(2,str(sunangle))
gp.setparameterastext(3,str(azimuth))
这段代码中的“新面孔” ：
getparameterastext ：是 Python 传递给 geoprocessor（我们称之为 gp）的一个
方法，允许工具提供输入参数。每次你运行这个工具时，都会看到一个对话框，提
示输入参数，这个方法允许你在接下来的程序中使用。索引 0 和 1 指第一个和第二
个参数。
setparameterastext ： 和 getparameterastext 相反，它传递第二个条目（比如
str(sunangle)的值）给指定的输出参数。前两个参数索引为 0 和 1，所以接下来输出
参数的索引为 2 和 3。
然后，我们需要将脚本加进工具（Making a script into a tool） ，那样才能在
ArcToolbox 或 ModelBuilder 或 Command line 中使用。如前面提到的那样，这个脚本
只能用于工具，包括输入/输出方法是 PythonWin 不能处理的，但这些是多数工具必
需的。
  在 ArcCatalog 里，定位到 Python 脚本保存的文件夹，最好和数据在一个盘里，右
键->新建 toolbox。当然也可以使用之前创建的 toolbox。
  在 ArcToolbox 里，右键 toolbox，选添加->scripts，填写如下图：

注意：脚本文件是一个脚本工具的参考！非常重要的一点，你使用脚本创建了
一个脚本工具，但是脚本本身并没有和工具一起保存，脚本工具作为 toolbox
的一部分保存在“*.tbx”文件中。你也可以右键脚本工具点击“编辑” ，出现
PythonWin 或其他任何 IDE 窗口，这看起来好像是脚本存在于工具中。所以，
记得移动时将脚本工具文件和脚本本身一起拷贝。
“下一步”后是参数配置页面，如下图设置各参数如表格所示：
Display Name  Data type  Direction  Default
Latitude  Double  input  0
Declination  Double  input  0
Sun Angle  Double  output  35
Azimuth  Double  output  300

  现在可以运行此 toolbox 了， 对话框提示你输入 latitude 和 declination。 工具是否正
确运行呢？如果是的话，你应该可以看到“Completed” ，你可能会看到有一黑色窗
体一闪而过，不用担心。那么，它是干什么的呢？
还记得结果是输出两个数字参数，那么，这些数字哪去了呢？
很好的问题，这仅能说明你能创建一个工具，但是不能想 ArcToolbox 那样运行。比
如输出一种数据，栅格或特征数据（.shp）之类的。
  但是它能在作为 Modelbuilder 工具正常运行， 下面我们将使用它的输出参数创建一
个 hillshade。
  在你的 toolbox 中新建一个 model，将刚才创建的脚本工具（script tool）拖进来。
  双击工具，输入参数，初始化之。打开“Sun Angle”和“Azimuth”发现它们还是
默认值，说明此脚本工具还没有运行。右键单击工具，选择 Run，然后发现两个输
出参数已经改变！

需要注意的是：latitude 范围是-90~90，declination 范围是-23.44~23.44。尝试输入
latitude -70，declination 23.5，你会得到什么？为什么？
  确保你已经得到值域内的太阳角和方位角，它们将是构建 hillshade 的输入参数。
首先， 添加 hillshade 工具， 双击指定一个 elevation 栅格数据 （这里我选择了 marbles
文件夹下的elevation） ， 用下拉条指定azimuth和altitude值为azimuth和sun angle。
运行，然后右键单击输出文件，选“Add to Display”在 ArcMap 里查看结果。
  探索 1 ：
如果你想通过日期获取太阳倾角，该怎么做呢？尝试下面代码，保存为
getnoonsunfromdata.py，现在有五个参数，如下表进行正确设置：
Display Name  Data type  Type  Direction  Default
Month  Long  Required  input  1
Date  Long  Required  input  1
Latitude  Double  Required  input  0
Sun Angle  Double  Derived  output  35
Azimuth  Double  Derived  output  300
import win32com.client,sys,math
gp=win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")
#构建两个函数，首先判断输入月/日的合法性，然后获取是一年中的第几天
def jdate(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])):
return jd[im - 1] + id
else:
print "Date not in month"
return 0
else:
print "Month should be between 1 and 12 inclusive"
return 0
#计算倾角的函数
def declin(day):

rad = math.pi / 180
xlong = 279.164 + 0.985647 * day
xanom = 356.381 + 0.985600 * day
g = xanom * rad
elong = xlong + 1.915 * math.sin(g) + 0.02 * math.sin(2 * g)
elong = elong * rad
oblecl = 23.44 * rad
return math.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 + decl
azimuth = 180
if sunangle > 90:
sunangle = 180-sunangle
azimuth = 0
gp.setparameterastext(3,str(sunangle))
gp.setparameterastext(4,str(azimuth))
我们在思考远一点，如何根据一天内的任何时间而不是正午时间，获取太阳角和方
位角呢？
感兴趣的可以从本文开头的链接中下砸相应代码学习。
&#61618;  探索 2 ：
如何在 PythonWin 里运行这个脚本？首先我们得明确几个点：①我们将把 hillshade
作为脚本的一部分使用，并为其提供输入参数：一个高程栅格（elevation raster） ；
②GetParameterastext 仅在用作工具时起作用。
下表为脚本工具的参数设置：
Display name  Data type  Type  Direction  Default
Month  Long  Required  Input
Date  Long  Required  Input
Latitude  Double  Required  Input
Workspace  Workspace  Required  Input
Elevation  Raster Dataset  Required  Input
Hillshade  String  Required  Input
代码如下：
import win32com.client, sys, math
gp = win32com.client.Dispatch("esriGeoprocessing.GPDispatch.1")
#函数定义，同上，下面仅给出函数名称：
def jdate(im, id):
def declin(day):
# 主程序，使用 sys.argv[]代替 getparameterastext()
month = int(sys.argv[1])

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 + decl
azimuth = 180
if sunangle > 90:
sunangle = 180-sunangle
azimuth = 0
gp.addmessage("about to try")
try:
gp.OverwriteOutput = 1
gp.addmessage("after overwriteoutputs setting, " + gp.workspace + "/" +
hillsh)
gp.CheckOutExtension("spatial")
gp.addmessage(gp.workspace + "/" + hillsh)
gp.HillShade_sa (elev, gp.workspace + "/" + str(hillsh), azimuth, sunangle)
gp.addmessage("done with hillshade")
gp.CheckInExtension("spatial")
except:
## print gp.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 中得到一个丰富的
消息，而不是“未知错误” 。

5.1  调试选择和消息
Python 优于 AML 的优点之一是它有更好的调试方法，调试程序有很多选择，这里
不能一一列举。
&#61656;  打印语句（Print statements ）
一开始就养成良好的调试方法是： 将变量的当前值或脚本的处理过程打印在屏幕上。
比如，对之前的一段脚本加以修改：
import arcgisscripting,sys
gp=arcgisscripting.create(9.3)
gp.Workspace="c:/prog/hmbarea"
gp.overwriteoutput=1
gp.copyfeatures_management("streams/arc", "streams.shp")
gp.Toolbox="Analysis"
gp.buffer("streams.shp", "stbuff200.shp", 200)
print "Finished Buffer.Now Cliping..."
gp.Clip("geol.shp", "stbuff200.shp", "geolstr200.shp")
可以看出成功运行脚本！
然而当在工具中运行时，print 语句不会产生错误，但也不会输出任何东西，因此，我们
用 gp.addmessage("Finished Buffer.Now Cliping...")代替 print 语句。 那么， 如果想无论在工
具中或 Pythonwin 中都可以显示消息，就可以这两句都写上。我喜欢的做法是定义一个
‘sendmsg’函数来输出消息：
def sendmsg(msg):
print msg
gp.addmessage(msg)
….
sendmsg("Finished Buffer.Now Cliping...")
&#61656;  获取工具消息 和 （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"

gp.buffer("stream.shp", "stbuff200.shp", 200)
gp.Clip("geol.shp", "stbuff200.shp", "geolstr200.shp")
except:
print gp.GetMessages()
运行可以看到以下错误提示：
顺便说一下 GetMessages 的集中形式和代表的含
义
② 运行下面简单的代码，并改变 x 的值为非零，查看结果：
x=0.
y=15.
try:
print y/x
except:
print "错误！0 不能是被除数！"
5.2PythonWin  的调试工具
PythonWin 提供了一些工具： ①单步执行代码②插入断点③观测变量， 或其他事情。
这些都可以在 PythonWin 的帮助系统里找到相关示例教程。这里我们演示单步执行代
码的过程。
以前面用过的太阳角计算脚本为例，点击运行之后，在 Debugging 下拉框里看到
调试选项。选择“Step-through in the debugger” 。然
后，一系列调试工具出现 ；一个黄色三角出现
在第一行代码前。 按下 （Step over） 按钮进入下一行， 继续执行； 点击 （Watch）
按钮，出现，在 Expression 下<New Item>里输入 lat 等可以查看当前变量值
； 当你遍历完程序， 点击 关闭调试； 尝试添加一个断点，
当定位到某一行时， 在 点 ， 然后运行此脚
本在“Run in the debugger”下，发现没有任何影响，但是修改第一行为“lat=10” ，运
GetMessages():所以信息
GetMessages(0):只显示消息
GetMessages(1):只显示警告
GetMessages(2):只显示错误

行发现在断点处停止， ，点击 按钮继续程序。得到结果
“Noon sun angle=76.56 Azimuth=0” 。
5.3  地理处理工具 举例
&#61656;  转换脚本（conversion script ）
既然我们已经会在帮助系统里寻找答案， 也掌握了几种程序调试方法， 现在让我们建立
一个脚本。
① 添加引用
import arcgisscripting
② 创建“sendmsg”函数
def sendmsg(msg):
print msg
gp.addmessage(msg)
③ gp.OverwriteOutput=1
④ 使用 Data Management toolbox 中的 Workspace 工具集中的 CreateFolder 工具在
“c:/prog”中创建新文件夹“woodside”
gp.CreateFolder_management("c:\prog","woodside")
⑤ 设置 gp 的工作文件夹为 woodside
gp.workspace="c:\prog\woodside"
⑥ 使用 Conversion 里的 To Raster 工具集中的 DEMtoRaster 工具将 c:/prog/pendata 下
的 woodside.dem 转换为 woodelev 栅格数据，不要提供任何可选参数值
gp.DEMtoRaster_conversion('c:\prog\pendata\woodside.dem','woodelev')
⑦ 在 第 一 行 前 添 加 try: ， 最 后 一 行 后 添 加 except: ， 最 后 加 上 ，
sendmsg(gp.GetMessages()).
&#61656;  要素转换为栅格（feature to raster conversion）
相信有了前面的联系，本例很简单了，将 d:/prog/pendata/landuse(polygon)根据
“LU-CODE”字段转换为“lugrid”栅格数据，分辨率为 30。
try:
import arcgisscripting
def sendmsg(msg):
print msg
gp.addmessage(msg)
gp.OverwriteOutput=1
gp.workspace="c:\prog\pendata"
gp.FeatureToRaster_conversion('landuse/polygon','LU-CODE','lugrid','30')
except:
sendmsg(gp.messages())
6  使用描述（Describe ）和存在（Exists ）获取数据信息
有很多情况下我们需要使用某个 GIS 数据的特征去处理其他数据。 比如， 在栅格运算中，

我们可能想用一个数据的边界去界定另个数据集，很像裁剪操作，或者检查一个数据的
拓扑错误。
6.1  描述
gp.Describe 方法适用于生成一个数据集的若干属性。Geoprocessor Programming Model
中的绿色区域部分， 有按等级划分的很多属性： 所有数据都有一个数据类型 （DataType）
和目录路径（CatalogPath） 。FeatureClass 属性包括所有表个属性和数据集属性；栅格数
据属性包括栅格波段属性和数据集属性，所以数据集属性在一定程度上是共有的
（Feature Class、Rasters、Coverages） 。
&#61656;  利用描述系统用一个数据边界裁剪另一个栅格数据，将裁剪后的土地利用图存在
woodside 文件夹下。
import arcgisscripting
gp = arcgisscripting.create(9.3)
gp.workspace="c:\prog\woodside"
dscRas=gp.Describe("woodelev")
#print dscRas.Extent
envel=str(dscRas.Extent.XMin)+''+str(dscRas.Extent.YMin)+''+str(dscRas.Extent.XM
ax)+' '+str(dscRas.Extent.YMax)

#print envel 因为从 9.3 版本之后，Extent 是一个 Object，所以直接使用
dscRas.Extent 在 clip 语句中会出现错误，而 clip 语句的第二个参数 Rectangle 是一
个 Envelope 类型数据，所以这里我使用了 envel 接收 Extent 的四个值。
try:
gp.clip_management("c:\prog\pendata\lugrid",envel,"luwood")
except:
print gp.GetMessages()
&#61656;  编 写 一 个 脚 本 输 出 “ c:/prog/hmbarea/tmcomp.bil ” 的 波 段 个 数 ， 命 名 为
countBands.py。在栅格数据集属性中寻找这个方法：
import arcgisscripting, sys
gp = arcgisscripting.create(9.3)
gp.workspace = "c:/prog/hmbarea"
dta = "tmcomp.bil"
dsc = gp.Describe(dta)
print dta + " " + dsc.DataType + ": " + str(dsc.BandCount) + " bands."
运行结果： “ tmcomp.bil RasterDataset: 3 bands. ”
6.2  存在（Exists ）
关于一个数据集的一个更基本的信息是其是否存在。 测试一个数据是否存在也
可帮助我们避免错误。 一个非常有用的技术是通过 gp.Exists 判断一个特定的数据是
否存在。
gp.Exists 能用在各种类型数据上， 比如在使用数据之前添加一个判断是否存在
的语句：
if gp.Exists(“streams.shp”)
gp.Buffer(“streams.shp”,”stbuff200.shp”,200)
在下面的代码中我们将使用这种方法判断数据是否存在，其可以成为
overwriteoutput 设置的替代语句。
&#61618;  判断一个字段属性是否存在。上面的方法无法测试一个数据的字段是否存在，但线
面的 listFields 枚举是个技巧：(listFields 返回一个 Python 属性对象列表)
if not gp.listFields(“streams.shp”,”st-code”):
我们将在后面使用这句判断字段是否存在。
6.3  在循环中使用描述和存在
这会让你感到脚本在处理多种数据方面的强大作用。我们将裁剪一整系列的栅
格数据，并存在新文件夹内。在这个处理过程中，我们也会在创建输出数据之前用
“存在工具”中一个很方面的方法检查其是否已经存在。在本例中，我们使用了
ArcInfo Workspace， 并把格网栅格存在这里， 当然你也可以将他们存在 Geodatabase
中。
代码如下：

import arcgisscripting
gp=arcgisscripting.create(9.3)
#gp.checkextension("spatial")
gp.workspace="c:/prog/hmbarea"
if not gp.Exists(gp.workspace+"/cities_polygon.shp"):
gp.FeatureclassToShapefile_conversion("c:/prog/pendata/cities/polygon",gp.workspace)
#因为没有找到 citypoly.shp 和 city 栅格数据，所以我使用 pendata 文件夹下的 cities
里的 polygon 转换的 shp 文件取代。
else:
print gp.workspace+"/cities_polygon.shp exists"
if not gp.Exists(gp.workspace+"/cities"):
gp.PolygonToRaster_conversion("cities_polygon.shp","CITY_CODE","cities","","","60")
else:
print gp.workspace+"/cities exists"
if not gp.Exists(gp.workspace+"hmbcity"):
gp.createArcInfoWorkspace(gp.workspace,"hmbcity")
else:
print gp.workspace+"/hmbcity exists"
try:
gp.mask="cities"
gp.checkoutextension("spatial")#检查有 spatial 的许可
raslist=gp.listrasters()
for ras in raslist:
outputname=gp.workspace+"/hmbcity/"+ras
print outputname
if not gp.exists(outputname):
try:
gp.ExtractByMask_sa(ras,gp.mask,outputname)#用掩膜裁剪
except:
gp.getmessages()
except:
gp.getmessages()
探索 1 ：selSanMateoSel.py
用 Select_Analysis 创建 San Mateo（ 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 很有用，你也可以裁
剪所有数据时设置通配符为“*” ，或者不适用通配符。
代码如下：

import arcgisscripting
gp=arcgisscripting.create(9.3)
gp.overwriteoutput=1
gp.workspace="c:/prog/pendata"
workspaceLoc="c:/prog"
try:
gp.Select_Analysis("cities/polygon","San_Mateo.shp",'"CITY-CODE"=26')
wsName="San_Mateo"
newWS=workspaceLoc+"/"+wsName
if not gp.exists(newWS):
gp.CreateFolder_management(workspaceLoc,wsName)
fList=gp.listdatasets("g*")
for f in fList:
if gp.exists(f+"/polygon"):
gp.Clip_analysis(f+"/polygon","San_Mateo.shp",newWS+"/"+f)
except:
print gp.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”,”V
EG-CODE”]
将每个 coverage 根据对应的字段转换为栅格，分辨率设置为 60。
代码如下：
import arcgisscripting
gp=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=1
if not gp.Exists("multiFeatRas"):#新建一个文件夹保存处理后数据
gp.CreateFolder_management(gp.workspace,"multiFeatRas")
for i in range(len(cov)):#遍历 coverage 数据
if gp.Exists(cov[i]+"/polygon"):
covTop=cov[i]+"/polygon"
elif gp.Exists(cov[i]+"/arc"):
covTop=cov[i]+"/arc"
gp.FeatureToRaster_conversion(covTop,fld[i],"multiFeatRas/"+cov[i]+"g","60")

panda小正太

学习学习！