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"

In a python script or a unix shell script? What is the pattern of the input and output file names? When you say "manually", do you mean those parameters to -scale_N are set manually in the dialog? Or are they the real max-min of the data? Do you have these 100 images as 100 layers in QGIS? — Spacedman, Oct 03 '19 at 08:03
Python script. The input is a uint16 geotiff with LZW compression and the output should be a uint8 geotiff without compression. I set the parameters in the command prompt. They are the values shown in Layer->Layer Properties->Band rendering (I don't think they are the real max-min of the data since it is uint16). I have 100 images in a directory. — bluetooth12, Oct 03 '19 at 14:53

Comrade Che · Answer 1 · 2019-10-03T12:26:29.377

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

score 1 · Answer 2 · answered Jul 14 '21 at 05:41

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)

Johan · Answer 3 · 2019-10-03T12:10:55.213

0

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)

edited Oct 03 '19 at 12:10

answered Oct 03 '19 at 11:51

Johan

1,141
1
10
19

This is great but how do I load the scales automatically which I currently obtain manually through QGIS GUI? – bluetooth12 Oct 03 '19 at 18:58
One approach would be to use gdal.Info() to determine the statistics (min, max, mean, stdev) of each file and use those values to determine your scale values. – Johan Oct 07 '19 at 07:18
See https://gis.stackexchange.com/a/146989/142299 for using gdal python API to get stats – Nick Jul 02 '20 at 03:36
This won't work as it is. You need to pass shell=True if you're passing a string to subprocess.check_call. Reference: https://stackoverflow.com/questions/18962785/oserror-errno-2-no-such-file-or-directory-while-using-python-subprocess-in-dj – ashnair1 Jan 28 '21 at 12:32

Transforming gdal_translate command into script

3 Answers3

prepare output array

Update metadata

Write as GTiff