利用ArcGIS Python批量拼接裁剪遥感影像(arcpy batch processing)

转载自http://blog.csdn.net/gisboygogogo/article/details/75195760


本篇文章将说明如何利用ArcGIS 10.1自带的Python IDLE进行遥感影像的批量拼接与裁剪。

1.运行环境:ArcGIS10.1 (安装传送门)、Python IDLE

2.数据来源:地理空间数据云 GDEMV2 30M分辨率数字高程数据

3.解决问题:制作山西省的DEM影像

如下图所示,以30M分辨率数字高程数据为例,影像皆是固定范围的经纬度保存在其服务器上,外在表现以小幅正方形影像。如果手动进行拼接,工作量会非常大且容易出错。

我们的目标:批处理,写一次代码,处理多幅影像


1.查找目标范围的经纬度信息,本文以山西为例,经纬度范围在N33-N41, E109-E114直接,所以先将这个范围内的影像下载后,解压后放在一个文件夹下。


2.拼接影像

开始菜单-ArcGIS-IDLE(Python GUI) 打开IDLE。File-New Window, 文件以 .py结尾,如MosaicToNewRasters.py

[python]  view plain  copy
  1. import arcpy  
  2. import os  
  3.   
  4. #指定工作目录,即存放影像的目录  
  5. arcpy.env.workspace = r"E:\Huangkun\arcpyData\shanxi\N34_N35"  
  6.   
  7. #指定该工作空间下的一副影像为基础影像,为后面的参数提取做准备  
  8. base = "ASTGTM2_N34E109_dem.tif"  
  9.   
  10. #以下一段代码是为执行拼接做参数准备  
  11. out_coor_system = arcpy.Describe(base).spatialReference #获取坐标系统  
  12. dataType = arcpy.Describe(base).DataType   
  13. piexl_type = arcpy.Describe(base).pixelType   
  14. cellwidth = arcpy.Describe(base).meanCellWidth #获取栅格单元的的宽度  
  15. bandcount = arcpy.Describe(base).bandCount #获取bandCount  
  16.   
  17. #打印一些信息  
  18. print out_coor_system.name  
  19. print dataType  
  20. print piexl_type  
  21. print cellwidth  
  22. print bandcount  
  23.   
  24. arcpy.CheckOutExtension("Spatial")  
  25.   
  26. #提取待拼接影像的文件名,且中间以;隔开,例如:a.tif;b.tif;c.tif  
  27. rasters = []  
  28. for ras in arcpy.ListRasters("*dem.tif"):    #for循环,对wrokspace下的所有以dem.tif结尾的影像进行过滤  
  29.     rasters.append(ras)  
  30.   
  31. ras_list = ";".join(rasters)                 #字符串拼接  
  32.   
  33. #打印出来,看看什么结果吧  
  34. print ras_list  
  35.   
  36.   
  37. #指定输出文件夹  
  38. outFolder = r"E:\Huangkun\arcpyData\shanxi\N34_N35\n34_n35_out"  
  39.   
  40. #执行拼接操作  
  41. arcpy.MosaicToNewRaster_management(ras_list, outFolder, "shanxi_n34_n35_dem_data.tif", out_coor_system, "16_BIT_SIGNED", cellwidth, bandcount, "LAST""FIRST")  
RUN-->Run Module F5 执行代码

这里用到了arcpy包下面的 MosaicToNewRaster_management()函数

得到如下结果


3.以山西省shp边界裁剪影像


再新建一个Window,文件命为BatchExtractByMask.py

[python]  view plain  copy
  1. import arcpy  
  2. import glob  
  3. import os  
  4.   
  5. arcpy.CheckOutExtension('Spatial')  
  6.   
  7. #指定先前拼接后的遥感影像所在目录  
  8. inws = r"E:\Huangkun\arcpyData\shanxi\shanxi_regular_dem\shanxi_province_dem"  
  9.   
  10. #指定裁剪后的影响存放目录  
  11. outws = r"E:\Huangkun\arcpyData\shanxi\shanxi_regular_dem\shanxi"  
  12.   
  13. #指定shp范围边界文件,即目标区域的边界  
  14. mask = r"E:\Huangkun\arcpyData\shanxi\shanxi_shp\shanxi.shp"  
  15.   
  16. #利用glob包,将inws下的所有tif文件读存放到rasters中  
  17. rasters = glob.glob(os.path.join(inws, "*.tif"))  
  18.   
  19. #循环rasters中的所有影像,进行按掩模提取操作  
  20. for ras in rasters:  
  21.     outname = os.path.join(outws, os.path.basename(ras).split(".")[0] + "_clp.tif")  #指定输出文件的命名方式(以被裁剪文件名+_clip.tif命名)  
  22.     out_extract = arcpy.sa.ExtractByMask(ras, mask)  #执行按掩模提取操作  
  23.     out_extract.save(outname)  #保存数据  

执行脚本代码,裁剪后的影像被存放在outw所指定的文件夹


4.通过ArcMap加载影像,最终结果:


这就实现了批量拼接与裁剪影像的工作,其实拼接与裁剪无所谓先后顺序。

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 程序猿惹谁了 设计师:白松林 返回首页