学习资源: 《面向ArcGIS的Python脚本编程》
英文版更详细 Python Scripting for ArcGIS
练习数据和资源 Esri官方资源

书一共分为四部分，第一部分是Python和地理处理的相关概念（第一章到第四章），第二部分是编写地理处理脚本（第五章到第九章），第三部分是执行地理处理任务（第十章到第十二章），最后一部分是 创建并使用脚本工具（十三到十四章）。核心主要是第二部分。

Part 2 Writing Scripts

Chapter 5 Geoprocessing using Python

引用工具函数的两种方法（推荐第一种）

跳过可选参数的三种方法

设置工作空间和环境

>>> from arcpy import env
>>> env.workspace = "C:/EsriPress/Python/Data/Exercise05"

>>> env.overwriteOutput = True
>>> env.outputCoordinateSystem = spatial_ref

坐标系类（SpatialReference类）

需要实例化之后使用

>>> prjFile = "C:/EsriPress/Python/Data/Exercise05/facilities.prj"
>>> spatial_ref = arcpy.SpatialReference(prjFile)
>>> arcpy.DefineProjection_management("hospitals", spatial_ref)
The result is <Result 'hospitals'>.

>>> print spatial_ref.name
The result is 
NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet

显示工具的提示信息

arcpy.Clip_analysis("parks.shp", "zip.shp", "Results/parks_Clip.shp")
print arcpy.GetMessages()
# 输出所有的提示信息
msgCount = arcpy.GetMessageCount()
print arcpy.GetMessage(msgCount-1)
# 输出指定的提示信息

Chapter 6 Exploring spatial data

使用工具前判断文件是否存在（arcpy.exists）

python自带的函数os.path.exist函数主要是判断文件是否存在，也就是必须要写出完整文件的后缀名，而Arcpy中的exists函数则是可以对不同类别的要素进行判断是否存在。

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/Exercise06"
if arcpy.Exists("cities.shp"):
    arcpy.CopyFeatures_management("cities.shp", "results/cities_copy.shp")

判断数据类型（Describe）

type = arcpy.Describe("cities.shp").shapeType
if type == "Polygon":
    # 需要使用的函数
else:
    print "The feature are not polygons"

数据库中的数据类型非常重要，因为数据库中的户数据元素没有文件扩展名。
文件地理数据库后缀名 .gdb
个人地理数据库后缀名.mdb
企业地理数据库后缀名.sde

>>>import arcpy
>>>element = "C:/Data/study.gdb/final"
>>>desc = arcpy.Describe(element)
>>> arcpy.Describe(element).datasetType
u'FeatureClass'
>>> arcpy.Describe(element).ShapeType
u'Point'

Describe函数返回一个具有很多属性的Describe对象，包含文件路径、目录、路径、文件名称等属性。

>>> mylayer = arcpy.Describe("cities")
>>> mylayer.dataType
#Running the code returns u'FeatureLayer'.
>>> mylayer.shapeType
The result is u'Point'.
----------
>>> mylayer.spatialReference.name
The result is u'GCS_North_American_1983'.
>>> mylayer.spatialReference.type
The result is u'Geographic'.
>>> mylayer.spatialReference.domain
The result is u'-400 -400 400 400'.

列表批处理

