Python 批量重采样、掩膜、坡度提取

DataCharm 2021-02-22 13:02:32
Python 提取 批量 采样 坡度


今日分享:

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

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 arcpy
in_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,glob
from arcpy import env
arcpy.CheckOutExtension("Spatial")
filepath=r"C:\\Users\\Admin\\Desktop\\GISPractice\\ resample"
env.workspace = filepath
os.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 arcpy
from arcpy import env
env.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 arcpy
from arcpy import env
from 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

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

本文分享自微信公众号 - DataCharm(shujumeili)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间: 2020-11-15

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

版权声明
本文为[DataCharm]所创,转载请带上原文链接,感谢
https://cloud.tencent.com/developer/article/1790264

  1. Python 3 entry, see this is enough
  2. 华为大佬打造的400集Python视频学起来,学完万物皆可爬
  3. 400 episodes of Python video created by Huawei boss
  4. django之csrf_exempt解决跨域请求的问题
  5. CSRF of Django_ Exempt solves the problem of cross domain requests
  6. 1.7 万 Star!一个简单实用的 Python 进度条库
  7. 17000 stars! A simple and practical Python progress bar library
  8. Python爬虫:设置Cookie解决网站拦截并爬取蚂蚁短租
  9. Python crawler: setting cookie to solve website interception and crawling ant short rent
  10. Python-Net编程
  11. Python net programming
  12. 学习Python数学英语基础重要吗?Python教程!
  13. Is it important to learn the basics of math and English in Python!
  14. Python数据分析常用库有哪些?Python学习!
  15. What are the common libraries for Python data analysis? Learn Python!
  16. win 创建python虚拟环境
  17. Creating Python virtual environment with win
  18. In order to automatically collect B station barrage, I developed a tool in Python
  19. 用Python编程语言来实现阿姆斯特朗数的检查
  20. Using python programming language to check Armstrong number
  21. Python中的解决中文字符编码的问题
  22. Solving the problem of Chinese character coding in Python
  23. Translation: practical Python Programming 02_ 01_ Datatypes
  24. Installation and use of Python and tensorflow in win10 environment (Python version 3.6, tensorflow version 1.6)
  25. Python series 46
  26. Linux安装Python3
  27. 【python接口自动化】- 正则用例参数化
  28. Python RestFul Api 设计
  29. filecmp --- 文件及目录的比较│Python标准库
  30. Installing python3 on Linux
  31. [Python] Matplotlib 圖表的繪製和美化技巧
  32. (資料科學學習手札108)Python+Dash快速web應用開發——靜態部件篇(上)
  33. 翻譯:《實用的Python程式設計》02_01_Datatypes
  34. 【python接口自动化】- 正则用例参数化
  35. 翻译:《实用的Python编程》02_02_Containers
  36. 两年Java,去字节跳动写Python和Go
  37. [Python interface automation] - regular use case parameterization
  38. Python restful API design
  39. 翻译:《实用的Python编程》02_02_Containers
  40. 两年Java,去字节跳动写Python和Go
  41. 翻译:《实用的Python编程》02_02_Containers
  42. Python基于粒子群优化的投资组合优化研究
  43. ubuntu部署django项目
  44. 兩年Java,去位元組跳動寫Python和Go
  45. 翻譯:《實用的Python程式設計》02_02_Containers
  46. 这样学习Python,爷爷都学会了!超简单Python入门
  47. [Python] 基于 jieba 的中文分词总结
  48. 【python】递归听了N次也没印象,读完这篇你就懂了
  49. [Python] 基于 jieba 的中文分词总结
  50. 人理解迭代,神则体会递归,从电影艺术到Python代码实现神的逆向思维模式
  51. [Python] 基於 jieba 的中文分詞總結
  52. Python属于后端开发还是前端开发?Python入门!
  53. 【python】递归听了N次也没印象,读完这篇你就懂了
  54. 一天快速入门python
  55. 学习Python对年龄有没有要求?30岁可以吗?
  56. 清华教授!12小时整理的最全Python教程(文末无偿分享)
  57. Filecmp -- comparison of files and directories
  58. Drawing and beautifying skills of [Python] Matplotlib chart
  59. Python + dash rapid web application development static components
  60. Translation: practical Python Programming 02_ 01_ Datatypes