0

I'm currently working with GDAL and numpy in Python to stitch 2 parts of different GeoTIFF files and produce a new GeoTIFF with the result.

I'm able to retrieve the files, transform it in numpy arrays and modify it as needed including the nodata value. The problem arises when I try to convert the numpy array to a GeoTIFF file. All the transformation goes OK but when I try to set the no data value the tif file produced is not readable by QGIS indicating that the file is not supported (same with other tif visors). If I delete the nodata value instructions works OK but the nodata value obviously is not included. Input files can be float32, uint16 or int16 but both files have the same variable type.

I'm missing something?

Here is the code used

import numpy
import tkinter
from osgeo import gdal, osr, ogr
from numpy import savetxt
from numpy import concatenate
from numpy import shape
from tkinter import filedialog as fd
filename1 = fd.askopenfilename(title="NORTH")
filename2 = fd.askopenfilename(title="SOUTH")
print("NORTH FILE =",filename1)
print("SOUTH FILE =",filename2)
filename3 = fd.asksaveasfile(initialfile = 'Untitled.tif', defaultextension=".tif",filetypes=[("All Files","*.*"),("Text Documents","*.txt")])
print("OUTPUT FILE=",filename3.name)
R01 = gdal.Open(filename1)
#------------------------------------------------------------------
# Here I rescue the nodatavalue from the input file
# -----------------------------------------------------------------
nulval= R01.GetRasterBand(1).GetNoDataValue()
RA1 = R01.ReadAsArray()
R02 = gdal.Open(filename2)
RA2 = R02.ReadAsArray()
print("Extract north & south")
NOR= RA1[0:10800,0:43200]
SUR=RA2[10800:21600,0:43200]
print("files divided")
print("Concatenate")
SUN1=concatenate((NOR,SUR))
print("concatenate ready")
if tipo == "float32": tipoarr= gdal.GDT_Float32
if tipo == "uint16": tipoarr= gdal.GDT_Uint16
if tipo == "int16": tipoarr= gdal.GDT_Int16
print("output tif file")
driver = gdal.GetDriverByName("GTiff")
SAL01=driver. Create(filename3.name,SUN1.shape[1],SUN1.shape[0],1,tipoarr)
SAL01.SetProjection(R01.GetProjection())
SAL01.SetGeoTransform(R01.GetGeoTransform())
#-------------------------------------------------------
# Here I set the same nodatavalue to the output file
# ------------------------------------------------------
SAL01.GetRasterBand(1).SetNoDataValue(nulval)
band = SAL01.GetRasterBand(1)
band.WriteArray(SUN1)
band.FlushCache()
band.ComputeStatistics(False)
print("END PROCESS")
1
  • Problem solved.
    – Eduardo
    Aug 3, 2023 at 13:04

1 Answer 1

0

Problem solved.

Appears that the files need to be closed (Var=None) at the end of the process to produce a valid GeoTIFF file.

End of the code attached.

Cheers!

band.WriteArray(SUN1)
band.FlushCache()
band.ComputeStatistics(False)
SAL01=None
R01=None
R02=None
print("FIN")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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