Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif.
From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it in QGIS. If I remember correctly, when you first add a raster, it often shows up all one color if you don't have a pre-existing color map.
import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
print 'Could not open image file'
sys.exit(1)
# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize
cropData = band1.ReadAsArray(0,0,cols,rows)
listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]
# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
print 'Could not create reclass_40.tif'
sys.exit(1)
outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)
for i in range(0, rows):
for j in range(0, cols):
if cropData[i,j] in listAg:
outData[i,j] = 100
elif cropData[i,j] in listNotAg:
outData[i,j] = -100
else:
outData[i,j] = 0
# write the data
outBand.WriteArray(outData, 0, 0)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
del outData
I managed to do it in the end, so I'll post the answer here in case anyone needs it later.
The problem was that my image contained float32 and sometimes even complex32 data, but QImage does not support it, so I had to convert it to uint8.
Another problem that came up was that my image was completely black, because the values were most of the time like 0.x, so it was considered as black. I tried to use ImageEnhance from the PIL library, but the result wasn't satisfying due to some unusual values that would trouble the results. So I just decided to multiply every value of the array by a factor, determined by the user. Depending on the values you have in your array, it is not unusual to set the factor really high, for example in my case I usually set it at 300.
The code looks like this:
# Import the libraries
from osgeo import gdal
from PIL import Image, ImageQt
from PyQt4 import QtGui
import numpy as np
# Read the dataset and setting up the numpy arrays
dsr = gdal.Open("path/to/file")
np_array = np.array(dsr.ReadAsArray())
np_array_uint8 = (np_array * factor).astype(np.uint8) # 'factor' is a user-determined value
# Convert to image
im8 = Image.fromarray(np_array_uint8)
imQt = QtGui.QImage(ImageQt.ImageQt(im8))
pxmap = QtGui.QPixmap.fromImage(imQt)
# Set the pixmap into a label
self.label = QtGui.QLabel(self.centralwidget)
self.label.setPixmap(QtGui.QPixmap(mr_img)) # You can scale it as well
There you go, in case you ever need it.
You can also save memory by setting everything to None if you need.
Hope this helps!
Best Answer
By using a raster with integer values (1, 100) and one equivalent condition (myarray >= 35, myarray <= 7), following code would work as expected:
After running above code, I corroborated in several pixels of raster of following image that myarray >= 35, myarray <= 7 condition worked as expected.