遥感数据处理系列

一些项目及科研中遇到的小需求，一方面记录自己的学习历程，另一方面帮助大家学习。本系列文章的开发环境为：ArcGIS 10.2.2 + Python 2.7、ENVI 5.3 + IDL 8.5

文章目录

  • 遥感数据处理系列
  • 前言
  • 一、计算栅格数据平均值
  • 1. 原理简介
  • 2. 代码
  • 二、栅格异常值处理
  • 总结
  • 后记

前言

ArcPy这个包也太重要了吧！如果没有IDL+Python+Matlab，我的实验又该如何展开？如果没有ArcPy，那可能就要用GDAL硬撕代码了。本文介绍如何处理单个栅格数据的平均值（最大值、最小值同理）。

一、计算栅格数据平均值

1. 原理简介

计算栅格数据平均值主要使用ArcPy的GetRasterProperties_management函数。

函数使用：

GetRasterProperties_management (in_raster, {property_type}, {band_index})

常用参数简介：

in_raster：要检索属性的栅格文件
	property_type：
	        MINIMUM —输入栅格中所有像元的最小值。
	        MAXIMUM —输入栅格中所有像元的最大值。
	        MEAN —输入栅格中所有像元的平均值。（本文主要使用该标签）

2. 代码

输入： 一个含有若干栅格数据的文件夹（本例为“.tif”格式）

输出： 一个文本文件，在当前目录下（与代码路径相同）

代码实例：

# -*- coding: UTF-8 -*-
import os
import glob
import arcpy
from arcpy.sa import *

'''
功能：
    计算输入文件夹inws内，所有栅格数据的平均值
'''

arcpy.CheckOutExtension("ImageAnalyst")  # 检查许可
arcpy.CheckOutExtension("spatial")

# 输入路径  应该注意，中文路径，会导致读不出文件；路径尽量不要有空格，写文件时会报错
inws = r"D:\LE-daily\2013"

OutputFile = open('2013-average.csv', 'w')

# 利用glob包，将inws下的所有tif文件读存放到rasters中
rasters = glob.glob(os.path.join(inws, "*.tif"))

# 循环rasters中的所有影像，进行“求平均值”操作
for ras in rasters:
    meanValueInfo = arcpy.GetRasterProperties_management(ras, 'MEAN')
    meanValue = meanValueInfo.getOutput(0)
    print os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n'
    OutputFile.write(os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n')

OutputFile.close()
print("All project is OK！")

以上代码适用于不含异常值的场景，那么，如果有异常值呢？那么多的背景值，或是计算时为标记异常情况设置的其它flag（例如：NDVI计算时，标记异常情况为-2…）

二、栅格异常值处理

栅格存在异常值如何处理？这里借助SetNull函数进行：读取栅格数据后，先进行一步去除异常值，接着再进行一步平均值获取，Perfect！

代码如下：

# -*- coding: UTF-8 -*-
import os
import glob
import arcpy
from arcpy.sa import *

'''
功能：
    计算输入文件夹inws内，所有栅格数据的平均值
'''

arcpy.CheckOutExtension("ImageAnalyst")  # 检查许可
arcpy.CheckOutExtension("spatial")

# 输入路径  应该注意，中文路径，会导致读不出文件；路径尽量不要有空格，写文件时会报错
inws = r"D:\LE-daily\2013"

OutputFile = open('2013-average.csv', 'w')

# 利用glob包，将inws下的所有tif文件读存放到rasters中
rasters = glob.glob(os.path.join(inws, "*.tif"))

whereClause = "VALUE = 0"  # 去除异常值

# 循环rasters中的所有影像，进行“求平均值”操作
for ras in rasters:
    outSetNull = SetNull(ras, ras, whereClause)  # 去除异常值

    meanValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'MEAN')
    meanValue = meanValueInfo.getOutput(0)
    print os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n'

    OutputFile.write(os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n')

OutputFile.close()
print("All project is OK！")

总结

ArcPy牛皮！毕业万岁！圣诞节快乐！感谢晶晶晶的圣诞帽！
计算最大值、最小值，参考上例，修改MEAN标签为MAXIMUM、MINIMUM即可。