3

The code below is a clunky attempt to get a raster numpy array loaded up for basic machine learning from the QGIS Python console. It usually works but has stability issues with the GetRasterBand() and I have asked a separate question. As an improvement I am trying to figure how to get the array values directly from the loaded rasters in QGIS. I.e. you open the raster layers in QGIS, do whatever is needed (e.g. training polygons) and then run the code.

#Part 1. Importing Python modules
import numpy as np
import gdal
from sklearn.cluster import KMeans
import os



#Part 2. Specify the image file to be classified and the class raster with the training data.

MyFolder = 'D:\\RemoteSensing\\Practical_7\\' #Change this to match the folder on your drive.  Note that we use \\ instead of \
MyImageRoot = 'S2_mini' #Specify the root image name. Do not nclude the extension, e.g. 'glacier_'.  If a color image, just write the file name without extension.
Nclasses = 4 #how many classes are in the image
MyClassRaster = 'KmeansTEST_4cls'
BandList = ['B02'] #Indicate the bands you want to load.  B


#Part 3 Load the image 
fileName = os.path.join(MyFolder + MyImageRoot + BandList[0] + '.tif')
ImageFile = gdal.Open(fileName)
band1 = ImageFile.GetRasterBand(1)
array1 = band1.ReadAsArray()

I am trying to get array1 by direct querry to a loaded band in the QGIS layers. I have found QgisProject:

layer = QgsProject.instance().mapLayersByName(BandList[0])

But I can't find how to access the raster values for the layer and have them as a numpy array for the next processing steps. I'm working with QGIS 3.2 and 3.4.

3

For accessing raster values for the layer in PyQGIS you can use QgsRasterBlock objects. It is exemplified in the following code for a raster named 'aleatorio':

layer = QgsProject.instance().mapLayersByName('aleatorio')

provider = layer[0].dataProvider()

extent = provider.extent()

rows = layer[0].height()
cols = layer[0].width()

# 1 is referring to band 1
block = provider.block(1, extent, cols, rows)

#accessing one particular value; e.g. value[4,3]
print (block.value(4,3))

After running above code, its value for (4,3) indexes is printed at Python Console (equal to 64). With Value Tool plugin it can be corroborated (following image).

enter image description here

  • After you have the raster block you can call QgsRasterBlock.data() to access the underlying pixel data in an efficient manner. It's possible to directly create a numpy array from this value. – ndawson Nov 17 '18 at 20:58

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.