4

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.

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.

2
  • Is it possible for you to use gdal or rasterio? I found some answers here: gis.stackexchange.com/questions/32995/…. This exports the data to a numpy array, which matplotlib understands. Commented Apr 3, 2019 at 8:57
  • Thanks - yes GDAL is possible, but I need the rasters in a QgsRasterLayer and display them in the main QGIS interface. That part is all ok. But if I took a GDAL approach I assume I'd have to open them a second time, which seems inelegant. Or have I got that wrong? I suspect the answer lies somewhere in QgsRasterLayer.dataProvider.block but so far it's still escaping me.
    – SJJ
    Commented Apr 3, 2019 at 13:03

3 Answers 3

3

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.

1
  • Ok thanks - am re-assured that others have struggled with the same problem. It's not just me then. As some of my rasters are large I am concerned about the time and memory implications of loading twice. But if its the only route, then it's the one I will take although I remain open to other suggestions. Thank you for your kind help. Sarj
    – SJJ
    Commented Apr 3, 2019 at 15:32
0

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

  1. 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.

0
from qgis.core import QgsProject    
import numpy as np    
import matplotlib.pyplot as plt    
rlayer = QgsProject.instance().mapLayersByName('Ta_')[0]    
provider = rlayer.dataProvider()    
extent = provider.extent()    
rows, cols = rlayer.height(), rlayer.width()    
block = provider.block(1, extent, cols, rows)  # !!!!!  
z = np.zeros((rows, cols))    
for i in range(rows):    
    for j in range(cols):    
          z[i,j] = block.value(i,j)    
z[z == provider.sourceNoDataValue(1)] = np.nan    
x = np.linspace(extent.xMinimum(), extent.xMaximum(), cols)    
y = np.linspace(extent.yMaximum(), extent.yMinimum(), rows)    
plt.close('all')    
plt.pcolor(x, y, z, cmap = 'jet')    
plt.show()
1
  • Kindly add some context how this code will help in solving the problem.
    – Padmanabha
    Commented Oct 11, 2023 at 16:49

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.