欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

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

程序员文章站 2022-03-10 20:54:14
本文实例为大家分享了python求多幅图像栅格值的平均值,供大家参考,具体内容如下本程序所采用的方法并不是最优方法,arcgis已经提供了相关的函数供调用。本程序仅供参考。程序说明:文件夹e://wo...

本文实例为大家分享了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"

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持。