arcpy批量生成可视域专题图

Dust_187J3X1

已于 2022-06-14 19:06:32 修改

阅读量619

点赞数 1

分类专栏： ArcGIS Arcpy 文章标签： python

于 2022-06-14 18:45:28 首次发布

本文链接：https://blog.csdn.net/qq_42072638/article/details/125283821

版权

ArcGIS 同时被 2 个专栏收录

3 篇文章 1 订阅

订阅专栏

Arcpy

2 篇文章 0 订阅

订阅专栏

arcpy批量生成可视域专题图

思路

使用 arcpy + data driven page ，针对每个观测点生成可视域专题图。

操作步骤

Step0.下载项目文件

项目文件地址
项目目录

Step1.制作观测点数据

向 ./viewshed/points/points.shp 中录入数据。

其中，distance代表观测距离（m），tower_h代表观测点相对于地面的高度，name为名称，longitude经度和latitude纬度需自行计算或输入。默认使用name字段作为导出的专题的名称。

Step2.制作DEM数据

下载并拼接DEM数据。可以在地理空间数据云中获取数据。
地理空间数据云下载DEM数据
拼接完成后替换./viewshed/dem目录下的所有文件。也可以在第三步中配置dir_dem目录（默认使用目录中的第一份数据）。

Step3.配置参数

编辑可视域.py文件，将根目录base_dir、点位目录dir_point、DEM目录dir_dem修改为自己的目录。避免转译，最好将目录中的斜杠\替换为反斜杠/。

import arcpy
from arcpy.sa import *
import os
import logging
import sys 
#这里只是一个对sys的引用，只能reload才能进行重新加载
stdi,stdo,stde=sys.stdin,sys.stdout,sys.stderr 
reload(sys) 
#通过import引用进来时,setdefaultencoding函数在被系统调用后被删除了，所以必须reload一次
sys.stdin,sys.stdout,sys.stderr=stdi,stdo,stde 
sys.setdefaultencoding('utf-8')

# 根目录
base_dir = "//192.168.1.6/viewshed"
# 点位目录
dir_point = "//192.168.1.6/viewshed/points"
# DEM目录
dir_dem = "//192.168.1.6/viewshed/dem"
# 缓冲区图层透明度
bufferTransparency = 30
# 可视域图层透明度
viewshedTransparency = 50

# 项目文件地址logging.warn
path_project = base_dir + "/viewshed_10_3.mxd"
# temp目录
dir_temp = base_dir + "/temp"
# temp目录:buffer
dir_temp_buffer = base_dir + "/temp/buffer.shp"
# temp目录:viewshed
dir_temp_viewshed = base_dir + "/temp/viewshed.shp"
# results目录
dir_results = base_dir + "/results"


# Step0: 预处理
# 查找必要数据
# 输入要素类
arcpy.env.workspace = dir_point
try:
    point = arcpy.ListFeatureClasses("","Point")[0]
    print("已找到观测点:" + str(point))
except Exception as result: 
    raise Exception("未找到观测点: %s" % result)
# 输入表面栅格
arcpy.env.workspace = dir_dem
try:
    dem = arcpy.ListRasters("*", "TIF")[0]
    print("已找到DEM:" + str(dem))
except Exception as result: 
    raise Exception("未找到DEM: %s" % result)

# 清空temp/result目录，避免异常
print("清空temp目录:开始")
arcpy.env.workspace = dir_temp
for tempRaster in arcpy.ListRasters("*", "TIF"):
    try:
        arcpy.Delete_management(tempRaster)
        print("删除"+str(tempRaster))
    except Exception as result: 
        raise Exception("删除%s失败: %s" % (str(tempRaster),result))
for tempFeature in arcpy.ListFeatureClasses():
    try:
        arcpy.Delete_management(tempFeature)
        print("删除"+str(tempFeature))
    except Exception as result: 
        raise Exception("删除%s失败: %s" % (str(tempFeature),result))
