GDAL and PyWavelets

[Jelmer Oosthoek]

I'm new in the field of wavelets. My background is geoscience and geoprocessing using python. I'm interested in wavelets since I saw its application in this publication: http://miplab.epfl.ch/pub/kalbermatten1101.pdf. See f.e. Figure 6.

I therefore would like to decompose elevation data into different scales.

In the PyWavelets zip I saw wp_2d.py and I'm currently trying to use GDAL instead of PIL to read an elevation data raster into a numpy array, then use pywt to decompose and save the results using GDAL.

As the raster files are big I already saw that pywt can run into memory errors.

I currently have something like this:

import numpy
from osgeo import gdalnumeric
from osgeo import gdal


from pywt import WaveletPacket2D


g = gdal.Open("elevation.tif")
arr = g.ReadAsArray()


wp2 = WaveletPacket2D(arr, 'db2', 'sym')


i = 1
prototype_ds = gdal.Open("elevation.tif")
for row in wp2.get_level(2, 'freq'):
    for node in row:
        gdalnumeric.SaveArray( node.data, "output" + str(i) + ".tif", prototype = prototype_ds )
        i += 1

But the problem is that the output1,2,3...16.tif files are 1/16th the size of elevation.tif. So I'm definitely doing something wrong here!

Concerning the memoryerror I already got the following reply from the PyWavelets creator:

- pywt can take single-precision float (32 bit float) numpy arrays as

input, which can halve the memory consumption compared to the default

double-precision float operations.

- even though WP2D constructs the packet tree lazily, it still keeps

the computed results for higher levels in the memory. What you could

try to do is to perform dwt2 decomposition manually only on selected

coefficient arrays while trowing away higher-level input arrays when

possible. It will be less convenient but you may be able to avoid

accumulating multi-level coefficient arrays in the memory.