python提取斜坡结构
介绍
开始是在帮师妹处理某个试验流程中发现需要进行斜坡结构的提取,后面百度找了教程一步一步的做,发现挺麻烦的,所有写了一段代码,所以里面文件夹名字可能emmmm,不重要,这些步骤主要还是python二开,比较简单,如果有什么写的不好的请大家多多包涵。这是我的第一篇博客,希望能有个好的开始把。
代码
import arcpy
from arcpy import env
from arcpy.sa import *
import math
def Xiepojiegou():
env.workspace = "E:\\gisforpyexcel\\mandata\\arcpyresult"
# 这是输入的坡向和倾向(倾角和倾向都是用产状点插值生成的)
inRaster1 = arcpy.Raster("aspect0516.tif")
inRaster2 = arcpy.Raster("倾向裁剪.tif")
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
# 进行栅格计算
outPlus1 = inRaster1 - inRaster2
outPlus2= Abs(outPlus1)
# Save the output
outPlus2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\pojianqinxiang.tif")
# 切换工作空间
env.workspace = "E:\\gisforpyexcel\\mandata\\斜坡结构"
# 这个是重分类的设置,首先不考虑倾角大小和坡向-倾向的绝对值大小,得出顺向坡、顺斜坡、横向坡、逆斜坡、逆向坡
a=[[0,30,0],[30,60,4],[60,120,5],[120,150,6],[150,210,7],[210,240,6],[240,300,5],[300,330,4],[330,360,0]]
outReclass1 = Reclassify("pojianqinxiang.tif", "Value", RemapRange(a))
outReclass1.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\reclass04567.tif")
env.workspace = "E:\\gisforpyexcel\\mandata\\arcpyresult"
# 对坡度和倾角大小进行判断,得出顺向飘倾斜坡、顺向层面坡、顺向伏倾坡
inRaster3 = arcpy.Raster("slope0516.tif")
inRaster4 = arcpy.Raster("倾角裁剪1.tif")
outPlus2 = inRaster3 - inRaster4
# Save the output
outPlus2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif")
# 对倾角以5°为间隔重分类,后面要用
outReclassqinjiao = Reclassify("倾角裁剪1.tif", "Value",
RemapRange([[0, 5, 0], [5, 90, 1]]))
outReclassqinjiao.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\Remapqinjiao05.tif")
print 3
env.workspace = "E:\\gisforpyexcel\\mandata\\斜坡结构"
# 读取输出的最大值和最小值,来判别关系
mmax = arcpy.GetRasterProperties_management("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif", "MAXIMUM")
mmin= arcpy.GetRasterProperties_management("E:\\gisforpyexcel\\mandata\\斜坡结构\\podujianqinjiao.tif","MINIMUM")
# 进行格式转换
max = int(math.ceil(float(mmax[0])))
min = int(math.floor(float(mmin[0])))
print min
print max
print 4
# 这个重分类我感觉我写的有问题,但是我不知道怎么同时运用RemapRange和RemapValue,实际应该是[min,0,1]、[0,0,2]、[0,max,3],但是会报错,因为range是[0,0],错误的范围
outReclass2 = Reclassify("podujianqinjiao.tif", "Value", RemapRange([[min,-0.000001,1],[-0.000001,0.000001,2],[0.000001,max,3]]))
outReclass2.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\reclass123.tif")
print 5
inRaster5 = arcpy.Raster("reclass04567.tif")
inRaster6 = arcpy.Raster("reclass123.tif")
outPlus3 = inRaster5 + inRaster6
# Save the output
outPlus3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\shunxiangpo12345678910.tif")
# 进行重分类的目的是为了之后的栅格合并不会出现值的交叉错误,利于区分
a = [[0, 1, 1], [1, 2, 2], [2, 3, 3],[3,10,0]]
outReclass3 = Reclassify("shunxiangpo12345678910.tif", "Value", RemapRange(a))
outReclass3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\shunxiangpo.tif")
inRaster7 = arcpy.Raster("shunxiangpo.tif")
inRaster8 = arcpy.Raster("reclass04567.tif")
outPlus4 = inRaster7 + inRaster8
# Save the output,这里得到的就是不考虑倾角关系下的其余7种,但是包含近水平状坡,后面排除
outPlus4.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\all1234567.tif")
inRaster9 = arcpy.Raster("Remapqinjiao05.tif")
inRaster10 = arcpy.Raster("all1234567.tif")
outPlus4 = inRaster9 * inRaster10
# Save the output
outPlus4.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\zzall01234567.tif")
a = [[-0.5, 0.5, 0], [0.5, 1.5, 1], [1.5, 2.5, 2],[2.5,3.5,3],[3.5,4.5,4],[4.5,5.5,5],[5.5,6.5,6],[6.5,7.5,7]]
outReclass3 = Reclassify("zzall01234567.tif", "Value", RemapRange(a))
outReclass3.save("E:\\gisforpyexcel\\mandata\\斜坡结构\\all01234567.tif")
arcpy.AddField_management("all01234567.tif", 'Class', "TEXT")
# 为了避免忘记,属性表加入字段名字
try:
cur = arcpy.UpdateCursor("all01234567.tif")
for row in cur:
# print row
if(row.Value==0):
row.Class = "近水平层状坡"
if(row.Value==1):
row.Class = "顺向飘倾坡"
if (row.Value == 2):
row.Class = "顺向层面坡"
if (row.Value == 3):
row.Class = "顺向伏倾坡"
if (row.Value == 4):
row.Class = "顺斜坡"
if (row.Value == 5):
row.Class = "横向坡"
if (row.Value == 6):
row.Class = "逆斜坡"
if (row.Value == 7):
row.Class = "逆向坡"
cur.updateRow(row)
del cur, row
except:
print arcpy.GetMessages()
注意事项
在进行处理之前要确保各要素的坐标一致,否则就会出现栅格点缺失的情况。
在各位进行代码复制的时候,修改文件夹的名字就可以了。我的基础数据是放在第一个文件夹里面。