[ArcPy] 遥感影像去黑边-第六届全国大学生GIS技能大赛试题

本文主要是介绍[ArcPy] 遥感影像去黑边-第六届全国大学生GIS技能大赛试题，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

完整试题：http://www.higis.cn/753

问题

去除遥感影像黑边

数据：

分析

【方法一】将边界进行矢量化，然后用掩膜提取
缺点：遥感影像边缘是锯齿形，矢量化边界精度不高

【方法二】边界像元值为0，将0值重分类为NoData
缺点：将内部像元值为0的像素也归为了NoData，会出现空洞
发现ArcGIS中有类似功能的工具：ArcToolbox–>Data Management Tools–>Raster–>Raster Dataset–>CopyRaster

【方法三】法一法二联合做就可以解决他们两个的缺点
法一的结果为1.tif，法二的结果为2.tif
栅格计算器：Con(IsNull(“2”) & ( ~IsNull(“1”)), “1”, “2” )
解释：如果2为NoData的地方1有值，那么就赋予1的值，其他地方赋予2

【方法四】使用ArcPy仅将边缘像元为0的地方赋予NoData，里面的地方不处理
缺点：处理速度慢，遥感图像n行m列，时间复杂度O(n*m)

【方法五】使用边界跟踪算法得出边界，存储成shp格式，再进行掩膜

ArcPy暴力剔除边缘的像元值

代码

# -*- coding:utf-8 -*-
# Author: PasserQi
# Time: 2017/12/18
# Func: 去边处理
# Desc: 从左边向右，从右边向左开始搜索，将边界像素替换成NoData，只会替换到边界的像素，不会替换到内部像元
# Attention: 若要打包成ArcGIS自定义工具，此文件下不能有中文
import arcpy,sys
from arcpy.sa import *arcpy.CheckOutExtension("sptial") #权限检查
# 得到初始参数
# rasterPath = sys.argv[1]
# outPath = sys.argv[2]
# BORDER_VALUE = int(sys.argv[3])
#
rasterPath = r'E:\user\Desktop\DeleteDarkSide\Data\Kejicheng0.TIF' #输入栅格
outPath = r'E:\user\Desktop\DeleteDarkSide\Data\q' #输出栅格
BORDER_VALUE = 0 #边界值if __name__ == '__main__':r = Raster(rasterPath) #打开栅格noDataValue = r.noDataValue #获取栅格中的NoData值array = arcpy.RasterToNumPyArray(r) #转成Numpy方便对每个像元进行处理bandNum,rowNum,colNum = array.shape #波段、行数、列数# 左边for i in range(0,rowNum):for j in range(0,colNum):if array[0][i][j]==BORDER_VALUE: #找到边界for w in range(0,bandNum):array[w][i][j] = noDataValue #赋值print "替换:%d %d" % (i,j)continueelif array[0][i][j]==noDataValue: #是无值continueelse: #是普通像素break #退出该行# 右边for i in range(0,rowNum):for z in range(0,colNum):j = colNum - 1 - zif array[0][i][j]==BORDER_VALUE: #找到边界for w in range(0,bandNum):array[w][i][j] = noDataValue  # 赋值print "替换:%d %d" % (i, j)continueelif array[0][i][j]==noDataValue: #是无值continueelse: #是普通像素break #退出该行#保存栅格lowerLeft = arcpy.Point(r.extent.XMin, r.extent.YMin) #左下角点坐标cellWidth = r.meanCellWidth #栅格宽度cellHeight = r.meanCellHeightnewRaster = arcpy.NumPyArrayToRaster(array,lowerLeft,cellWidth,cellHeight,noDataValue) #转换成栅格newRaster.save(outPath) #保存

将Py脚本制定成自定义工具

执行: 去黑边 "E:\user\Desktop\[ArcPy] 去除遥感影像黑边\Data\Kejicheng0.TIF" "E:\user\Desktop\[ArcPy] 去除遥感影像黑边\result\tmp.tif" 0
开始时间: Thu Dec 21 20:15:58 2017
正在运行脚本 去黑边...
Completed script 去黑边...
成功 在 Thu Dec 21 20:19:02 2017 (经历的时间: 3 分 4 秒)

注：可以看到这个方法时间要很久

结果

找出边界后掩膜

# -*- coding:utf-8 -*-
# Author: PasserQi
# Time: 2017/12/19
# Func: 提取出遥感图像的边界
# Desc: 使用边界跟踪算法，此实例使用八连通邻域，输入的边界为
import arcpy,sys,threading
from arcpy.sa import *rasterPath = r'E:\user\Desktop\[ArcPy] 去除遥感影像黑边\Data\Kejicheng0.TIF' #输入栅格
outPath = r'E:\user\Desktop\DeleteDarkSide\Data\q' #输出栅格
BORDER_VALUE = 0 #边界值directions = [ [0,1],[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1],[1,0],[1,1]] #方向数组
def following(dir,boundary): #内边界start_i,start_j = boundary[0] #找出第一个x,y = boundary[-1] #找出最后一个# print boundary[0] + boundary[-1]if len(boundary)!=1 and start_i==x and start_j==y: #闭合了returnx = x + directions[dir][0] #移动y = y + directions[dir][1] #移动flag = False #标识是否找到下一个边界if x<0 or y<0 or x>=rowNum or y>=colNum: #出边界flag = False #没有找到，继续找elif array[0][x][y]!=BORDER_VALUE and array[0][x][y]!=noDataValue: #有像素flag = Trueboundary.append([x,y] )if x==start_i and y==start_j:return#没有找到下一个继续搜索if flag==True: #找到了下一个边界像素if dir%2==0:dir = (dir+7)%8else:dir = (dir+6)%8else: #未找到dir = (dir+1)%7print len(boundary)following(dir,boundary)if __name__ == '__main__':sys.setrecursionlimit(1000000)  # 例如这里设置为一百万r = Raster(rasterPath) #打开栅格noDataValue = r.noDataValue #获取栅格中的NoData值array = arcpy.RasterToNumPyArray(r) #转成Numpy方便对每个像元进行处理bandNum,rowNum,colNum = array.shape #波段、行数、列数# 找出图像最左上方的像素flag = Falsefor i in range(0,rowNum):for j in range(0,colNum):if array[0][i][j]!=BORDER_VALUE and array[0][i][j]!=noDataValue:start_i = istart_j = jflag = Truebreakif flag:breakboundary = [] #边界数组boundary.append([start_i,start_j] )print boundary# 开始搜索dir = 7 #起始搜索方向threading.stack_size(200000000) #堆大小设置thread = threading.Thread(target=following(dir,boundary))thread.start()

这篇关于[ArcPy] 遥感影像去黑边-第六届全国大学生GIS技能大赛试题的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---