print("清空temp目录:完成")

print("清空result目录:开始")
arcpy.env.workspace = dir_results
for tempRaster in arcpy.ListRasters("*", "All"):
    arcpy.management.Delete(tempRaster)
    print("删除"+str(tempRaster))
for tempFeature in arcpy.ListFeatureClasses():
    arcpy.management.Delete(tempFeature)
    print("删除"+str(tempFeature))
print("清空result目录:完成")

# Step1: 生成缓冲区
# 与要缓冲的输入要素之间的距离
distanceField = "distance"
# 指定将在输入要素的哪一侧进行缓冲
sideType = "FULL"
# 指定线输入要素末端的缓冲区形状
endType = "ROUND"
# 指定移除缓冲区重叠要执行的融合类型
dissolveType = "NONE"
# 融合输出缓冲区所依据的输入要素的字段列表
dissolveField = ""
# 指定是使用平面方法还是测地线方法来创建缓冲区
method = "PLANAR"

print("创建缓冲区:开始")
buffer = arcpy.Buffer_analysis(os.path.join(dir_point, point), dir_temp_buffer, distanceField, sideType, endType,dissolveType, dissolveField,method)
print("创建缓冲区:完成")


# Step2: 计算可视域
# 地表以上 (AGL) 输出栅格
outAGL = ""
# 可见性分析类型
analysisType = "FREQUENCY"
# 表面高程值中不确定项（均方根错误，或称 RMSE）的数量
verticalError = ""
# 用于识别对于每个观察点都可见的区域的输出表
outAnalysisRelationTable = ""
# 空气中可见光的折射系数
refractCoeff = "0.13"
# 该值将以表面单位指示要添加到各目标像元 z 值的垂直距离
surfaceOffset = ""
# 此值用于定义观察点或折点的表面高程
observerElevation = ""
# 该值将以表面单位指示要添加到观察点高程的垂直距离
observerOffset = "tower_h"
# 此值用于定义确定可见性的起始（最小）距离
innerRadius = ""
# 内半径参数的距离类型
innerIs3D = "False"
# 此值用于定义确定可见性的最大距离
outerRadius = "distance"
# 外半径参数的距离类型
outerIs3D = "True"
# 该值定义水平扫描范围的起始角度。此值以度为单位，介于 0 至 360.0 度之间，其中 0 指向北。默认值为 0
horizStartAngle = ""
# 该值定义水平扫描范围的终止角度。此值以度为单位，介于 0 至 360.0 度之间，其中 0 指向北。默认值为 360.0
horizEndAngle = ""
# 该值定义扫描的（位于水平面上）垂直角上限。该值应以度为单位，介于 0 到 90.0 之间，可为整数或浮点。默认值为 90.0
vertUpperAngle = ""
# 该值定义扫描的（位于水平面下）垂直角下限。该值应以度为单位，介于 -90.0 到 0 之间，可为整数或浮点。默认值为 -90.0
vertLowerAngle = ""
# 选择用于计算可见性的方法
analysisMethod = "ALL_SIGHTLINES"

print("计算可视域:开始")
viewshedResult = Viewshed2(os.path.join(dir_dem, dem), os.path.join(dir_point, point), outAGL, analysisType,verticalError, outAnalysisRelationTable, refractCoeff,surfaceOffset, observerElevation, observerOffset,innerRadius, innerIs3D, outerRadius, outerIs3D,horizStartAngle, horizEndAngle, vertUpperAngle,vertLowerAngle, analysisMethod)
print("计算可视域:转shapefile")
viewshed = arcpy.RasterToPolygon_conversion(viewshedResult, dir_temp_viewshed, "SIMPLIFY", "VALUE")
print("计算可视域:完成")


