ArcGIS空间分析实验教程-部分实验案例(ArcPy实现)

Tips

如果是直接复制代码运行，请注意查看以下提醒：

• 1.本代码运行时，需要输入对应资源所在的绝对路径，均只需输入至...\Chp7\Ex1,如：D:\桌面\ArcPy实验代码抽查\Chp7\Ex1。
• 2.# -*- coding: utf-8 -*- 这一行注释不能删除，否则会报编码错误。

数据分享

7-1 市区择房分析

# -*- coding: utf-8 -*-
import arcpy
print "输入示例：D:\桌面\Chp7\Ex1"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')  # 获取资源路径
arcpy.Buffer_analysis('Marketplace.shp', 'Marketplace_Buff', "YUZHI_")  # Process: Marketplace.shp缓冲区
arcpy.Buffer_analysis('school.shp', 'school_Buff', "750 Meters")  # Process: school.shp缓冲区
arcpy.Buffer_analysis('famous place.shp', 'famous_place_Buff', "500 Meters")  # Process: famous place.shp缓冲区
arcpy.Intersect_analysis(['Marketplace_Buff.shp', 'school_Buff.shp', 'famous_place_Buff.shp'], 'Intersect')  # Process: 相交
arcpy.Select_analysis('network.shp', 'network_select', "TYPE = 'ST'")  # Process: 筛选字段为ST的记录
arcpy.Buffer_analysis('network_select.shp', 'network_select_Buff', "200 Meters")  # Process: network.shp缓冲区
arcpy.Erase_analysis('Intersect.shp', 'network_select_Buff.shp', '适宜区')  # Process: 擦除

7-2 最短路径问题分析与应用

不确定是否正确喔~

# -*- coding: utf-8 -*-
import arcpy
print "输入示例：D:\桌面\Chp7\Ex2"
arcpy.env.workspace = raw_input("请输入数据所在路径：").decode("utf-8")
arcpy.TraceGeometricNetwork_management('city.mdb\\city\\city_Net', 'city_Net_t', 'city.mdb\\city\\place', "FIND_PATH")  # Process: 追踪几何网络
arcpy.CopyFeatures_management("city_Net_t\\net", '几何路径')  # Process: 复制要素

8-2 寻找最佳路径

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp8\Ex2"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
a = Reclassify(r'Road.mdb\river', "Value", "0 1;1 2;2 5;3 8;4 10")  # Process: river重分类
b = Slope('Road.mdb\dem')  # Process: dem坡度
c = Reclassify(b, "Value", '0.0 5.730143 1;5.730143 11.460286 2;11.460286 17.190429 3;17.190429 22.920572 4;22.920572 28.650715 5;28.650715 34.380858 6;34.380858 40.111001 7;40.111001 45.841144 8;45.841144 51.571287 9;51.571287 57.30143 10') # Process: Slope_dem重分类
d = FocalStatistics('Road.mdb\dem', "Rectangle 3 3 CELL", "MEAN")  # Process: dem焦点统计
e = Reclassify(d, "Value", '179.462219 210.7653318 1;210.7653318 242.0684446 2;242.0684446 273.3715574 3;273.3715574 304.6746702 4;304.6746702 335.977783 5;335.977783 367.2808958 6;367.2808958 398.5840086 7;398.5840086 429.8871214 8;429.8871214 461.1902342 9;461.1902342 492.493347 10')  # Process: FocalSt_dem重分类
f = a + c * 0.6 + e * 0.4  # Process: 栅格计算
g = CostDistance('Road.mdb\startPot.shp', f, "", 'backlink')  # Process: 成本距离
CostPath('Road.mdb\endPot.shp', g, 'backlink').save('成本路径')  # Process: 成本路径

8-3 熊猫密度分布图

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp8\Ex3"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
xm = 'Xmpoint.shp'
arcpy.env.extent = '15166438 -134313 15220978 -97358'  # 设置处理范围
a = EucAllocation(xm, "5000", "", "500", "ID")  # Process: source_path欧氏分配
arcpy.AddField_management(a, "Area", "LONG")  # Process: EucAllo_shp添加字段Area
arcpy.CalculateField_management(a, "Area", "[COUNT] * 500 *500")  # Process: EucAllo_shp计算字段Area
arcpy.JoinField_management(xm, "ID", a, "Value", "Area")  # Process: EucAllo_shp连接字段
arcpy.AddField_management(xm, "power", "LONG")  # Process: Xmpoint.shp添加字段power
arcpy.CalculateField_management(xm, "power", "3.1415927*5000*5000 / [Area]")  # Process: Xmpoint.shp计算字段power
b = KernelDensity(xm, "power", "500")  # Process: Xmpoint.shp核密度分析
Times(b, 10000000).save("熊猫分布图")  # Process: 乘

