ArcPy实验代码抽查

小提醒

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

• 1.本代码运行时，需要输入对应资源所在的绝对路径，均只需输入至...\Chp7\Ex1,如：D:\桌面\ArcPy实验代码抽查\Chp7\Ex1
• 2.在运行代码之前，需要确定输入的绝对路径下是否存在result_py文件夹，若存在，则将其删除
• 3.本文代码均已测试，若存在问题，请评论留言或私信，谢谢！

数据分享

7-1 市区择房分析

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 7-1 市区择房分析.pyimport arcpy
import ossource_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')  # 获取资源路径
target_path = source_path + '\\result_py'
os.mkdir(target_path)
arcpy.env.scratchWorkspace = target_path
arcpy.env.workspace = target_path
arcpy.Buffer_analysis(source_path + '\Marketplace.shp', 'Market_Buffer', "YUZHI_", "FULL", "ROUND", "ALL", "", "PLANAR")  # Process: Marketplace.shp缓冲区
arcpy.Buffer_analysis(source_path + '\school.shp', 'school_Buffer', "750 Meters", "FULL", "ROUND", "ALL", "", "PLANAR")  # Process: school.shp缓冲区
arcpy.Buffer_analysis(source_path + r'\famous place.shp', 'famous_Buffer', "500 Meters", "FULL", "ROUND", "ALL", "", "PLANAR")  # Process: famous place.shp缓冲区
arcpy.Intersect_analysis(['Market_Buffer.shp', 'school_Buffer.shp', 'famous_Buffer.shp'], 'intersect', "ALL", "", "INPUT")  # Process: 相交
arcpy.Select_analysis(source_path + r'\network.shp', 'network_Sel', "TYPE = 'ST'")  # Process: 筛选字段为ST的记录
arcpy.Buffer_analysis('network_Sel.shp', 'net_S_Buff200', "200 Meters", "FULL", "ROUND", "ALL", "", "PLANAR")  # Process: network.shp缓冲区
arcpy.Erase_analysis('intersect.shp', 'net_S_Buff200.shp', '适宜区', "")  # Process: 擦除

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

略…(不会诶)

8-2 寻找最佳路径

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-2 寻找最佳路径.pyimport arcpy
from arcpy import Raster
from arcpy.sa import *
import os# 重分类——相等间隔
def divide(raster, times=10):MinMUM = float(str(arcpy.GetRasterProperties_management(raster, "MINIMUM", "")))MamMUM = float(str(arcpy.GetRasterProperties_management(raster, "MAXIMUM", "")))minNum = round(MinMUM, 6)maxNum = round(MamMUM, 6)total = maxNum - minNuminterval = total / timesclass_table = ''end = minNumfor i in range(1, times + 1):start, end = end, end + intervalclass_table += "{} {} {};".format(start, end, i)class_table = class_table[:-1]return class_tablesource_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
os.mkdir(target_path)
arcpy.env.scratchWorkspace = target_path
arcpy.env.workspace = target_path
Reclassify(source_path + r'\Road.mdb\river', "Value", "0 1;1 2;2 5;3 8;4 10", "DATA").save('Reclass_rive')  # Process: river重分类
Slope(source_path + '\Road.mdb\dem', "DEGREE", "1", "PLANAR", "METER").save('Slope_dem')  # Process: dem坡度
Reclassify('Slope_dem', "Value", divide("Slope_dem"), "DATA").save('Reclass_Slop')  # Process: Slope_dem重分类
FocalStatistics(source_path + '\Road.mdb\dem', "Rectangle 3 3 CELL", "MEAN", "DATA", "90").save('FocalSt_dem')  # Process: dem焦点统计
Reclassify('FocalSt_dem', "Value", divide('FocalSt_dem'), "DATA").save('Reclass_Foca')  # Process: FocalSt_dem重分类
(Raster('Reclass_rive') + Raster('Reclass_Slop') * 0.6 + Raster('Reclass_Foca') * 0.4).save('rastercacu')  # Process: 栅格计算
CostDistance(source_path + '\Road.mdb\startPot.shp', 'rastercacu', "", 'huisuo', "", "", "", "", "").save('CostDis_star')  # Process: 成本距离
CostPath(source_path + '\Road.mdb\endPot.shp', 'CostDis_star', 'huisuo', "EACH_CELL", "Id", "INPUT_RANGE").save('成本路径')  # Process: 成本路径

