2

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"
2
  • 2
    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 3 '19 at 8: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 3 '19 at 14:53
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)
4
  • This is great but how do I load the scales automatically which I currently obtain manually through QGIS GUI? – bluetooth12 Oct 3 '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 7 '19 at 7:18
  • See gis.stackexchange.com/a/146989/142299 for using gdal python API to get stats – Nick Jul 2 '20 at 3: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: stackoverflow.com/questions/18962785/… – Ashwin Nair Jan 28 at 12:32
0

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
0

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)

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.