import arcpy
from arcpy import env
env.workspace = "C:/Data"
# 列出目录下所有文件名
fclist  =  arcpy.ListFeatureClasses( )
print fclist
# [u'floodzone.shp', u'roads.shp', u'streams.shp', u'wetlands.shp', u'zipcodes.shp']
# 列出以w开头的文件名
fclist  =  arcpy.ListFeatureClasses( "w*" ）
# 列出点要素
fclist = arcpy.ListFeatureClasses( "", "point" )
# 列出栅格数据
rasterlist = arcpy.ListRasters( "", "tif" ）

import arcpy
from arcpy import env
env.workspace = "D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise06"
fc_list = arcpy.ListFeatureClasses()
for fc in fc_list:
    desc = arcpy.Describe(fc)
    print desc.basename + "is a " + desc.shapeType +"feature class."
#另外一种 print "{0} is a {1} feature class".format(desc.basename, desc.shapeType)

列表元组和字典

元组不可修改，元组为()，列表为[]，字典{}，引用是均用[]

cities = ["Austin", "Baltimore", "Cleveland", "Denver"]
states = ["Texas", "Maryland", "Ohio", "Colorado"]
>>>  states[cities.index( "Cleveland" )]

或者直接用字典查询

>>> statelookup = {"Austin": "Texas", "Baltimore": "Maryland", "Cleveland": "Ohio", "Denver": "Colorado"}
>>> statelookup["Cleveland"]
>>> statelookup.key()
>>> statelookup.values()
# 以列表的形式返回
>>>  zoning.items( )
The result is [( 'IND',  'Industrial' ),  ( 'RES',  'Residential' )].

Chapter 7 Manipulating spatial data(Arcpy.da)

游标实际上就是类似SQL的一种查询功能，或者说条件筛选。通过游标搜索得到的记录会输出到一个字段列表中。

• SearchCursor 搜索游标检索行
• InsertCursor 插入行
• UpdateCursor 更新和删除行

搜索游标

import arcpy
from arcpy import env
env.workspace = "D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise07"

fc = "airports.shp"
cursor = arcpy.da.SearchCursor(fc,["NAME"])
for row in cursor:
    print "Airport name = {0}".format(row[0])
del row
del cursor  
# Airport name = New Stuyahok
# Airport name = Ted Stevens Anchorage International
# Airport name = Iliamna
# Airport name = Koyukuk

---

• 使用with，无论游标成功运行还是出现异常，都可以保证数据库锁的关闭和释放，并重新迭代

import arcpy
fc = "C:/Data/study.gdb/roads"
with arcpy.da.SearchCursor( fc, ["STREETNAME"] ) as cursor:
    for row in cursor:
        print "Streetname = {0}".format( row[0] )

---

• SQL语句中，字段格式有(“”)([])，AddFieldDelimiters函数

cursor = arcpy.da.SearchCursor(fc, ["NAME"], '"TOT_ENP" > 100000')
cursor = arcpy.da.SearchCursor(fc,["NAME"],'"NAME" = \'Ketchikan\'')
cursor = arcpy.da.SearchCursor(fc, ["NAME"], '"FEATURE" = \'Seaplane Base\'')
# 注意此时对于条件的书写，不同的文件由不同的格式
delimitedField = arcpy.AddFieldDelimiters(fc, "COUNTY")
cursor = arcpy.da.SearchCursor(fc, ["NAME"], delimitedField + " = 'Anchorage Borough'")

插入游标

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/Exercise07"
fc = "Results/airportsCopy.shp"
cursor = arcpy.da.InsertCursor(fc, "NAME")
cursor.insertRow(["New Airport"])
del cursor

更新游标

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/Exercise07"
fc = "Results/airports.shp"
delimfield = arcpy.AddFieldDelimiters(fc, "STATE")
cursor = arcpy.da.UpdateCursor(fc, ["STATE"], delimfield + " <> 'AK'")
for row in cursor:
    row[0] = "AK"
    cursor.updateRow(row)
del row
del cursor
# 删除
fc = "Results/airports.shp"
cursor = arcpy.da.UpdateCursor(fc, ["TOT_ENP"])
for row in cursor:
    if row[0] < 100000:
        cursor.deleteRow()
del row
del cursor

表和字段

• ValidateTableName( name, {workspace} ）将无效字符替换为下划线(_)

# 将shp文件转移到地理数据库
import arcpy
import os
from arcpy import env
env.workspace = "C:/Data"
outworkspace = "C:/Data/test/study.gdb"
fclist  =  arcpy.ListFeatureClasses( )
for shapefile in fclist:
    fcname  =  arcpy.Describe( shapefile ).basename
    newfcname  =  arcpy.ValidateTableName( fcname )
    outfc = os.path.join( outworkspace, newfcname )
    arcpy.CopyFeatures_management( shapefile,  outfc)

---

• 新增字段检测

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/Exercise07"
fc = "Results/airports.shp"
newfield = "NEW CODE"
fieldtype = "TEXT"
fieldname = arcpy.ValidateFieldName(newfield)
fieldlist = arcpy.ListFields(fc)
fieldnames = []
for field in fieldlist:
    fieldnames.append(field.name)
    if fieldname not in fieldnames:
        arcpy.AddField_management(fc, fieldname, fieldtype, "", "", 12)
        print "New field has been added."
    else:
        print "Field name already exists."

处理文本文件

f = open("D:/Archive01/03Coding/ArcGIS+Python/mytext.txt")
f.read()#读取全部
f.seek(0)
f.read(5)
f.readline()
f.readlines()#以列表的形式返回所有的行
f.close()
# 写入行
f= open("D:/Archive01/03Coding/ArcGIS+Python/mytext.txt","r+")
f.write("Triangulated\nIrregular\nNode" )
f.read()
f.close()
# 修改指定行
f= open("D:/Archive01/03Coding/ArcGIS+Python/mytext.txt","r+")
lines = f.readlines()
lines[2] = "Network"
f.writelines(lines)
# 遍历行
f  =  open( "C:/Data/mytext.txt" )
for line in f.readlines( ):
    <function>
f.close( )

Chapter 8 Working with geometries（处理几何图形）

几何对象的属性
– SHAPE@XY will return a tuple of x,y coordinates representing the feature’s centroid
– SHAPE@LENGTH will return the feature’s length as a double
– SHAPE@ will return the full geometry object

读取几何

import arcpy
from arcpy import env
env.workspace = "D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise08"
# 求周长
fc = "rivers.shp"
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@LENGTH"])
length = 0
for row in cursor:
    length = length + row[0]
units = arcpy.Describe(fc).spatialReference.linearUnitName
print str(length) + " " + units
# 输出XY坐标
fc = "dams.shp"
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@XY"])
for row in cursor:
    x, y = row[0]
    print("{0}, {1}".format(x, y))

读取线和面（single part）

• 线和多边形，访问其中的每一个要素都会返回一个点对象组。想要处理这些数组，需要两次迭代，第一个for循环遍历行对象（即一条线和多边形记录），需要输入OID（对象标识符）字段，用来指明每一个点对象数组的起始位置和结束位置。访问每一个行对象都会获得一个几何对象，该对象是一个点对象数组。getPart函数可以用来获取几何图形的第一部分的点对象数组。第二个for循环用来遍历每一个线或多边形中的点，并输出x、y坐标值。

fc = "rivers.shp"
cursor = arcpy.da.SearchCursor(fc, (["OID@", "SHAPE@"]))
for row in cursor:
    print("Feature {0}: ".format(row[0]))
    for point in row[1].getPart(0):
        print("{0}, {1}".format(point.X, point.Y))

• 注意到getPart()函数索引值为零，只返回几何对象的第一部分，对于singlepart要素类只有第一部分。如果没有索引值，则会返回一个包含所有点对象数组的数组。

多要素（点线面）

• 多个要素只有一个字段表，

# 判断要素类用shapeType属性，Point, Polyline, Polygon, MultiPoint, and MultiPatch（三维数据）几种返回值。
fc = "Hawaii.shp"
print arcpy.Describe(fc).shapeType
# 判断属性
fc = "Hawaii.shp"
cursor = arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"])
for row in cursor:
    if row[1].isMultipart:
        print("Feature {0} is multipart and has {1} parts.".format(row[0],str(row[1].partCount)))
    else:
        print("Feature {0} is single part.".format(row[0]))

• 多边形要素类的xy点

cursor = arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"])
for row in cursor:
    print("Feature {0}: ".format(row[0]))
    partnum = 0
    for part in row[1]:
        print("Part {0}: ".format(partnum))
        for point in part:
            print("{0}, {1}".format(point.X, point.Y))
        partnum += 1

内环

写入几何图形

• 创建线文件

# Write geometries
import fileinput
import string
import os
env.overwriteOutput = True
outpath = "D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise08"
newfc = "Results/newpolyline.shp"
# 创建几何图形
arcpy.CreateFeatureclass_management(outpath, newfc, "Polyline")
# 坐标点文件
infile = "D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise08/coordinates.txt"
# 
cursor = arcpy.da.InsertCursor(newfc, ["SHAPE@"])
array = arcpy.Array()
# 第一种创建线文件
# 将文档中的内容分到ID，X,Y三个字段。
for line in fileinput.input(infile):
    ID, X, Y = string.split(line," ")
    array.add(arcpy.Point(X, Y))
# 插入到新文件中
cursor.insertRow([arcpy.Polyline(array)])
fileinput.close()
del cursor
# 第二种创建面文件
point = arcpy.Point( )
for line in fileinput.input( infile ):
    point.ID, point.X, point.Y = line.split( )
    array.add(point)
polygon = arcpy.Polygon(array)
cursor.insertRow([polygon])
fileinput.close()
del cursor

challenge

• 通过坐标点创建

fc = "newpoly2.shp"
arcpy.CreateFeatureclass_management("C:/Data", fc, "Polygon")
cursor = arcpy.da.InsertCursor(fc, ["SHAPE@"])
array = arcpy.Array()
coordlist =[[0, 0], [0, 1000], [1000, 1000], [1000, 0]]

for x, y in coordlist:
    point = arcpy.Point(x,y)
    array.append(point)
polygon = arcpy.Polygon(array)
cursor.insertRow([polygon])
del cursor

• 统计每个singlepart的周长和面积

fc = "Hawaii.shp"
newfc = "Results/Hawaii_single.shp"
arcpy.MultipartToSinglepart_management(fc, newfc)
spatialref = arcpy.Describe(newfc).spatialReference
unit = spatialref.linearUnitName
cursor = arcpy.da.SearchCursor(newfc, ["SHAPE@"])
for row in cursor:
    print ("{0} {1}".format(row[0].length, unit))
    print ("{0} square {1}".format(row[0].area, unit))

• 创建最大长方形包围圈

fc = "Hawaii.shp"
newfc = "envelope8.shp"
desc = arcpy.Describe(fc)
spatialref = desc.spatialReference
extent = desc.extent
arcpy.CreateFeatureclass_management("D:/Archive01/03Coding/ArcGIS+Python/EsriPress/Python/Data/Exercise08", newfc, "Polygon", "", "", "", spatialref)
cursor = arcpy.da.InsertCursor(newfc, ["SHAPE@"])
array = arcpy.Array()
array.append(extent.upperLeft)
array.append(extent.upperRight)
array.append(extent.lowerRight)
array.append(extent.lowerLeft)
polygon = arcpy.Polygon(array)
cursor.insertRow([polygon])
del cursor

Chapter 9 Working with rasters（处理栅格图形 Arcpy.sa）

ListRasters

rasterlist = arcpy.ListRasters("*","IMG")
for raster in rasterlist:
    print raster

栅格数据基本特征

raster = "tm.img"

desc = arcpy.Describe(raster)

print "Raster base name is: " + desc.basename
print "Raster data type is: " + desc.dataType
print "Raster file extension is: " + desc.extension
print "Raster spatial reference is: " + desc.spatialReference.name
print "Raster format is: " + desc.format
print "Raster compression type is: " + desc.compressionType
print "Raster number of bands is: " + str(desc.bandCount)
print "Raster meanCellHeight is: " + desc.meanCellHeight
print "Raster meanCellWidth is: " + desc.meanCellWidth
# 输出栅格图层的基本属性

处理栅格对象

# 引入栅格数据，之后进行计算
elevraster = arcpy.Raster("elevation")
outraster3 = elevraster * 3.281
outraster3.save("elev_ft")
# 使用工具
outraster = arcpy.sa.Slope("elevation")
outraster2 = arcpy.sa.Aspect("elevation")
outraster.save("Aspect")
# 多图层计算
elev = arcpy.Raster("elevation")
lc = arcpy.Raster("landcover.tif")
slope = Slope(elev)
aspect = Aspect(elev)
goodslope = ((slope > 5) & (slope < 20))
goodaspect = ((aspect > 150) & (aspect < 270))
goodland = ((lc == 41) | (lc == 42) | (lc ==43))
outraster = (goodslope & goodaspect & goodland)
outraster.save("C:/EsriPress/Python/Data/Exercise09/Results/final")

转载自：https://blog.csdn.net/leeloo666/article/details/72539364