Pythonは、複数の画像のバリア値の平均値を計算します。
本論文の例では、Pythonが複数の画像のバリア値の平均値を求めていることを共有しています。
本プログラムで採用した方法は最適な方法ではなく、ARCGISは呼び出しのために関連する関数を提供しました。このプログラムは参考のみです。
プログラムの説明:
フォルダE://work/EVI_Data_tifに保存されているのは、ある地域の2000年から2010年までのEVI画像で、その中の年ごとに13枚です。毎年13枚の画像の各格子を加算して平均値を求め、その年のtifを生成するのが目的です。例えば、2000年の13枚の画像を加算して平均値を求めて20000.tifを生成し、各バリアの値は13枚の画像に対応するグリッド値を加算して得られた平均値である。結果はE:\work\resultに保存されます。ソースファイルの組織方式は2000年を例にして、ファイル名は順次20006.tif、20007.tif、20008.tif、…、200018 tifです。
本プログラムで採用した方法は最適な方法ではなく、ARCGISは呼び出しのために関連する関数を提供しました。このプログラムは参考のみです。
プログラムの説明:
フォルダE://work/EVI_Data_tifに保存されているのは、ある地域の2000年から2010年までのEVI画像で、その中の年ごとに13枚です。毎年13枚の画像の各格子を加算して平均値を求め、その年のtifを生成するのが目的です。例えば、2000年の13枚の画像を加算して平均値を求めて20000.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"