一堆的遥感影像里求平均（逐个像元），求几个平均。。。

废话不多说，直接上代码

# -*- coding: utf-8 -*-

"""

Created on Sun Apr 14 14:25:16 2019

@author: Administrator

"""

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('F://hourly_maps_raster//liang//2018-01-01_00.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 = 'F:\hourly_maps_raster\liang'

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

for filename in filenames:

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

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

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("F:\\hourly_maps_raster\\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")

求平均修改版

可以批量求平均，但是删减了nodata的条件，你需要对自己的数据清晰明了，没有nodata值

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 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__":

rows,cols,geotransform,projection,noDataValue = Readxy('D:/minxinan/o3idw/0.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 = 'D:/minxinan/o3idw'

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

for ii in range(4):

for filename in filenames[int(ii*(len(filenames)/4)):int((ii+1)*(len(filenames)/4))]:

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

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

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):

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

WriteGTiffFile("D:/minxinan/o3ave/"+str(ii)+".tif", rows, cols, average,geotransform,projection, -9999, GDT_Float32) #写入结果文件