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.

• GeoMonkey: Take a look here: gis.stackexchange.com/questions/290776/…
• Aaron: Could you please edit your question to include what approach you tried converting from numpy to tiff and the problem you encountered?

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.

enter image description here

• sameer_nubia: Hi, What is the important of this line (myarray >= 35, myarray <= 7), ?