Transforming gdal_translate command into script

Question

I want to run the command below over 100 images in a directory. The problem is that scale_1, scale_2 and scale_3 are currently obtained manually in the QGIS GUI under Layer > Layer Properties > Band rendering. How do I obtain scale_1, scale_2 and scale_3 automatically?

gdal_translate -ot Byte -b 1 -b 2 -b 3 -scale_1 98 854 -scale_2 196 1028 -scale_3 182 1102 -r cubic "20170729_182758_0f28_3B_AnalyticMS_SR.tif" "output.tif"

Answer 1

See Python API. Also this code may help you understand how to use gdal.Translate in python:

import gdal
from gdalconst import GA_ReadOnly

def convert_raster_to_LZW_compression(input_raster):
    basename = os.path.basename(input_raster)
    basename_without_extention, extention = os.path.splitext(basename)
    dirname = os.path.dirname(input_raster)
    output_raster = os.path.join(dirname, f'{basename_without_extention}_LZW_CONVERTED{extention}')

    input_dataset = gdal.Open(input_raster, GA_ReadOnly)
    InfoOptions = gdal.InfoOptions(format='json')
    input_dataset_info = gdal.Info(input_dataset, options=InfoOptions)

    try:
        compression = input_dataset_info['metadata']['IMAGE_STRUCTURE']['COMPRESSION']
    except KeyError:
        compression = None

    if compression == 'LZW':
        print(f'Raster compression is {compression}')
        input_dataset = None  # закрываем файлы
        return 'skip'
    else:
        translateOptions = gdal.TranslateOptions(creationOptions=['COMPRESS=LZW'])
        gdal.Translate(output_raster, input_dataset, options=translateOptions)
        input_dataset = None  # закрываем файлы

        if not os.path.isfile(output_raster):
            print(f'Ошибка! Файл не конвертирован в LZW сжатие:\n{output_raster}')
            return False
        else:
            output_dataset = gdal.Open(output_raster, GA_ReadOnly)
            output_dataset_info = gdal.Info(output_dataset, options=InfoOptions)
            output_dataset = None  # закрываем файлы
            out_compression = output_dataset_info['metadata']['IMAGE_STRUCTURE']['COMPRESSION']
            if out_compression == 'LZW':
                return True
            else:
                print(f'Ошибка! Файл был сохранён, но не в LZW сжатие:\n{output_raster}\nУдаляю файл...', end='')
                os.remove(output_raster)
                print('OK')
                return False

Answer 2

I recommend to use rasterio if you want to process in python. I didn't test the following code with multiband raster, but I hope the code will work with your data.

This code stretch raster values 0 to 254 range, set nodata 255, dtype Byte (uint8), and save without compression.

import rasterio
import numpy as np
import numpy.ma as ma

with rasterio.open('data/R022_32TQQ_20200309_B02.tif') as src:
    array = src.read()
    profile = src.profile

# prepare output array
output_array = np.empty(array.shape, dtype=np.uint8)

for idx, band_array in enumerate(array):
    masked_array = ma.masked_array(band_array,
                                   mask=(band_array==profile['nodata']))
    min_ = masked_array.min()
    max_ = masked_array.max()
    
    # stretch value in 0-254 range
    output_array[idx] = ((band_array - min_) / (max_ - min_) * 254).round().astype(np.uint8)

    # set 255 at nodata pixel
    output_array[idx][masked_array.mask] = 255

# Update metadata
profile.pop('compress', None)
profile.update(nodata=255, dtype='uint8')

# Write as GTiff
with rasterio.open('tmp/out.tif', 'w', **profile) as dst:
    dst.write(output_array)

Answer 3

In python you can use subprocess to call gdal_translate from the command line:

import subprocess
import os

basecmd = "gdal_translate -ot Byte -b 1 -b 2 -b 3 -scale_1 98 854 -scale_2 196 1028 -scale_3 182 1102 -r cubic"
indirectory = '/path/to/in/directory'
outdirectory = '/path/to/out/directory'
infiles = [x for x in os.listdir(indirectory) if x.endswith('.tif')]

for fn in infiles:
  infn = os.path.join(indirectory, fn)
  outfn = os.path.join(outdirectory, fn)
  if not os.path.exists(outfn):
    cmd = basecmd + " " + infn + " " + outfn
    subprocess.check_call(cmd)