10

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?

2
  • Possible duplicate:gis.stackexchange.com/questions/34082/…
    – Fran Raga
    Commented Jul 13, 2016 at 12:53
  • 1
    I saw this answer, thanks. I was hoping there was a way to do it without saving to disk. If not, that's OK I guess, but it seems like a wasted step in many instances.
    – Matt D.
    Commented Jul 13, 2016 at 15:17

2 Answers 2

2

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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