Writing numpy array to raster file

I'm new to GIS, and I'm lost with something I feel ought to be relatively simple. I have some code that converts infrared images of Mars into thermal inertia maps, which are then stored as 2D numpy arrays. I've been saving these maps as hdf5 files but I'd really like to save them as raster images so that I can process them in QGIS. I've gone through multiple searches to find how to do this but with no luck. I've tried following the instructions in the tutorial at http://www.gis.usu.edu/~chrisg/python/ but the files I produce using his example code open as plain grey boxes when I import them to QGIS. I feel like if someone could suggest the simplest possible procedure to a simplified example of what I'd like to do then I might be able to make some progress. I have QGIS and GDAL, I'd be very happy to install other frameworks that anyone could recommend. I use Mac OS 10.7.

So if for example I have a numpy array of thermal inertia that looks like:

``````TI = ( (0.1, 0.2, 0.3, 0.4),
(0.2, 0.3, 0.4, 0.5),
(0.3, 0.4, 0.5, 0.6),
(0.4, 0.5, 0.6, 0.7) )
``````

And for each pixel I have the latitude and longitude:

``````lat = ( (10.0, 10.0, 10.0, 10.0),
( 9.5,  9.5,  9.5,  9.5),
( 9.0,  9.0,  9.0,  9.0),
( 8.5,  8.5,  8.5,  8.5) )
lon = ( (20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5) )
``````

Which procedure would people recommend to convert this data into a raster file that I can open in QGIS? Or is there something I'm misunderstanding about the capabilities of GIS.

-
Which slide on the tutorial are you referring to? – R.K. Oct 22 '12 at 11:22

One possible solution to your problem: Convert it into a ASCII Raster, documention for which is here. This should be fairly easy to do with python.

So with your example data above, you'd end up with the following in a .asc file:

``````ncols 4
nrows 4
xllcorner 20
yllcorner 8.5
cellsize 0.5
nodata_value -9999
0.1 0.2 0.3 0.4
0.2 0.3 0.4 0.5
0.3 0.4 0.5 0.6
0.4 0.5 0.6 0.7
``````

This successfully adds to both QGIS and ArcGIS, and stylised in ArcGIS it looks like this:

Addendum: While you can add it to QGIS as noted, if you try and go into the properties for it (to stylise it), QGIS 1.8.0 hangs. I'm about to report that as a bug. If this happens to you too, then there are plenty of other free GIS's out there.

-
That's fantastic, thanks. And I imagine that having written my array as an ascii file I could convert it into a binary format using a pre-written conversion function. – EddyThe B Oct 22 '12 at 19:03
FYI, I didn't have the hanging issue with QGIS, I have version 1.8.0 too. – EddyThe B Oct 22 '12 at 19:10

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)

band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

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
``````
-
+1 for flushing -- was banging my head against the wall trying to figure out how to 'save' the thing! – badgley Aug 7 at 16:46

Thanks everyone for your help. I've finally hit upon this solution, which I gained from this discussion (http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html). I like it because I can go straight from a numpy array to a tif raster file. I'd be very grateful for comments that could improve on the solution. I'll post it here in case anyone else searches for a similar answer.

``````import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

array = np.array(( (0.1, 0.2, 0.3, 0.4),
(0.2, 0.3, 0.4, 0.5),
(0.3, 0.4, 0.5, 0.6),
(0.4, 0.5, 0.6, 0.7),
(0.5, 0.6, 0.7, 0.8) ))
# My image array
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
( 9.5,  9.5,  9.5,  9.5),
( 9.0,  9.0,  9.0,  9.0),
( 8.5,  8.5,  8.5,  8.5),
( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5),
(20.0, 20.5, 21.0, 21.5) ))
# For each pixel I know it's latitude and longitude.
# As you'll see below you only really need the coordinates of
# one corner, and the resolution of the file.

xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()]
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
# That's (top left x, w-e pixel resolution, rotation (0 if North is up),
#         top left y, rotation (0 if North is up), n-s pixel resolution)
# I don't know why rotation is in twice???

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
# Anyone know how to specify the
# IAU2000:49900 Mars encoding?
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system
# to the file
output_raster.GetRasterBand(1).WriteArray(array)   # Writes my array to the raster
``````
-
The "rotation is in twice" to account for the effect of a rotated bit of y on x and the rotated bit of x on y. See lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html which tries to explain the interrelationships between the "rotation" parameters. – Dave X Aug 11 '14 at 4:52