Python中如何計算多幅圖像柵格值的平均值

這篇文章主要為大家展示了“Python中如何計算多幅圖像柵格值的平均值”，內(nèi)容簡而易懂，條理清晰，希望能夠幫助大家解決疑惑，下面讓小編帶領(lǐng)大家一起研究并學(xué)習(xí)一下“Python中如何計算多幅圖像柵格值的平均值”這篇文章吧。

程序說明：

文件夾E://work//EVI_Data_tif中存放的是某地區(qū)2000-2010年的EVI圖像，其中每個年份共13幅。目的是將每年的13幅圖像的每個柵格相加求均值，生成相應(yīng)年份的tif。例如，將2000年的13幅圖像相加求均值生成2000.tif，里面的每個柵格的值就是13幅圖像對應(yīng)柵格值相加得到的均值。結(jié)果存放于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):#向磁盤寫入結(jié)果文件
    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():#遍歷文件，讀取數(shù)據(jù)，算出均值
    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 #柵格值和，二維數(shù)組
    average= [[0.0]*cols]*rows# 存放平均值，二維數(shù)組
    filesum=np.array(filesum)#轉(zhuǎn)換類型為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','') #獲得除去擴(kuò)展名的文件名，比如201013.tif，purename為201013               
                if purename[:4] == '2010':              #判斷年份
                    filedata = [[0.0]*cols]*rows
                    filedata = np.array(filedata)
                    filedata = Read(filepath)           #將2010年的13幅圖像數(shù)據(jù)存入filedata中
                    count+=1
                    np.add(filesum,filedata,filesum)    #求13幅圖像相應(yīng)柵格值的和
                    #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) #寫入結(jié)果文件           

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"

以上是“Python中如何計算多幅圖像柵格值的平均值”這篇文章的所有內(nèi)容，感謝各位的閱讀！相信大家都有了一定的了解，希望分享的內(nèi)容對大家有所幫助，如果還想學(xué)習(xí)更多知識，歡迎關(guān)注億速云行業(yè)資訊頻道！