聚合相鄰像素Python/GDAL和Numpy無插值

Question

考慮我們有一個2000 x 2000像素的圖像和像素大小爲10 x 10米。像素值也是範圍從0.00到10.00的浮點數。這是圖像A.

我想通過在非空間聚合四個相鄰像素，將圖像A的大小調整爲像素大小爲20 x 20米（圖像B）的四分之一尺寸（即1000 x 1000像素）重疊塊，從圖像的左上角開始，而圖像B中每個像素的值將是其算術平均值的結果。

我寫了下面的代碼使用來自stackoverflow的幾個來源;但由於某些原因，我不理解所得到的圖像（圖像B）並不總是正確書寫，並且它不能被我想進一步處理它的任何軟件（即ArcGIS，ENVI，ERDAS等）讀取。

我希望得到任何幫助

問候 季米特里斯

import time 
import glob 
import os 
import gdal 
import osr 
import numpy as np 

start_time_script = time.clock() 

path_ras='C:/rasters/' 

for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')): 
rasterfile_name=str(rasterfile[rasterfile.find('IMG'):rasterfile.find('.tif')]) 

print 'Processing:'+ ' ' + str(rasterfile_name) 

ds = gdal.Open(rasterfile,gdal.GA_ReadOnly) 
ds_xform = ds.GetGeoTransform() 

print ds_xform 

ds_driver = gdal.GetDriverByName('Gtiff') 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(26716) 

ds_array = ds.ReadAsArray() 

sz = ds_array.itemsize 

print 'This is the size of the neighbourhood:' + ' ' + str(sz) 

h,w = ds_array.shape 

print 'This is the size of the Array:' + ' ' + str(h) + ' ' + str(w) 

bh, bw = 2,2 

shape = (h/bh, w/bw, bh, bw) 

print 'This is the new shape of the Array:' + ' ' + str(shape) 

strides = sz*np.array([w*bh,bw,w,1]) 

blocks = np.lib.stride_tricks.as_strided(ds_array,shape=shape,strides=strides) 

resized_array = ds_driver.Create(rasterfile_name + '_resized_to_52m.tif',shape[1],shape[0],1,gdal.GDT_Float32) 
resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*2,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*2)) 
resized_array.SetProjection(srs.ExportToWkt()) 
band = resized_array.GetRasterBand(1) 

zero_array = np.zeros([shape[0],shape[1]],dtype=np.float32) 

print 'I start calculations using neighbourhood' 
start_time_blocks = time.clock() 

for i in xrange(len(blocks)): 
    for j in xrange(len(blocks[i])): 

     zero_array[i][j] = np.mean(blocks[i][j]) 

print 'I finished calculations and I am going to write the new array' 

band.WriteArray(zero_array) 

end_time_blocks = time.clock() - start_time_blocks 

print 'Image Processed for:' + ' ' + str(end_time_blocks) + 'seconds' + '\n' 

end_time = time.clock() - start_time_script 
print 'Program ran for: ' + str(end_time) + 'seconds' 

Answer 1

import time 
import glob 
import os 
import gdal 
import osr 
import numpy as np 

start_time_script = time.clock() 

path_ras='C:/rasters/' 

for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')): 
rasterfile_name=str(rasterfile[rasterfile.find('IMG'):rasterfile.find('.tif')]) 

print 'Processing:'+ ' ' + str(rasterfile_name) 

ds = gdal.Open(rasterfile,gdal.GA_ReadOnly) 
ds_xform = ds.GetGeoTransform() 

print ds_xform 

ds_driver = gdal.GetDriverByName('Gtiff') 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(26716) 

ds_array = ds.ReadAsArray() 

sz = ds_array.itemsize 

print 'This is the size of the neighbourhood:' + ' ' + str(sz) 

h,w = ds_array.shape 

print 'This is the size of the Array:' + ' ' + str(h) + ' ' + str(w) 

bh, bw = 2,2 

shape = (h/bh, w/bw, bh, bw) 

print 'This is the new shape of the Array:' + ' ' + str(shape) 

strides = sz*np.array([w*bh,bw,w,1]) 

blocks = np.lib.stride_tricks.as_strided(ds_array,shape=shape,strides=strides) 

resized_array = ds_driver.Create(rasterfile_name + '_resized_to_52m.tif',shape[1],shape[0],1,gdal.GDT_Float32) 
resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*2,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*2)) 
resized_array.SetProjection(srs.ExportToWkt()) 
band = resized_array.GetRasterBand(1) 

zero_array = np.zeros([shape[0],shape[1]],dtype=np.float32) 

print 'I start calculations using neighbourhood' 
start_time_blocks = time.clock() 

for i in xrange(len(blocks)): 
    for j in xrange(len(blocks[i])): 

     zero_array[i][j] = np.mean(blocks[i][j]) 

print 'I finished calculations and I am going to write the new array' 

band.WriteArray(zero_array) 

end_time_blocks = time.clock() - start_time_blocks 

print 'Image Processed for:' + ' ' + str(end_time_blocks) + 'seconds' + '\n' 

end_time = time.clock() - start_time_script 
print 'Program ran for: ' + str(end_time) + 'seconds'