# Step3: 绘制专题图
# 查找必要图层并替换数据源
mxd = arcpy.mapping.MapDocument(path_project)
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    if(layer.isServiceLayer):
        print("ServiceLayer:"+layer.name)
    if(layer.isRasterLayer):
        print("RasterLayer:"+layer.name)
        layer.visible = False
    if(layer.isFeatureLayer):
        print("FeatureLayer:"+layer.name)
        if(layer.supports("DATASOURCE") and layer.supports("workspacePath")):
            print("图层dataSource:"+layer.dataSource)
            if(layer.dataSource.endswith(str(point))):
                lyr_point = layer
            elif(layer.dataSource.endswith(str(buffer))):
                lyr_buffer = layer
            elif(layer.dataSource.endswith(str(viewshed))):
                lyr_viewshed = layer
            else:
                layer.visible = False
if(lyr_point):
    print("找到观测点图层:"+str(lyr_point.dataSource))
    lyr_point.replaceDataSource(lyr_point.workspacePath,"SHAPEFILE_WORKSPACE")
else:
    raise Exception("未找到观测点图层")
if(lyr_buffer):
    print("找到buffer图层:"+str(lyr_buffer.dataSource))
    lyr_buffer.replaceDataSource(lyr_buffer.workspacePath,"SHAPEFILE_WORKSPACE")
    print("设置buffer图层透明度:",bufferTransparency)
    lyr_buffer.transparency = bufferTransparency
else:
    raise Exception("未找到buffer图层")
if(lyr_viewshed):
    print("找到Viewshed图层:"+str(lyr_viewshed.dataSource))
    lyr_viewshed.replaceDataSource(lyr_viewshed.workspacePath,"SHAPEFILE_WORKSPACE")
    print("设置Viewshed图层透明度:",viewshedTransparency)
    lyr_viewshed.transparency = viewshedTransparency
else:
    raise Exception("未找到Viewshed图层")
try:
    mxd.save()
except Exception as result: 
    raise Exception("项目打开期间无法保存，请新建空白项目运行此脚本: %s" % result)

# 开始导出
print("导出专题图:开始")
mxd = arcpy.mapping.MapDocument(path_project)
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    if(layer.supports("DATASOURCE") and layer.supports("workspacePath")):
        if(layer.dataSource.endswith(str(point))):
            lyr_point = layer
        if(layer.dataSource.endswith(str(buffer))):
            lyr_buffer = layer
        if(layer.dataSource.endswith(str(viewshed))):
            lyr_viewshed = layer

for pageNum in range(1, mxd.dataDrivenPages.pageCount+1):
    mxd.dataDrivenPages.currentPageID = pageNum
    mapName = mxd.dataDrivenPages.pageRow.getValue(mxd.dataDrivenPages.pageNameField.name)
    if (mxd.dataDrivenPages.pageNameField.type =="Double"):
        sql = '"%s" = %s' %(mxd.dataDrivenPages.pageNameField.name,str(mapName))
    elif (mxd.dataDrivenPages.pageNameField.type =="Integer"):
        sql = '"%s" = %s' %(mxd.dataDrivenPages.pageNameField.name,str(mapName))
    elif (mxd.dataDrivenPages.pageNameField.type =="Single"):
        sql = '"%s" = %s' %(mxd.dataDrivenPages.pageNameField.name,str(mapName))
    elif (mxd.dataDrivenPages.pageNameField.type =="SmallInteger"):
        sql = '"%s" = %s' %(mxd.dataDrivenPages.pageNameField.name,str(mapName))
    else :
        sql = '"%s" = \'%s\'' %(mxd.dataDrivenPages.pageNameField.name,str(mapName))

    lyr_point.definitionQuery = lyr_buffer.definitionQuery = sql
    print("导出专题图:%s,查询条件为：%s"% (str(mapName),sql))
    arcpy.mapping.ExportToJPEG(mxd, dir_results +"/"+ str(mapName) +".jpg")

print("导出专题图:完成")

# 重置筛选条件
lyr_point.definitionQuery = ""
lyr_buffer.definitionQuery = ""
mxd.save()
del point,buffer,dem,viewshed,mxd