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
Best Answer
Can be done fairly easy. First, create the dataset:
You then need to set the geotransform using
ds.SetGeoTransformwhere the argument is a six element tuple:(upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)Set the projection using
ds.SetProjection, argument is the way string, then Write the array for the first band:That should be enough to get you going. If the purpose is just to reproject the data, then you could look at using an in memory dataset