Delete raster values and compress multiple raster files

Question

I have hundreds of raster files. They have a "band 1" and contain values from 0 - 9000, however only values between 0 and 255 are of interest. In order to make the files smaller I'm thinking to delete all data with values >255 and then use LZW compression to make it even smaller. I'm not sure about how to do this though.

I have made a python script like this:

myDir = '/output/'
layers = QgsProject.instance().mapLayers()
pipe = QgsRasterPipe()
for layer in layers.values():
    print(layer.name())
    opts = ["COMPRESS=LZW"]
    extent = layer.extent()
    width, height = layer.width(), layer.height()
    renderer = layer.renderer()
    provider=layer.dataProvider()
    provider.setUserNoDataValue(1, [QgsRasterRange(256,999999)])
    crs = layer.crs().toWkt() 
    pipe.set(provider.clone())
    pipe.set(renderer.clone())
    file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
    file_writer.setCreateOptions(opts)
    file_writer.writeRaster(pipe,
                            width,
                            height,
                            extent,
                            layer.crs())

The visible file in my qgis workspace gets values 0-255 and appears to be correct, but the saved file in /output/ folder has now 3 bands and not correct values. If I manually save the file in my workspace without changing anything, it turns out to be correct. But the python-saved files turns out incorrect.

What am I doing wrong? (I use QGIS 3.4)

Answer 1

Eventually I figures out how to do it with pyQgis:

myDir = '/output/'
style = '/stile_file.qml'

layers = QgsProject.instance().mapLayers()
for layer in layers.values():
    layer.loadNamedStyle(style)
    opts = ["COMPRESS=deflate"]
    extent = layer.extent()
    width, height = layer.width(), layer.height()

    renderer = layer.renderer()
    provider=layer.dataProvider()
    provider.setUserNoDataValue(1, [QgsRasterRange(256,999999)])
    crs = layer.crs().toWkt()

    pipe = QgsRasterPipe()
    pipe.set(provider.clone())
    pipe.set(renderer.clone())
    file_writer = QgsRasterFileWriter(myDir + layer.name() + '.tif')
    file_writer.setCreateOptions(opts)
    file_writer.writeRaster(pipe,
                            width,
                            height,
                            extent,
                            layer.crs())

Answer 2

You can accomplish this with gdal (which should already be installed if you have QGIS on your machine). Furthermore, I would suggest changing the raster data type so the size of the raster file is smaller. To store values from 0 to 9000, the data type of your raster must be at least an (un)signed 16-bit integer. However, if you want to store values from 0 to 255, an unsigned 8-bit integer is enough and will save some space. Take into account that the NoData value (which is going to be assigned to those values you "delete") has to be inside this range.

Here is a small snippet that takes an input raster file and creates a smaller version by adding LZW compression and changing the raster dataset. Note that changing the data type can't be done in place and a copy must be created.

from osgeo import gdal

src = 'path/to/input.tif'
dst = 'path/to/output.tif'
nodata = 1

# Open input raster and get misc. information
in_ds = gdal.Open(src)
xsize = in_ds.RasterXSize
ysize = in_ds.RasterYSize
bands = in_ds.RasterCount
sr = in_ds.GetProjection()
gt = in_ds.GetGeoTransform()

# Create empty raster with same properties as the original (except the data type)
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(dst, xsize, ysize, bands, gdal.GDT_Byte, options=['COMPRESS=LZW'])
out_ds.SetProjection(sr)
out_ds.SetGeoTransform(gt)

# Write the original values truncated to the 0-255 range
for i in range(bands):
    in_band = in_ds.GetRasterband(i+1)
    out_band = in_ds.GetRasterBand(i+1)
    data = in_band.ReadAsArray()
    data[data > 255] = nodata
    out_band.WriteArray(arr)
    out_band.FlushCache()

# Save output file
out_ds = None

I'll leave it to you to adapt this snippet to work with multiple files but a good start is the os.walk() function if you have nested folders. Furthermore, you might want to delete the original raster after creating the smaller copy.