8-4 GDP区域分布图的生成与对比

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp8\Ex4"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
arcpy.env.mask = r"GDP.mdb\bound.shp"  # 设置掩膜
GDP = 'GDP.mdb\GDP.shp'
IDW2 = Idw(GDP, "GDP", "500", "2")  # Process: 反距离权重法
IDW5 = Idw(GDP, "GDP", "500", "5")  # Process: 反距离权重法 (2)
Abs(IDW2 - IDW5).save('IDW2_IDW5')  # Process: 栅格计算
Spr0 = Spline(GDP, "GDP", "500", "REGULARIZED", "0")  # Process: 样条函数法
Spr01 = Spline(GDP, "GDP", "500", "REGULARIZED", "0.01")  # Process: 样条函数法 (2)
Abs(Spr0 - Spr01).save('Spr0_Spr01')  # Process: 栅格计算 (2)
Spt0 = Spline(GDP, "GDP", "500", "TENSION", "0")  # Process: 样条函数法 (3)
Spt5 = Spline(GDP, "GDP", "500", "TENSION", "5")  # Process: 样条函数法 (4)
Abs(Spt0 - Spt5).save('Spt0_Spt5')  # Process: 栅格计算 (3)
Abs(IDW2 - Spr01).save('IDW2_Spr01')  # Process: 栅格计算 (4)

8-5 山顶点提取

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp8\Ex5"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
a = FocalStatistics('dem', "Rectangle 21 21 CELL", "MAXIMUM")  # Process: 焦点统计
b = Test(a - 'dem', 'value=0')  # Process: 条件测试
c = Con(b, 1)  # Process: 条件函数
arcpy.RasterToPoint_conversion(c, '山顶点')  # Process: 栅格转点

11-1 山脊线提取

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp11\Ex1"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
zd = Test('dem' - FocalStatistics('dem', "Rectangle 11 11 CELL", "MEAN"), 'value>0')  # 求正地形
a = Fill('dem')  # 填洼
b = FlowDirection(a)  # 流向
c = FlowAccumulation(b)  # 流量
d = Test(c, 'value=0')  # 条件测试
Con(zd*d, 1).save("山脊线")  # 条件函数

11-1 山谷线提取

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp11\Ex1"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
fd = Test('dem' - FocalStatistics('dem', "Rectangle 11 11 CELL", "MEAN"), 'value<0')  # 求负地形
a = Abs(Minus(5000, "dem"))  # 作差取绝对值
b = FlowDirection(a)  # 流向
c = FlowAccumulation(b)  # 流量
d = Test(c, 'value=0')  # 条件测试
Con(fd*d, 1).save("山谷线")  # 条件函数

11-2 地形鞍部点提取

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp11\Ex2"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
# 变量z(正)开头为山脊线提取的前部分步骤，变量f(负)开头为山谷线提取的前部分步骤
zdx = Test('dem' - FocalStatistics('dem', "Rectangle 11 11 CELL", "MEAN"), 'value>0')  # 求正地形
za = Fill('dem')  # 填洼
zb = FlowDirection(za)  # 流向
zc = FlowAccumulation(zb)  # 流量
zd = Test(zc, 'value=0')  # 条件测试
fa = Abs(Minus(5000, "dem"))  # 作差取绝对值
fb = FlowDirection(fa)  # 流向
fc = FlowAccumulation(fb)  # 流量
fd = Test(fc, 'value=0')  # 条件测试
c = Con(zd*fd*zdx, 1)  # 条件函数
arcpy.RasterToPoint_conversion(c, "鞍部点")  # 栅格转点

11-3 河谷网络的提取

# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
print "输入示例：D:\桌面\Chp11\Ex3"
arcpy.env.workspace = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')  # 工作空间
a = Fill('dem')  # 填洼
b = FlowDirection(a)  # 流向
c = FlowAccumulation(b)  # 流量
d = Test(c, "value>100")  # 条件测试
e = Con(d, 1)  # 条件函数
StreamToFeature(e, b, '河谷')  # 栅格河网矢量化