查看原文
其他

【GIS实验课】05 坡度提取的不确定性分析(分享批量重采样、批量掩膜、坡度批量提取代码)

不许人间见白头 数读城事 2022-04-25
今日分享:

后台回复“批量”可以获取批量重采样、批量掩膜、批量坡度提取和批量分区统计的代码,不过你们懂得。



01主要内容


本次实验下载的是GDEMV2 30M分辨率数字高程数据,利用Python提取不同分辨率的DEM,基于上述不同分辨率DEM提取每种地貌类型的平均坡度,最后以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图,并对散点图进行拟合,通过回归算法求得回归方程的系数及常数项。


02具体要求


1.以30m空间分辨率的DEM数据为基础数据,重采样为40、50、60、70、80、90、100、110、120 m共10组不同分辨率的DEM。

2. 基于不同分辨率DEM提取每种地貌类型的平均坡度。

3. 以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图,并对散点图进行拟合,通过回归算法求得回归方程的系数及常数项。


03具体步骤
◐ 1. 使用ArcPy进行处理

1.1 将五景DEM数据镶嵌起来然后利用ArcPy进行批量重采样,具体代码如下所示:

import arcpyin_raster = r"C:\Users\Admin\Desktop\GIS Practice\ dem_mosaic.tif"out_raster_workspace = "C:\\Users\\Admin\\Desktop\\GISPractice\\ resample\\"for n in range(40,130,10): out_raster =out_raster_workspace + "resample_" + str(n) + ".tif" arcpy.Resample_management(in_raster, out_raster,str(n),"CUBIC")print "OK!"


1.2 将重采样得到10组不同分辨率的DEM,利用行政区的矢量边界,编写Python代码进行批量剪裁,具体代码如下所示:

import arcpy,os,globfrom arcpy import envarcpy.CheckOutExtension("Spatial")filepath=r"C:\\Users\\Admin\\Desktop\\GISPractice\\ resample"env.workspace = filepathos.chdir(filepath)Rasters = glob.glob("*.tif")outfile="C:\\Users\\Admin\\Desktop\\GISPractice\\test2\\clip"inMaskData = "C:\\Users\\Admin\\Desktop\\GISPractice\\hefei\\He.shp"for Raster in Rasters: inRaster =Raster arcpy.gp.ExtractByMask_sa(Raster, inMaskData, outfile+"\\"+Raster.split(".")[0]+'_clip.tif')print "ok"

批量剪裁结果:

 图1批量剪裁结果

1.3 将上述批量剪裁完的不同分辨率的DEM数据进行批量提取坡度,具体的Python代码如下所示:

import arcpyfrom arcpy import envenv.workspace = "C:\\Users\\Admin\\Desktop\\GISPractice\\ clip" Rasters =arcpy.ListRasters("*","tif")arcpy.CheckOutExtension("3d")for inRaster in Rasters: outMeasurement = "DEGREE" zFactor =1 outRaster = "C:\\Users\\Admin\\Desktop\\GISPractice\\ slope\\" +\ inRaster.split(".")[0] +"_Slope.tif" arcpy.Slope_3d(inRaster,outRaster, "DEGREE", zFactor)print "you are very wise!"

1.4 基于上述不同分辨率DEM提取每种地貌类型的平均坡度(批量分区统计),具体代码如下所示:

import arcpyfrom arcpy import envfrom arcpy.sa import *env.workspace ="C:\\Users\\Admin\\Desktop\\GIS Practice\\ slope"ValueRaster=arcpy.ListRasters("*","tif")arcpy.CheckOutExtension("Spatial")for inValueRasterin ValueRaster: inZoneData="C:\\Users\Admin\\Desktop\\GISPractice\\quanguodimao__WGS1984\\hefeidimao_prj.shp" zoneField = "GEOMOR_ID" outTable="C:\\Users\Admin\\Desktop\\GISPractice\\ZonalSta\\"+inValueRaster.split(".")[0]+".dbf" outZSaT=ZonalStatisticsAsTable(inZoneData,zoneField, inValueRaster,outTable,"NODATA","MEAN")print "you arevery wise!"


1.5 利用EXCEL统计不同分辨率DEM下提取的每种地貌类型的平均坡度,如下表所示:

表1|不同分辨率DEM下提取的每种地貌类型的平均坡度

以DEM分辨率为横坐标、区域平均坡度为纵坐标做不同地貌类型的散点图,并对散点图进行拟合,通过回归算法求得回归方程的系数及常数项(使用的工具是excel)。

图2不同地貌类型的散点图


