Writing numpy array to raster file

I'm new to GIS.

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?

• Which slide on the tutorial are you referring to?
– R.K.
Oct 22 '12 at 11:22
• I love Python, but see soo many answers like these below that would be one liners in R. May 30 '20 at 12:32
• Has anyone tested Create method here with AAIGrid option? OR even AIG. It does not create the output_raster object. No error is thrown either. Metadata looks fine to create an Create object but I always get NoneType. GTiff works just fine.. Here is the code snippet I tried: ``` xmin,ymin,xmax,ymax = [lons.min().values,lats.min().values,\ lons.max().values,lats.max().values] ncols =lons.shape[0] nrows = lats.shape[0] xres = (xmax-xmin)/float(ncols) yres = (ymax-ymin)/float(nrows) geotransform=[xmin,xres,0,ymin,0,yres ] driver= gdal.GetDriverByName('AAIGrid') metadata = driver.GetMetadata() outp Sep 9 '20 at 3:45
• @geotheory, but then you'd have to write it in R, and everyone knows: R<-[i for i not in tuitive] Apr 8 at 1:09

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! Aug 7 '15 at 16:46
• I had to add `outDs = None` to get it saved Nov 20 '17 at 11:44

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

output_raster.FlushCache()
``````
• 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. Aug 11 '14 at 4:52
• This post is really useful, thanks. In my case, however, I am getting a tif file that is completely black when I open it as an image outside ArcGIS. My spatial reference is the British National Grid (EPSG=27700), and the units are meters. Mar 16 '17 at 9:13
• I have posted a question here: gis.stackexchange.com/questions/232301/… Mar 16 '17 at 9:39
• Did you find out how to set IAU2000:49900 Mars encoding? Mar 18 '19 at 7:38
• I am trying to use this solution. I successfully create a geotiff, but it is projected in the wrong place. I think it's because I'm using a south polar projection? Currently EPSG:3976 WGS 84 / NSIDC Sea Ice Polar Stereographic South. Can anyone help advise how I would do differently in this case. Mar 23 at 22:56

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. Oct 22 '12 at 19:03
• FYI, I didn't have the hanging issue with QGIS, I have version 1.8.0 too. Oct 22 '12 at 19:10
• Hi, Is there any link where i can find how to convert my numpy array to ASCII, so that final file i get geotiff to plot. Nov 2 '20 at 10:54

There is also a nice solution in the official GDAL/OGR Cookbook for Python.

This recipe creates a raster from an array

``````import gdal, ogr, os, osr
import numpy as np

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

cols = array.shape[1]
rows = array.shape[0]
originX = rasterOrigin[0]
originY = rasterOrigin[1]

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()

def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
reversed_arr = array[::-1] # reverse array so the tif looks like the array
array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster

if __name__ == "__main__":
rasterOrigin = (-123.25745,45.43013)
pixelWidth = 10
pixelHeight = 10
newRasterfn = 'test.tif'
array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
[ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
[ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
[ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
[ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
``````
• This recipe is good but there is an issue with the final tiff file. The lat-lon values of pixels are not proper. Aug 16 '19 at 13:50
• You might be seeing weird incompatibilities between ESRI WKT and OGC WKT: gis.stackexchange.com/questions/129764/… Aug 17 '19 at 17:13
• One thing I encountered is, the way mentioned by you will definitely change the array to a raster with ease. But we need to georeference this raster with the top left and bottom right coordinated using gdal_translate. One way to do that is following the two steps : 1) First find the top-left and bottom right lat-lon values through gdalinfo 2) Then, through gdal_translate utilise the geotiff (generated with the above mentioned approach of converting array to raster) to georeference it with the top-left and bottom right lat-lon coordinates. Aug 18 '19 at 7:59

An alternative to the approach suggested in the other answers is to use the `rasterio` package. I had issues generating these using `gdal` and found this site to be useful.

Assuming you have another tif file(`other_file.tif`) and a numpy array (`numpy_array`) that has the same resolution and extent as this file, this is the approach that worked for me:

``````import rasterio as rio

with rio.open('other_file.tif') as src: