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"

来源：https://blog.csdn.net/hnyzwtf/article/details/44851681