2

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)

2
  • I would suggest changing the raster data type as well. An unsigned byte can store values from 0 to 255, which is exactly your value range. Unfortunately, I am not familiar with the QGIS python bindings so can't help you with the code, unless. However, if you are open to use gdal in your Python script, I might help you out. – Marcelo Villa-Piñeros Aug 28 '20 at 17:15
  • @MarceloVilla-Piñeros, I'll take whatever help you have to offer. Just need it to be automated for all files in the folder/working tree, to avoid spending too much time. What would be the preferred raster data type? – agh Aug 31 '20 at 2:04
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())
0

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.

1
  • 1
    Thank you very much Marcelo! I eventually solved it with pyQgis and posted my answer. So now future users have two alternatives to solve this issue. – agh Sep 4 '20 at 3:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.