How to create raster layer directly from numpy array in QGIS, without saving to file?

Question

I have managed to create a numpy array from a raster file using commands from the gdal library, and now I would like to load that raster as a layer without saving it first. It seems like there should be a simple way to do this. Is there?

Answer 1

I found some kind of a solution to create a temporary layer but it involves QGIS Temporary Output, so it writes to a temporary file anyway. I don't think it's possible to create a raster layer in memory with PyQGIS api.

Anyway, here's some code which uses QgsRasterDataProvider and QgsRasterBlock to manually set raster values from a numpy array:

import numpy as np

#creating a constant raster which will be overwritten later
extent = QgsRectangle()
extent.setXMinimum(457552.5)
extent.setYMinimum(248698.5)
extent.setXMaximum(459804.5)
extent.setYMaximum(251031.5)
crs = QgsCoordinateReferenceSystem('EPSG:2180')
params = {
    'EXTENT': extent,
    'TARGET_CRS': crs,
    'PIXEL_SIZE': 10,
    'NUMBER': 0.5,
    'OUTPUT_TYPE': Qgis.DataType.Byte,
    'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
}
r = processing.run('qgis:createconstantrasterlayer', params)['OUTPUT']
rlayer = QgsRasterLayer(r, 'temp', 'gdal')
provider = rlayer.dataProvider()

w = rlayer.width()
h = rlayer.height()

# dataType = Qgis.DataType.Byte
dataType = provider.dataType(1)
block = QgsRasterBlock(dataType, w, h)

#this is where you would use your numpy array to set value of each and every pixel
#arr is a random 2-dimensional array with pixel values
arr = np.random.rand(w, h)
for i in range(0, w):
    for j in range(0, h):
        block.setValue(i, j, int(arr[i][j] * 255))

# alternative using setData
# f = lambda x: int(x * 255)
# data = bytearray(np.array(map(f, arr)))
# block.setData(data)

provider.setEditable(True)
provider.writeBlock(block, band=1)
provider.setEditable(False)
provider.reload()

#this is needed for the min and max value to refresh in the layer panel
stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent)
min = stats.minimumValue
max = stats.maximumValue

renderer = rlayer.renderer()
myType = renderer.dataType(1)
GrayRenderer = QgsSingleBandGrayRenderer(provider,1)
contrastEnhancement = QgsContrastEnhancement.StretchToMinimumMaximum
myEnhancement = QgsContrastEnhancement()
myEnhancement.setContrastEnhancementAlgorithm(contrastEnhancement, True)
myEnhancement.setMinimumValue(min)
myEnhancement.setMaximumValue(max)
rlayer.setRenderer(GrayRenderer)
rlayer.renderer().setContrastEnhancement(myEnhancement)

QgsProject.instance().addMapLayer(rlayer)

Answer 2

This was not possible with QGIS2, but we will be able to edit raster in QGIS3 : https://qgis.org/api/classQgsRasterBlock.html#a69e7601b2c9b997805dd523438395b3f