本文实例为大家分享了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。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
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