NumPy array to Raster file (GeoTIFF)

Question

I converted a raster file (one band grey scale) to an numpy array and set a range to it. I Did like that:

import gdal
import numpy
ds = gdal.Open('H:/GDAL_docs/a.tif')
myarray = numpy.array(ds.GetRasterBand(1).ReadAsArray())
selection = numpy.logical_or(myarray >= 0.3557, myarray <= 0.0753)

That works fine. Now I have the problem to convert the numpy array varriable back to a raster file (my favorite format would be GeoTIFF) is that possible? I also struggle to install rasterio to my osgeo4w64 instalation.

Answer 1

By using a raster with integer values (1, 100) and one equivalent condition (myarray >= 35, myarray <= 7), following code would work as expected:

from osgeo import gdal, osr
import numpy

ds = gdal.Open('/home/zeito/pyqgis_data/aleatorio.tif')
band = ds.GetRasterBand(1)
myarray = numpy.array(band.ReadAsArray())
selection = numpy.logical_or(myarray >= 35, myarray <= 7)

new_array = [ [0 for i in range(band.XSize)] for j in range(band.YSize)]

for i, item in enumerate(myarray):
    for j, element in enumerate(item):
        if selection[i][j] == True:
            new_array[i][j] = myarray[i][j]
        else:
            new_array[i][j] = -999

geotransform = ds.GetGeoTransform()
wkt = ds.GetProjection()

# Create gtif file
driver = gdal.GetDriverByName("GTiff")
output_file = "/home/zeito/pyqgis_data/raster_geotiff.tif"

dst_ds = driver.Create(output_file,
                       band.XSize,
                       band.YSize,
                       1,
                       gdal.GDT_Int16)

new_array = numpy.array(new_array)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(-999)
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

ds = None
dst_ds = None

After running above code, I corroborated in several pixels of raster of following image that myarray >= 35, myarray <= 7 condition worked as expected.