8-3 熊猫密度分布图

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-3 熊猫密度分布图.pyimport arcpy
from arcpy.sa import *
import osdef set_extent(feature):# 左下右上arcpy.env.extent = featureextent = str(arcpy.env.extent).split()[:4]return '{} {} {} {}'.format(float(extent[0]) - 5000, float(extent[1]) - 5000, float(extent[2]) + 5000, float(extent[3]) + 5000)source_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
Xmpoint = source_path + '\Xmpoint.shp'  # 因为多次用到这个文件，所以把路径定义为常量
os.mkdir(target_path)
arcpy.env.scratchWorkspace = target_path
arcpy.env.workspace = target_path
arcpy.env.extent = set_extent(Xmpoint)  # 设置处理范围
# Process: 检查属性表字段是否含有"area", "power"，有则删除
fields = arcpy.ListFields(Xmpoint)
for field in fields:if field.name.lower() in ["area", "power"]:arcpy.DeleteField_management(Xmpoint, field.name)
EucAllocation(Xmpoint, "5000", "", "500", "ID", '', '', "PLANAR", "", '').save('EucAllo_shp')  # Process: source_path欧氏分配
arcpy.AddField_management('EucAllo_shp', "Area", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")  # Process: EucAllo_shp添加字段Area
arcpy.CalculateField_management('EucAllo_shp', "Area", "[COUNT] * 500 *500", "VB", "")  # Process: EucAllo_shp计算字段Area
arcpy.JoinField_management(Xmpoint, "ID", 'EucAllo_shp', "Value", "Area")  # Process: EucAllo_shp连接字段
arcpy.AddField_management(Xmpoint, "power", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")  # Process: Xmpoint.shp添加字段power
arcpy.CalculateField_management(Xmpoint, "power", "3.1415927*5000*5000 / [Area]", "VB", "")  # Process: Xmpoint.shp计算字段power
KernelDensity(Xmpoint, "power", "500", "", "SQUARE_MAP_UNITS", "DENSITIES", "PLANAR").save('KernelD')  # Process: Xmpoint.shp核密度分析
Times('KernelD', 10000000).save("熊猫分布图")  # Process: 乘

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

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-4 GDP区域分布图的生成与对比.pyimport arcpy
from arcpy.sa import *
from arcpy import Raster
import ossource_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
os.mkdir(target_path)
GDP = source_path+'\GDP.mdb\GDP.shp'
arcpy.env.scratchWorkspace = target_path
arcpy.env.workspace = target_path
arcpy.env.mask = source_path + r"\GDP.mdb\bound.shp"  # 设置掩膜
Idw(GDP, "GDP", "500", "2", "VARIABLE 12", "").save('IDW2')  # Process: 反距离权重法
Idw(GDP, "GDP", "500", "5", "VARIABLE 12", "").save('IDW5')  # Process: 反距离权重法 (2)
Abs((Raster('IDW2') - Raster('IDW5'))).save('IDW2_IDW5')  # Process: 栅格计算
Spline(GDP, "GDP", "500", "REGULARIZED", "0", "12").save('Spr0')  # Process: 样条函数法
Spline(GDP, "GDP", "500", "REGULARIZED", "0.01", "12").save('Spr01')  # Process: 样条函数法 (2)
Abs((Raster('Spr0') - Raster('Spr01'))).save('Spr0_Spr01')  # Process: 栅格计算 (2)
Spline(GDP, "GDP", "500", "TENSION", "0", "12").save('Spt0')  # Process: 样条函数法 (3)
Spline(GDP, "GDP", "500", "TENSION", "5", "12").save('Spt5')  # Process: 样条函数法 (4)
Abs((Raster('Spt0') - Raster('Spt5'))).save('Spt0_Spt5')  # Process: 栅格计算 (3)
Abs((Raster('IDW2') - Raster('Spr01'))).save('IDW2_Spr01')  # Process: 栅格计算 (4)

8-5 山顶点提取

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-5 山顶点提取.pyimport arcpy
from arcpy.sa import *
from arcpy import Raster
import ossource_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
os.mkdir(target_path)
dem = source_path + '\dem'
arcpy.env.scratchWorkspace = target_path
arcpy.env.workspace = target_path
arcpy.env.extent = dem
arcpy.env.cellSize = dem
arcpy.env.mask = dem
Contour(dem, 'Contour_dem15', "15", "0", "1", "CONTOUR", "")  # Process: 15m等值线
Contour(dem, 'Contour_dem75', "75", "0", "1", "CONTOUR", "")  # Process: 75m等值线 (2)
Hillshade(dem, "315", "45", "NO_SHADOWS", "1").save('HillSha_dem')  # Process: 山体阴影
FocalStatistics(dem, "Rectangle 21 21 CELL", "MAXIMUM", "DATA", "90").save('FocalSt_dem')  # Process: 焦点统计
Test((Raster('FocalSt_dem') - Raster(dem)), 'value=0').save('F_dem')  # Process: 栅格计算
Con('F_dem', '1', '', "Value=1").save('raster_peak')  # Process: 条件函数
arcpy.RasterToPoint_conversion('raster_peak', '山顶点', "VALUE")  # Process: 栅格转点