image processing - Aggregating neighbouring pixels Python / GDAL and Numpy without interpolation -
consider have image of 2000 x 2000 pixels , pixel size 10 x 10 meters. pixel values float numbers ranging 0.00 - 10.00. image a.
i resize image a quarter of dimensions (i.e. 1000 x 1000 pixels) pixel size 20 x 20 meters (image b) spatially aggregating 4 neighbouring pixels in non-overlapping blocks, starting top-left corner of image, while value each pixel in image b result of arithmetic average.
i have written following code using several sources stackoverflow; reason not understand resulting image (image b) not written , not readable of software want process further (i.e. arcgis, envi, erdas etc).
i appreciate help
best regards dimitris
import time import glob import os import gdal import osr import numpy np start_time_script = time.clock() path_ras='c:/rasters/' 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 size of neighbourhood:' + ' ' + str(sz) h,w = ds_array.shape print 'this size of array:' + ' ' + str(h) + ' ' + str(w) bh, bw = 2,2 shape = (h/bh, w/bw, bh, bw) print 'this new shape of 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() in xrange(len(blocks)): j in xrange(len(blocks[i])): zero_array[i][j] = np.mean(blocks[i][j]) print 'i finished calculations , going write 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'
import time import glob import os import gdal import osr import numpy np start_time_script = time.clock() path_ras='c:/rasters/' 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 size of neighbourhood:' + ' ' + str(sz) h,w = ds_array.shape print 'this size of array:' + ' ' + str(h) + ' ' + str(w) bh, bw = 2,2 shape = (h/bh, w/bw, bh, bw) print 'this new shape of 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() in xrange(len(blocks)): j in xrange(len(blocks[i])): zero_array[i][j] = np.mean(blocks[i][j]) print 'i finished calculations , going write 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'
Comments
Post a Comment