Блог

У меня есть кратный растр с разрешением в пикселях 10 X 10 метров, а значение растра в формате float, и я хочу сохранить исходное значение растра только в формате float.

Я хочу рассчитать усреднение пикселей окна 5 x 5, чтобы удалить шум с изображений и создать гладкий растр в Python. Я разработал приведенный ниже код из разных источников stackoverflow, но он не работает и выдает ошибку.

Я был бы признателен за любую помощь.

Приведенный ниже код:

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

start_time_script = time.clock()

path_ras=r'D:Firm_SMF1A/'

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(4726)

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 = 5,5

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')
  

Массаж ошибок:

 This is the size of the neighbourhood: 4
This is the size of the Array: 106 144
This is the new shape of the Array: (21.2, 28.8, 5, 5)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-fbee5d0233b0> in <module>
     21 strides = sz*np.array([w*bh,bw,w,1])
     22 
---> 23 blocks = np.lib.stride_tricks.as_strided(ds_array,shape=shape,strides=strides)
     24 
     25 resized_array = ds_driver.Create(rasterfile_name   '_resized_to_52m.tif',shape[1],shape[0],1,gdal.GDT_Float32)

c:python37libsite-packagesnumpylibstride_tricks.py in as_strided(x, shape, strides, subok, writeable)
    101         interface['strides'] = tuple(strides)
    102 
--> 103     array = np.asarray(DummyArray(interface, base=x))
    104     # The route via `__interface__` does not preserve structured
    105     # dtypes. Since dtype should remain unchanged, we set it explicitly.

c:python37libsite-packagesnumpycore_asarray.py in asarray(a, dtype, order)
     83 
     84     """
---> 85     return array(a, dtype, copy=False, order=order)
     86 
     87 

TypeError: 'float' object cannot be interpreted as an integer