Python计算多幅图像栅格值的平均值_Python

本文实例为大家分享了Python求多幅图像栅格值的平均值，供大家参考，具体内容如下

本程序所采用的方法并不是最优方法，ARCGIS已经提供了相关的函数供调用。本程序仅供参考。

程序说明：

文件夹E://work//EVI_Data_tif中存放的是某地区2000-2010年的EVI图像，其中每个年份共13幅。目的是将每年的13幅图像的每个栅格相加求均值，生成相应年份的tif。例如，将2000年的13幅图像相加求均值生成2000.tif，里面的每个栅格的值就是13幅图像对应栅格值相加得到的均值。结果存放于E:\work\result。源文件组织方式为：以2000年为例，文件名依次为 20006.tif，20007.tif，20008.tif，……，200018.tif。

				?

									import os

									import os.path

									import gdal

									import sys

									from gdalconst import *

									from osgeo import gdal

									import osr

									import numpy as np

									#coding=utf-8

									def WriteGTiffFile(filename, nRows, nCols, data,geotrans,proj, noDataValue, gdalType):#向磁盘写入结果文件

									    format = "GTiff"  

									    driver = gdal.GetDriverByName(format)

									    ds = driver.Create(filename, nCols, nRows, 1, gdalType)

									    ds.SetGeoTransform(geotrans)

									    ds.SetProjection(proj)

									    ds.GetRasterBand(1).SetNoDataValue(noDataValue)

									    ds.GetRasterBand(1).WriteArray(data)    

									    ds = None

									def File():#遍历文件，读取数据，算出均值

									    rows,cols,geotransform,projection,noDataValue = Readxy('E://work//EVI_Data_tif//20006.tif')

									    #获取源文件的行，列，投影等信息，所有的源文件这些信息都是一致的

									    print 'rows and cols is ',rows,cols

									    filesum = [[0.0]*cols]*rows #栅格值和，二维数组

									    average= [[0.0]*cols]*rows# 存放平均值，二维数组

									    filesum=np.array(filesum)#转换类型为np.array

									    average = np.array(average) 

									    print 'the type of filesum',type(filesum)

									    count=0

									    rootdir = 'E:\work\EVI_Data_tif'

									    for dirpath,filename,filenames in os.walk(rootdir):#遍历源文件

									        for filename in filenames:

									            if os.path.splitext(filename)[1] == '.tif':#判断是否为tif格式

									                filepath = os.path.join(dirpath,filename)

									                purename = filename.replace('.tif','') #获得除去扩展名的文件名，比如201013.tif，purename为201013               

									                if purename[:4] == '2010':              #判断年份

									                    filedata = [[0.0]*cols]*rows

									                    filedata = np.array(filedata)

									                    filedata = Read(filepath)           #将2010年的13幅图像数据存入filedata中

									                    count+=1

									                    np.add(filesum,filedata,filesum)    #求13幅图像相应栅格值的和

									                    #print str(count)+'this is filedata',filedata

									    print 'count is ',count    

									    for i in range(0,rows):

									        for j in range(0,cols):

									            if(filesum[i,j]==noDataValue*count):        #处理图像中的noData

									                average[i,j]=-9999

									            else: 

									                average[i,j]=filesum[i,j]*1.0/count     #求平均

									    WriteGTiffFile("E:\\work\\result\\2010.tif", rows, cols, average,geotransform,projection, -9999, GDT_Float32) #写入结果文件           

									def Readxy(RasterFile): #读取每个图像的信息     

									    ds = gdal.Open(RasterFile,GA_ReadOnly)

									    if ds is None:

									        print 'Cannot open ',RasterFile

									        sys.exit(1)

									    cols = ds.RasterXSize

									    rows = ds.RasterYSize

									    band = ds.GetRasterBand(1)

									    data = band.ReadAsArray(0,0,cols,rows)

									    noDataValue = band.GetNoDataValue()

									    projection=ds.GetProjection()

									    geotransform = ds.GetGeoTransform()

									    return rows,cols,geotransform,projection,noDataValue

									def Read(RasterFile):#读取每个图像的信息

									    ds = gdal.Open(RasterFile,GA_ReadOnly)    

									    if ds is None:

									        print 'Cannot open ',RasterFile

									        sys.exit(1)

									    cols = ds.RasterXSize

									    rows = ds.RasterYSize

									    band = ds.GetRasterBand(1)

									    data = band.ReadAsArray(0,0,cols,rows)  

									    return data    

									if __name__ == "__main__":

									    print"ok1"

									    File()   

									    print"ok2"