[Py] 利用 ArcPy 進行批次轉換 (DEM -> SHD)

如果需要準備大量不同方向、仰角的 Hillshade 的話,可以參考以下範例:

本文可以見到:

  • 使用 ArcPy 時需要先確認相關模組與授權是否合法。
  • 使用 arcpy.HillShade_3d() 生成 Hillshade 檔。
  • 使用 arcpy.Slope_3d() 生成 Slope 檔。

程式碼如下:

# -*- coding: utf-8 -*-
"""
HP(HillshadeProducer).py
Created on Mon Oct 31 16:09:45 2016

@author: WayneSun
"""
import time
Ts=time.time()
from sys import exit,argv
from os import getcwd,mkdir
import os.path as oph
print("This is a Python 2.7 script.\nRequire 3D extension in ArcGIS")
print("Usage : in command line tool:\npython HP.py DEM.ext OutputDir(optional)")
try:
    print("Initializing arcpy ...")
    import arcpy
    from arcpy import env
    print("Done! Time : {}s".format(time.time()-Ts))
except ImportError:
    print("Make sure you have installed \
    ArcGIS and use the python includes the package.")
    exit(1)

#Check licence
#http://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-functions/checkoutextension.htm    
print("3D analyst toolbox : {}".format(arcpy.CheckExtension("3D")))
print("Licence of 3D analyst toolbox : {}".format(arcpy.CheckOutExtension("3D")))   

# Set environment settings
env.overwriteOutput = True


inRaster = oph.abspath(argv[1])
print("inRaster:{}".format(inRaster))
inBasename=oph.basename(inRaster).split(".")[0]
if not oph.isfile(inRaster) :
    print("File : {0}\ndoes not exist. Try another file.")
    exit(0)
outFolder = oph.join(getcwd(),"OutputRaster")
if not oph.exists(outFolder):
    mkdir(outFolder)


env.workspace = getcwd()

# Set local variables

#outRaster = "C:/output/outhillshd02"
#azimuth = 180
#altitude = 75
modelShadows = "SHADOWS"
zFactor = 1
for azi in range(0,360,10):
    for alt in range(10,91,10):
        Tss=time.time()        
        outRaster = oph.join(outFolder,"{0}-az{1}-al{2}.tif".format(inBasename,azi,alt))
        print("Now processing : {}".format(outRaster))
        # Execute HillShade
        arcpy.HillShade_3d(inRaster, outRaster, azi, alt, modelShadows, zFactor)
        print("Time: {}s".format(time.time()-Tss))
        
#Slope
Tss=time.time()
outRaster=oph.join(outFolder,"{}-slp".format(inBasename))
#Slope_3d (in_raster, out_raster, {output_measurement}, {z_factor})
arcpy.Slope_3d(inRaster, outRaster, "DEGREE",1.)
print("Slope processed! Time : {}s".format(time.time()-Tss))

print("\n\nAll done! Time : {}s".format(time.time()-Ts))

 

請多多指教!

這個網站採用 Akismet 服務減少垃圾留言。進一步瞭解 Akismet 如何處理網站訪客的留言資料