遥感数据处理系列
一些项目及科研中遇到的小需求,一方面记录自己的学习历程,另一方面帮助大家学习。本系列文章的开发环境为:ArcGIS 10.2.2 + Python 2.7、ENVI 5.3 + IDL 8.5
ArcPy批量计算栅格数据平均值
GLDAS数据下载及处理(NC转TIF)
ArcGIS批量裁剪栅格数据
ArcPy批量栅格重采样
ArcPy批量裁剪栅格数据
前言
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即可。
后记
写博客的初衷是分享我的一些经验,同时也方便自己在其他电脑上进行数据处理。帮了很多人,但评论区小伙伴也有遇到问题的,那么:知识付费,我的时间和经验正好可以解决你的问题。