Search code examples
pythonnumpygeotiff

How to convert numpy list to geoTIFF in python?


I haven't found answer on this simple question. Please help. How to convert Qcal (numpy list) to TIFF image ? Everything I've found doesn't really work.

import math
import numpy as np
from osgeo import gdal

substr1 = 'RADIANCE_MULT_BAND_10'
substr2 = 'RADIANCE_ADD_BAND_10'
substr3 = 'K1_CONSTANT_BAND_10'
substr4 = 'K2_CONSTANT_BAND_10'

RADIANCE_MULT_BAND_10 = 1
RADIANCE_ADD_BAND_10 = 1
K1_CONSTANT_BAND_10 = 1
K2_CONSTANT_BAND_10 = 1

with open('LC08_L1TP_180028_20170623_20170630_01_T1_MTL.txt') as file:
    for line in file:
        if substr1 in line:
            startIndex = line.find('=')
            RADIANCE_MULT_BAND_10 = float((line[startIndex+2:]))
        if substr2 in line:
            startIndex = line.find('=')
            RADIANCE_ADD_BAND_10 = float((line[startIndex + 2:]))
        if substr3 in line:
            startIndex = line.find('=')
            K1_CONSTANT_BAND_10 = float((line[startIndex + 2:]))
        if substr4 in line:
            startIndex = line.find('=')
            K2_CONSTANT_BAND_10 = float((line[startIndex + 2:]))

ds = gdal.Open("B10.tif")
Qcal = np.array(ds.GetRasterBand(1).ReadAsArray()) # Quantized and calibrated standard product pixel values (DN)
for i in range(Qcal.shape[0]):
   for j in range(Qcal.shape[1]):
       Qcal[i][j] = RADIANCE_MULT_BAND_10 * Qcal[i][j] + RADIANCE_ADD_BAND_10
       Qcal[i][j] = K2_CONSTANT_BAND_10 / math.log1p(K1_CONSTANT_BAND_10/Qcal+1)

Solution

  • Do you want to modify the existing image?

    If so you should open it with the update option like this

    gdal.Open("B10.tif", gdal.GA_Update)
    

    Next do your modifications to your np array as you are doing. You are actually just editing the numpy array Qcal and not the actual rasterband in the tif.

    Now to save your modifications into the same band you can do the following

    ds.GetRasterBand(1).WriteArray(Qcal)
    

    This is writing the updated Qcal array into the rasterband of the tif.

    If you want to save as a new image

    You can create a copy of the existing image and save the Qcal array into it like this

    driver = gdal.GetDriverByName('Gtiff')
    dst_ds = driver.CreateCopy("example.tif", ds, 1)
    dst_ds.GetRasterBand(1).WriteArray(Qcal)