2

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.

2
  • Take a look here: gis.stackexchange.com/questions/290776/…
    – GeoMonkey
    Commented Jan 8, 2019 at 12:57
  • Could you please edit your question to include what approach you tried converting from numpy to tiff and the problem you encountered?
    – Aaron
    Commented Jan 8, 2019 at 15:05

1 Answer 1

6

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

1
  • Hi, What is the important of this line (myarray >= 35, myarray <= 7), ? Commented Nov 3, 2020 at 13:19

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.