Python計算多幅圖像柵格值的平均值
本文實例為大傢分享瞭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"
以上就是本文的全部內容,希望對大傢的學習有所幫助,也希望大傢多多支持WalkonNet。
推薦閱讀:
- python多線程方法詳解
- 使用python進行nc轉tif的3種情況解決
- python 使用GDAL實現柵格tif轉矢量shp的方式小結
- 基於Python實現nc批量轉tif格式
- Python讀取hdf文件並轉化為tiff格式輸出