Python calculates the average of grid values of multiple images

In this paper, we share the average value of Python for grid values of multiple images for your reference. The specific contents are as follows

The method used in this program is not the optimal method, ARCGIS has provided related functions for calling. This procedure is for reference only.

Procedure description:

The folder E://work//EVI_Data_tif stores EVI images of a certain area from 2000 to 2010, including 13 images in each year. The purpose is to average each grid of 13 images per year to generate tif for the corresponding year. For example, 13 images from 2000 are added and averaged to generate 2000. tif, in which the value of each grid is the average value obtained by adding the corresponding grid values of 13 images. The results are stored in E:\ work\ result. The source files are organized as follows: Take 2000 as an example, and the file names are 20006. tif, 20007. tif, 20008. tif, …, 200018. tif in turn.

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):# Write the result file to disk 
    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():# Traverse the file, read the data, and calculate the average value 
    rows,cols,geotransform,projection,noDataValue = Readxy('E://work//EVI_Data_tif//20006.tif')
    # Get the row, column, projection and other information of the source file. All the information of the source file is 1 To 
    print 'rows and cols is ',rows,cols
    filesum = [[0.0]*cols]*rows # Grid values and, 2 Dimensional array 
    average= [[0.0]*cols]*rows#  Store the average value, 2 Dimensional array 
    filesum=np.array(filesum)# The conversion type is 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):# Traverse the source file 
        for filename in filenames:
            if os.path.splitext(filename)[1] == '.tif':# Determine whether it is tif Format 
                filepath = os.path.join(dirpath,filename)
                purename = filename.replace('.tif','') # Get the file name with the extension removed, such as 201013.tif , purename For 201013               
                if purename[:4] == '2010':              # Judgement year 
                    filedata = [[0.0]*cols]*rows
                    filedata = np.array(filedata)
                    filedata = Read(filepath)           # Will 2010 Year 13 Image data is stored in filedata Medium 
                    count+=1
                    np.add(filesum,filedata,filesum)    # Seek 13 The sum of the corresponding grid values of the images 
                    #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):        # Process the in the image noData
                average[i,j]=-9999
            else: 
                average[i,j]=filesum[i,j]*1.0/count     # Average 
    WriteGTiffFile("E:\\work\\result\\2010.tif", rows, cols, average,geotransform,projection, -9999, GDT_Float32) # Write the result file            

def Readxy(RasterFile): # Read the information of each image      
    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):# Read the information of each image 
    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"