Raster are matrices with geospatial reference. In order to make a finite difference model or any other model compatible with QGIS we have to process the matrix with its parameters to be represented as a spatially referenced raster. The process involves some Python scripting on the QGIS console, processing time is low and the process can be applied for many matrices without visual representation of intermediate files that saves computational time and resources.

This tutorial show the complete process of matrix construction, parameter calculation and spatially representation in the Python console of QGIS.

Tutorial

Code

This is the code that generates the raster and represent it as geotiff:

import os
import numpy as np
from osgeo import gdal, gdal_array, osr

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)))
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)))

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]

ruta = "C:\Users\Saul\Documents\Ih_HowtoConvertaMatrixinaGeospatialTiffwithPyQGIS"
raster = gdal.GetDriverByName('GTiff').Create(os.path.join(ruta,"myraster.tif"),ncols,nrows,1,gdal.GDT_Float32)
raster.SetGeoTransform(geotransform)

srs=osr.SpatialReference()
srs.ImportFromEPSG(4326)

raster.SetProjection(srs.ExportToWkt())
raster.GetRasterBand(1).WriteArray(array)

iface.addRasterLayer(os.path.join(ruta,"myraster.tif"),"MyRaster","gdal")