2

This is really a question of how to better understand QGIS API documentation (C++ or Python) but I'll use this example:

How do I use QgsRasterInterface.histogram? I was told in this thread I can use QgsRasterInterface.histogram to calculate a list of unique values and their counts in a raster layer.

My best attempt so far:

layer = iface.activeLayer()        #the active layer is a raster layer
provider = layer.dataProvider()

print(provider.hasHistogram(1, 127))   #returns false, no histogram yet

provider.histogram(1)    #for band 1

and QGIS crashes. No errors or other information, just a crash. I see that QgsRasterInterface has a method .initHistogram. Maybe I need to initialize the histogram first?:

enter image description here

The first argument of .initHistogram is a QgsRasterHistogram. I am unsure how to make this QgsRasterHistogram object. The documentation for this class lists no methods. Though it does have some attributes which are inherited from QgsRasterInterface, and none of the attributes have an description. The QgsRasterHistogram has a short description which makes me think I am on the right path:

enter image description here

By doing a quick test I can see I have not constructed the QgsRasterHistogram properly.

histogram = QgsRasterHistogram()
print(histogram.extent)     #prints empty extent

The documentation is mostly empty for QgsRasterHistogram and when I try to use it in the PyQGIS console, I am not prompted to add any arguments or anything.

Do I need to use .initHistogram to construct a histogram for this raster layer? Or is there another reason QGIS crashes when I use provider.histogram?

6

If you are executing your Python code inside the QGIS Python console, there is already an existing object called iface that you can use.

For stand-alone scripts, see Why does iface.activeLayer() not work in a PyQGis standalone app?

l = iface.activeLayer()
p = l.dataProvider()
p.initHistogram(QgsRasterHistogram(),1,100)
h = p.histogram(1)

returns

<qgis._core.QgsRasterHistogram object at 0x0000016D963189D8>

Afterwards accessing some histogram values, such as h.minimum and h.maximum works but when running

h.histogramVector

QGIS crashes, so I've created a ticket: https://issues.qgis.org/issues/21885 but it's well possible that there's something missing from the Python code here ...

| improve this answer | |
  • Thank you for the example. One more question, what does histogramVector do? There is no information in the documentation. – Erich Purpur Apr 19 '19 at 13:45
  • 1
    @ErichPurpur based on Nyall's answer to your previous question, I expect that some property of the histogram has to give you access to the computed distribution and histogramVector seems to be the most likely candidate – underdark Apr 19 '19 at 15:31
  • h.histogramVector works for me (3.9 nightly, Windows 10, using IPyConsole plugin from git). It returns a list of length h.bins containing count integers. – foobarbecue Oct 9 '19 at 22:28
  • 2
    Ah, it was fixed. OP's bug report got ported to here: github.com/qgis/QGIS/issues/29700 – foobarbecue Oct 9 '19 at 22:32
0

I'm adding an alternate solution to my question. This is a workaround to the .initHistogram() method as it seems to be broken. It is actually a solution I adapted from another StackOverflow thread of mine. Here is the original thread: PyQGIS Raster Band Stats - Access pixel values of raster layer

Basically what I am trying to do is to re-create the properties > symbology tab, but programmatically. Using the symbology tab allows you to see the color, raster values, and a histogram of those raster values for a given raster layer. Here is an example of what it looks like: enter image description here

The block() method from QgsRasterBandStats class allows you to work with the individual pixels in a raster layer, creating a QgsRasterBlock object. Here is a code example of reading a raster layer, using the block() method, and accessing the pixels. Also, I'm using a really tiny clip of a raster layer, for the sake of this example. If you are iterating over a raster with many thousands or millions of pixels, this could take a while:

layer = iface.activeLayer()
provider = layer.dataProvider()
stats = provider.bandStatistics(1, QgsRasterBandStats.All) 

extent = layer.extent()
rows = layer.height()
columns = layer.width()
block = provider.block(1, extent, columns, rows)

values = []
for row in range(rows):
    values.append([])
print(f"Values: {values}")

for row in range(rows):
    for column in range(columns):
        values[row].append(block.value(row, column))
print(f"Values: {values}")

flattened_list = []
for list in values:
    for element in list:
        flattened_list.append(element)

unique_values = []
for item in flattened_list:
    if item not in unique_values:
        unique_values.append(item)
print(f"Unique Values: {unique_values}")

for item in unique_values:
    count = 0
    for element in flattened_list:
        if element == item:
            count += 1
    print(f"Value: {item}, Count: {count}")

Here is an example of what the output should look like: enter image description here

| improve this answer | |

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.