ArcPy实验代码抽查(11-3)

ArcPy实验代码抽查

小提醒

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

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

数据分享

7-1 市区择房分析

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 7-1 市区择房分析.py

import arcpy
import os

source_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')  # 获取资源路径
target_path = source_path + '\\result_py'
os.mkdir(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 最短路径问题分析与应用

抽到了我也不会o(╥﹏╥)o

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 7-2 最短路径问题分析与应用.py

import arcpy
import os
import shutil
import sys

reload(sys)
sys.setdefaultencoding('utf8')  # 解决编码错误
print "输入路径示例：D:\\桌面\\应用数据\\Chp7\\Ex2"
path = raw_input("请输入数据所在路径：").decode("utf-8")
place = "city.mdb\\city\\place"
arcpy.env.workspace = path
point_path = path + "\\points"
result_path = path + "\\result"
if os.path.exists(point_path):  # 测试用，自动删除前面运行时所建的文件夹
    shutil.rmtree(point_path)
    os.mkdir(point_path)
else:
    os.mkdir(point_path)
if os.path.exists(result_path):
    shutil.rmtree(result_path)
    os.mkdir(result_path)
else:
    os.mkdir(result_path)
reName = 1  # 目的地重复时，命名用
name = []  # 判断目的地是否重复时，使用
group = 1  # 输出图层组时，使用
destinations = [row[0] for row in arcpy.da.SearchCursor(place, "NAME")]  # 目的地的名称
destinations.remove("HOME")  # 删除出发点
for destination in destinations:
    points = "HOME_{}.shp".format(destination)
    if points in name:  # 因为海星超市有三个，所以的检验，重复的在名字后面加1
        points = "HOME_{}{}.shp".format(destination, reName)
        reName += 1
    name.append(points)
    arcpy.Select_analysis(place, point_path + "\\" + points, "[NAME] = 'HOME' or [NAME] = '{}' and [OBJECTID]= {}".format(destination, group))
    arcpy.TraceGeometricNetwork_management('city.mdb\\city\\city_Net', 'city_Net_{}'.format(group), point_path + "\\" + points, "FIND_PATH", "", "", "", "", "", "NO_TRACE_ENDS", "NO_TRACE_INDETERMINATE_FLOW", "", "", "AS_IS", "", "", "", "AS_IS")  # Process: 追踪几何网络
    arcpy.CopyFeatures_management("city_Net_{}\\net".format(group), result_path+'\\' + points)  # Process: 复制要素
    group += 1
print "ok!"
# 不知道为什么，HOME_海星超市.shp 这个是为空の，点感觉没问题吖！就是输出为空...
# 有解决到这个问题的大佬，记得告诉我，我是实在解决不了o(╥﹏╥)o

8-2 寻找最佳路径

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# 8-2 寻找最佳路径.py

import arcpy
from arcpy import Raster
from arcpy.sa import *
import os

# 重分类——相等间隔
def divide(raster, times=10):
    minNum = round(float(str(arcpy.GetRasterProperties_management(raster, "MINIMUM", ""))), 6)
    maxNum = round(float(str(arcpy.GetRasterProperties_management(raster, "MAXIMUM", ""))), 6)
    total = maxNum - minNum
    interval = total / times
    class_table = ''
    end = minNum
    for i in range(1, times + 1):
        start, end = end, end + interval
        class_table += "{} {} {};".format(start, end, i)
    class_table = class_table[:-1]
    return class_table

source_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
os.mkdir(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 熊猫密度分布图.py

import arcpy
from arcpy.sa import *
import os

def set_extent(feature):
    # 左下右上
    arcpy.env.extent = feature
    extent = 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.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区域分布图的生成与对比.py

import arcpy
from arcpy.sa import *
from arcpy import Raster
import os

source_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.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 山顶点提取.py

import arcpy
from arcpy.sa import *
from arcpy import Raster
import os

source_path = raw_input("请输入数据所在的文件绝对路径：").decode('utf-8')
target_path = source_path + '\\result_py'
os.mkdir(target_path)
dem = source_path + '\dem'
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: 栅格转点
posted @   槑孒  阅读(219)  评论(0)  编辑  收藏  举报
编辑推荐：
· AI与.NET技术实操系列（二）：开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识（上）
· 浏览器原生「磁吸」效果！Anchor Positioning 锚点定位神器解析
阅读排行：
· 全程不用写代码，我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· 记一次.NET内存居高不下排查解决与启示
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布：重大改进与新特性概览！
点击右上角即可分享
微信分享提示
AI FOR CODE 大赛