Tuesday 24 December 2019

python - 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?



Answer



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: raster version of above


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.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...