表2不同地貌拟合函数表达式

由图2和表2可以看出平均坡度值与分辨率呈现很强的线性关系且是正相关关系,可以以线性方程来描述这种关系方程式为y = ax + b,其中x为DEM分辨率,y为不同DEM分辨率下的平均坡度。低海拔丘陵的斜率最大,低海拔洪积平原的斜率最小,斜率绝对值之差为0.3888。从整体上看,按照拟合曲线的斜率,可大致将上述地貌类型分为两类:(1)斜率较大类:低海拔丘陵、低海拔冲积洪积台地、低海拔冲积平原、低海拔冲积扇平原;(2)斜率略小类:低海拔小起伏山地、低海拔冲积台地、低海拔洪积平原、低海拔湖积平原、湖泊。对某市各地貌类型拟合曲线的斜率求均值作为该研究区的平均曲率,求得a = -0.32493。


◐ 2. 用Model Builder进行处理

具体模型如下所示:

 3|模型示意

在Model Builder中拖入各种数据进行建模,先加入包含不同分辨率DEM数据的文件夹clip,然后插入栅格迭代器,并设置工作空间或栅格目录为带有迭代号的文件夹clip,接着加入按掩模提取工具,将某市区域提取出来,然后加入Slope工具和分区统计工具,在分区统计工具设置中,输入要素区域数据为某市地貌矢量数据,使用地貌数据的ID字段对每种分辨率下的坡度数据进行统计,输出文件的名称为:%名称%.tif,统计类型为均值,运行该模型,则可以得到相应数据类型下的平均坡度。


Tips:

在编写ArcPy代码进行DEM数据的批量重采样的时候出现了报错,经过排查发现主要原因是因为out_raster = out_raster_workspace +"resample_" + str(n) + ".tif"这一句代码出现了错误,我们对DEM数据进行重采样,从30米到120米一共有10景DEM数据,输出的每个DEM的名称肯定是不一样的,都是根据DEM数据的分辨率来进行命名,采用的Python语句是:for n in range(40,130,10),而问题就是出现这里,这里面的n是表示数字,所以在下面的代码中需要写成str(n),因为如果不这样写的话,这个n会被认定为一个无效字符。

除此之外,在利用矢量边界对不同分辨率的DEM进行批量剪裁的时候出现了错误,在这之前我也编写ArcPy做过不少批量剪裁,不过是用不同的矢量边界去裁剪同一个栅格,遍历矢量数据的语法是:Features=arcpy.ListFiles(“*shp”),但是本次需要的是用同一个矢量边界去批量剪裁多个栅格数据,所以遍历数据的语句则改为: Rasters =glob.glob("*.tif"),在编写代码的时候我导入的库有:arcpy、os、glob,如果只是导入arcpy则程序无法执行,通过多次的调试代码终于运行成功!!!


后台回复“批量”可以获取批量重采样、批量掩膜、坡度批量提取和批量分区统计的代码,emmmmmm,不过你们懂得==


作者|不许人间见白头

排版|Moon

校阅|数读菌、不许人间见白头


那今天就到这里结束啦,欢迎留言讨论。文中的图片文字未经许可不要随便“引用”。

如果可以的话,希望能够转发分享,点个在看并且点个,给个赞赏~~也欢迎规范转载~

也希望大家和我多留言互动啊!(据说这样可以增加我的推送在你的订阅号里出现的概率)

文章推荐 【数据分享】中国生态服务功能重要性与中国生态系统敏感性【双评价学习笔记】农业:环境评价【GIS实验课】04 GIS数据的批处理方法【GIS实验课】03 城镇住房选址适宜区评价分析【ENVI实验课】02利用ENVI进行农业耕作用地变化监测【ENVI实验课】01基于像元二分模型的植被覆盖度反演【GeoDa基础】GeoDa平台下的江西省县级统计数据分析【GIS实验课】02 超市选址适宜性评价【数据整理】2019年行政区划调整的数据更新与分享【双评价学习笔记】农业:气候评价(协同克里金方法)【GIS实验课】01山区土壤侵蚀评价及其地形的关系——以皖西大别山为例【GIS基础】基本空间分析工具:缓冲区、网络分析、相交分析、密度分析

需要你的“分享”和“在看”


END>

如需全文转载文章、投稿或者合作

可添加微信

(回复超慢!!!)

(不要添加我问各种问题,我大概率不会的==)

(入群请一定要备注入群)

(添加后会在晚上非工作时间通过,请稍安勿躁)


公众号


微博

▼ 点击阅读原文,使用关键词搜索历史文章

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存