Converting QgsRaster object to array

Question

So I'm a newbie, so sorry if this is an idiot question.

I'm trying to write a plugin for QGIS3 for my organisation, and I'd like to output a standalone figure (i.e. not in the QGIS window) as the final output. This way I can show rasters, line plots and vectors, in multiple subplots in a single figure and the user can save the figure directly if they want to.

I'm using Matplotlib as that's already available in pyQGIS, so I want to stick with that for logistic reasons.

My problem is that I have some nice QgsRasterLayer rasters (and QgsVectorLayer vectors too) that I'd like to display as some of the subplots but I don't know how to convert them into a form suitable display in Matplotlib. I can display them into the main QGIS window but that's not what I'm after. I need to plot them with Matplotlib.

Can anyone give me some suggestions? It's probably really obvious, but I'm just not finding the solution as my python experience is limited.

Answer 1

Barry Rowlingson build an plugin for QGIS2 called "Rasterlang" (https://www.maths.lancs.ac.uk/~rowlings/Software/Qgis/Plugins/rasterlang/), inside there is a modul named layers.py with a function called layerAsArray():

def layerAsArray(layer):
    """ read the data from a single-band layer into a numpy/Numeric array.
    Only works for gdal layers!
    """

    gd = gdal.Open(str(layer.source()))
    array = gd.ReadAsArray()
    return array

Like you I was concerned about loading data twice, but didn't see a way through QgsRasterLayer. Provider seems to be always gdal, but you can't "cast" it back and use ReadAsArray() on the provider.

Answer 2

I'm facing the same dilemma, you want to create statistical graphs of the image.
You can create a histogram using PyQGis:

1. layer= iface.activeLayer()

2. provider=layer.dataProvider()

   You can ask the provider to generate the histogram, specifying the band, the number of cuts and if you want the extension (to create the histogram of a portion of the image)

3. histogram = provider.histogram(1,30)

4. y=histogram.histogramVector

   The histogramVector is the frequency of the defined cuts (30 line 3)/ How to obtain the x-axis values? You can request the minimum and maximum values of the histogram with histogram.minimum and histogram.maximum

Then use numpy lynspace to generate them

5. x=np.linspace(histogram.minimum, histogram.maximum,30,endpoint=False)

Then graphs with matplotlib bar

For the other graphics I think the only thing left to do is to use GDAL as an array, it's faster. Another way to read the data and save it in an array would be using QgsRasterBlock but it doesn't make sense, it's probably slower than reading the image again with GDAL.