10

I am working on a plugin for Qgis to calculate spatial Kernel density maps. I have all the calculations working, all I am missing is a way to turn a Numpy Array, with density values into a multiband raster layer.

Do I have to create a geotiff on a temp file using Gdal and then load it?

Or is there a direct way to create the layer from data in memory?

if so, how to do it?

6

Here is the code that I use to convert an array to gdal raster saving it to the disk, "param" is a dicitionary containing gdal parameters (check the gdal documentation) and "array" is a numpy array. Than you can instantiate a QgsMapLayer with your file as source. You have to create the geotiff in the disk.

    from osgeo import gdal as osgdal  # Adapt the import to fit yor environement.

    driver = osgdal.GetDriverByName(param['out_format'])

    dataset = driver.Create(
            param['dst_filename'],
            param['x_pixels'],
            param['y_pixels'],
            1,
            osgdal.GDT_Float32,
            )

    dataset.SetGeoTransform((
            param['xmin'],           #0
            param['pixel_size'],     #1
            0,                       #2
            param['ymin'],           #3
            0,                       #4
            param['pixel_size']))    #5

    out_srs = osr.SpatialReference()
    out_srs.ImportFromEPSG(param['SRID'])

    dataset.SetProjection(out_srs.ExportToWkt())
    dataset.GetRasterBand(1).WriteArray(array.T)  # Remove "T" if it's inverted.
    dataset = None
| improve this answer | |
  • Thanks, what is pixel size? (xmax-xmin)/x_pixels? – fccoelho Sep 25 '12 at 11:07
  • what are the 0s in the SetGeoTransform? – fccoelho Sep 25 '12 at 11:22
  • For all the questions check here - SetGeoTransform: gdal.org/… – Pablo Sep 25 '12 at 13:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.