资源描述
精品文档
Python脚本使用详解
目录
写在前面的话 2
前言 2
一、PYTHON语言基础 3
1数学运算符 3
2字符串操作 4
3模块的使用(Modules) 5
4使用def构建函数 6
5流程控制结构:If,While,For 6
6简单输入和输出 9
二、ARCGIS&PYTHON 10
1如何创建地理处理对象(geoprocessor object) 10
2获取地理处理帮助 10
2.1举例:如何使用Geoprocessor Programming Model中的Lists 11
3使用地理处理工具——Toolboxes和Aliases 11
4在建模中使用脚本(Scripts in ModelBuilder) 12
5 在PythonWin里调试地理处理脚本 18
5.1 调试选择和消息 19
5.2PythonWin的调试工具 20
5.3地理处理工具举例 21
6使用描述(Describe)和存在(Exists)获取数据信息 21
6.1描述 22
6.2存在(Exists) 23
6.3在循环中使用描述和存在 23
7在Python脚本中使用地图代数(Map Algebra) 26
8数据管理和指针(Data Management and Cursors) 27
8.1数据管理(Data Management) 27
8.2指针(Cursors) 28
附录1:地理处理脚本中输入&输出方法指南 31
附录2:其他 32
写在前面的话
一直想学习ArcGIS中的Python脚本,大四下半学期终于有了时间,可是想找到这么一本好的教材不容易。茫茫互联网,终于找到了旧金山州立大学Jerry Davis教授的个人主页,对其中《Geoprocessing Scripts With Python》如获至宝,独乐乐不如众乐乐,现在将其教程翻译并结合自己的学习情况给出总结。希望能够给更多想学习Python的同学一个参考。
另外,在我刚开始接触Python时,是看了台湾辅仁大学一位老师的视频课件,在此致谢。
我想从两个大部分总结:一、Python语言基础;二、ArcGIS&Python。其中第一部分参考了《Python精要参考(第二版)》、《Python编程金典(读书笔记)》等书籍文献。对于多数读者来说,可能或多或少有一些编程基础,所以理解起来应该不成问题。
文中多数数据来自Jerry Davis教授的主页,放在“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中引入了全新的Python Window来增强内嵌的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
Print newaz%360
90
取模的用途之一——方位角加180后逆转方向
5**2
25
25**0.5
5.0
没有sqrt()功能,除非添加math模块
2字符串操作
注:使用Python帮助:有超过30种内置方法来处理字符,请到Sequence Types下的String Methods寻找帮助!
字符串是一串字母,比如’San Francisco’,字符串下标从0开始。学习字符串语法的最好方法是自己动手尝试,下标展示之:
输入
结果
Notes
print 'zhulj'.capitalize()
Zhulj
s.capitalize()即将capitalize()方法用于s
s='zhulj'
print s.capitalize()
print s[0]
z
Strings可以像一个字母列表一样处理,第一个字母下标为0,某个字符段可以用1:3来格式化:从第1个的开头到第3个的开头,不包括下标为3的字母;s[-1]表示倒过来第一个,相当于s[len(s)-1]
s1=s[1]
print s1
h
print s[-2:]
lj
print s[2:3]
u
print s[2:4]
ul
print s[2:],s[:5]
ulj zhulj
s2=s.upper()
print s2
ZHULJ
我们可以将字符串方法的结果赋给新的变量
s3=s+s2
print s3
zhuljZHULJ
字符串组合用“+”
print s*3
zhuljzhuljzhulj
字符串重复用“*”,后为重复次数
selstr='"elev">1000'
print selstr
"elev">1000
字符串可以使用单引号或双引号,跨行时用双引号。
othersel=”’elev’>1000”
print othersel
‘elev’>1000
print s.isupper()
False
一些方法返回值为布尔型(True或False),一些返回索引值(下标值)
print s2.isupper()
True
p='d:/work/lu.shp'
print p.find(‘.’)
10
print p.find(‘/’)
2
plist=p.split('/')
print plist
['d:', 'work', 'lu.shp']
你可以用split()方法解析出不同的字符串片段,并创建一个列表(List),我们可以使用其中不同的元素
print plist[0]
d:
print plist[1]
work
p2='d:\\work\\soil.shp'
print p2
d:\work\soil.shp
反斜线“\”和某些字母一起有特殊用法,如\n为换行,“\”为转义字符,如“\\”则表示“\”
print 'Jerry\'s Kids'
Jerry's Kids
print 'Jerry\'s\nKids'
Jerry's
Kids
p3=r'd:\work\soil.shp'
print p3
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 <Module名>放在程序顶部,比如:
import arcgisscripting
当然,这不必成为你程序的第一行,但必须在使用它里面方法之前。
当要引用多个模块是,中间用逗号分隔,比如:
import arcgisscripting,sys,string,os,math
我们也可以自己为频繁使用的方法创建Module,下面,我们开始体验内置的Modules。
math和random模块
很多常用的数学计算功能都可以通过math找到,比如三角计算或对数计算,如果要使用复杂数字,就使用cmath模块。
和之前一样,通过以下表格来体现模块的使用:
输入
结果
Notes
import math
print math.log10(100)
2.0
以10为底的对数
print math.log(100)
4.60517018599
自然对数
print math.pi
3.14159265359
π是一个静态常量,所以不需要括弧
pi=math.pi
print pi
3.14159265359
如果不想总是输入“math.pi”可以将其赋给一变量
pi
3.1415926535897931
不需要print即可查看变量值
print math.sin(radians)
print math.cos(radians)
print math.tan(radians)
三角函数的计算是弧度制
degrad=pi/180
45*degrad
0.78539816339744828
度转化为弧度
sin=math.sin
sin(45*degrad)
sin(90*degrad)
0.70710678118654746
1.0
即使功能函数(像sin)都可以赋给一个变量
math.e
2.7182818284590451
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.28468923956189
不同的方式,相同的结果
import random
random.random()
0.27281588185756478
random()方法,每次结果都不同,值域为[0.0,1.0)
rnd=random.random
rnd()
0.4456393835072503
mu=50
s=10
print random.gauss(mu,s)
46.528202194
4使用def构建函数
有点像Module,但更简单,函数是一个自己定义功能,用在之后的代码中,并且提供任何你想要使用的参数。这个函数从此可像变量那样在程序中使用,结合例子更容易理解。接下来的代码定义了一个将度转换为弧度的简单函数,同时也定义了一个弧度转换为度的函数,它们和Excel内置的函数类似。
import math
def radians(angdeg):
return angdeg*math.pi/180
def degrees(angrad):
return angrad*180/math.pi
print math.sin(radians(45))
print degrees(math.acos(0.5))
运行之,得到结果:0.707106781187 60.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
if sunangle>90:
sunangle=180-sunangle
azimuth=0
print sunangle,azimuth
上面的例子中lat和decl强制赋了值。
有三种流程控制操作:
if 仅在一个特定情形下才执行语句;
while 当一种情形存在下,持续执行语句
for 遍历一系列值
这些语法和def有些相似:初始语句后加顿号、需要执行的语句块有缩进。
这三个结构的一些重要的公共特征:
①if、while、for语句均以冒号结尾,接下来是缩进的代码块,用于if、while、for定义的情形。在脚本编写窗口,你会发现,你在一行末尾打上冒号后,下一行自动缩进,在接下来的一行按下退格键取消缩进。
②如果你只需做一件事情,你可以在冒号后面同一行添加简短的语句,比如:
if x>0: print ‘x 比0大’
print ‘下一行不要缩进了。。。’
if(continued)
接下来,我们会探索一下另一个方便的模块:os.path:
开始之前,在d:/下创建一个“testfolder”文件夹,然后新建一个“test.txt”文件;
尝试以下代码段,确保print语句前有缩进。
import os.path
if os.path.exists("d:/testfolder/test.txt"):
print "测试文件夹存在"
print "txt文件存在"
elif os.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”代表米,所以可以通过读取这个文件判断。
下面的脚本演示了上述内容:
import fileinput
infile=r"c:\prog\pendata\woodside.dem"
firstline=fileinput.input(infile)[0]
unitchar=firstline[539]
unit="(unknown:not a 7.5' DEM?)"
if unitchar=="1":unit="feet"
if unitchar=="2":unit="meters"
print "\nElevation in"+" "+unit
fileinput.close()
输出结果:Elevation in feet
while(continued)
Ø 运行下面的代码,说明了一种while循环:
x=1
while x<10:
print x
x=x+1
屏幕依次输出1~9
Ø 下面说明一下“==”(等于)的概念:
x=5
z=x==4
print z
输出False
x=5
z=x==5
print z
输出True
“==”是逻辑运算符之一,其他有“<”(小于)、“>”(大于)、“>=”(大于或等于)、“<=”(小于或等于)、“<>”(不等于)。
使用逻辑运算符计算得到结果为布尔型:true(1)和false(0)。
下面例子简单体会一下布尔型表达式:
x=1
while x<10:
print x
x=x+1
表达式“x<10”结果是true或false,所以这样允许我们在计算完一种情况时运行一系列代码。许多情况下我们需要使用条件代码。
while循环的一个优点是允许我们跳过整个部分,如果条件不满足初始情况。由于while提供一种容易结束循环的方法,我们甚至用它代替if语句。当循环一个数据集时(GIS中很常用的工作)while循环很有用。在后面地理处理中我们会接触一些例子。
for
Ø 尝试下面代码,演示了for循环:
for x in [1,2,3,4]:(注:[1,2,3,4]用range(1,5)代替是一样的)
msg="hello world"
print str(x)+" "+msg(注:当我们希望一个数字x和字符串连接时,必须先对数字进行格式转换即str(x),否则系统报错)
Ø 下面的代码创建并输出指定文件夹内shp文件名列表(每个都以‘.shp’结尾)
import os
ws="c:/prog/hmbarea"
ilist=os.listdir(ws)#创建一个列表保存工作文件夹内的文件
fcs=[]#创建一个空列表,保存结尾为‘.shp’的文件
for i in ilist:
if i.endswith(".shp"):
fcs.append(i)
for fc in fcs:
print fc(输出结果如右图所示)
Ø 下面这个例子的循环较多次数。变量mu是算术平均数,s是标准差——这两个是random.gauss()用到的参数,可以尝试改变n的值查看结果!
import random
mu=50
s=10
z5=mu-s*1.96
z95=mu+s*1.96
count=0
n=100000
for i in range(n):
x=random.gauss(mu,s)
if x<z5 or x>z95:count=count+1
print float(count)/n(每次运行的结果都不同,但都在0.05左右,即统计结果在5%左右)
6简单输入和输出
我们现在利用前面计算太阳角代码的片段,之前是直接指定参数值,现在我们有很多种方法提供输入参数,现在我们用sys方法。尝试下面的代码,点击(run运行)之后,在对话框中函数自变量(Arguments)一栏填入:40 23.44,如下图:
import sys
lat=float(sys.argv[1])
decl=float(sys.argv[2])#使用arguments(argv)方法从“Arguments”一栏中获取输入参数,并指定一个浮点型转换将字符型输入值传递给lat和decl
sunangle=90-lat+decl
azimuth=180
if sunangle>90:
sunangle=180-sunangle
azimuth=0
print "正午太阳角="+str(sunangle)
print "方位角="+str(azimuth)(结果:正午太阳角=73.44 方位角=180)
二、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,填写如下图:
l 注意:脚本文件是一个脚本工具的参考!非常重要的一点,你使用脚本创建了一个脚本工具,但是脚本本身并没有和工具一起保存,脚本工具作为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,
展开阅读全文