It seems that gdal can caused in some crashes in QGIS, specific in used for raster layers (such as links : QGIS crashes on close if data was read through GDAL in Python Console or GDAL/QGIS broken Ubuntu 14.04 or Why close QGIS for extract pixel values by GDAL?. So Is it possible to get pixel values only with pyqgis and without gdal?
-
1is docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/… what you're after?– Michael StimsonJul 12, 2017 at 21:49
-
1When you have it working can you answer your own question with some of your code, I'm sure that other users might have a similar problem and would like to see how it's done.– Michael StimsonJul 12, 2017 at 21:56
3 Answers
Yes, it is possible to get pixel values without gdal; only PyQGIS. Easier way is with a QgsRasterBlock object. This is the code:
layer = iface.activeLayer() #in my case, a 20x20 raster
provider = layer.dataProvider()
extent = provider.extent()
rows = layer.height()
cols = layer.width()
block = provider.block(1, extent, cols, rows)
for i in range(rows):
for j in range(cols):
print block.value(i,j),
print
After running the code at Python Console of QGIS, I got:
27.0 17.0 48.0 60.0 32.0 53.0 73.0 5.0 94.0 35.0 18.0 31.0 35.0 63.0 54.0 26.0 79.0 17.0 3.0 98.0
32.0 44.0 7.0 9.0 54.0 52.0 12.0 19.0 75.0 26.0 27.0 28.0 36.0 58.0 25.0 83.0 83.0 67.0 70.0 6.0
27.0 17.0 44.0 63.0 82.0 76.0 64.0 69.0 42.0 52.0 81.0 4.0 46.0 37.0 97.0 86.0 88.0 44.0 17.0 16.0
36.0 70.0 53.0 91.0 34.0 9.0 90.0 47.0 39.0 23.0 30.0 27.0 36.0 89.0 99.0 48.0 23.0 55.0 56.0 77.0
75.0 96.0 75.0 64.0 36.0 72.0 55.0 87.0 35.0 56.0 23.0 6.0 97.0 58.0 97.0 54.0 14.0 93.0 9.0 54.0
41.0 58.0 52.0 3.0 96.0 45.0 32.0 85.0 32.0 98.0 9.0 26.0 63.0 20.0 83.0 10.0 67.0 24.0 58.0 19.0
61.0 5.0 15.0 17.0 53.0 57.0 57.0 59.0 79.0 43.0 7.0 65.0 83.0 65.0 32.0 59.0 22.0 24.0 25.0 22.0
66.0 60.0 3.0 94.0 53.0 93.0 10.0 8.0 24.0 16.0 92.0 95.0 79.0 42.0 97.0 45.0 82.0 87.0 31.0 8.0
2.0 40.0 34.0 5.0 6.0 37.0 9.0 83.0 28.0 11.0 54.0 84.0 35.0 30.0 46.0 24.0 49.0 15.0 30.0 36.0
30.0 45.0 65.0 64.0 18.0 49.0 78.0 70.0 83.0 64.0 48.0 73.0 69.0 2.0 69.0 1.0 30.0 20.0 45.0 74.0
7.0 93.0 76.0 96.0 60.0 84.0 7.0 81.0 79.0 6.0 47.0 80.0 9.0 61.0 43.0 19.0 53.0 100.0 28.0 49.0
77.0 4.0 79.0 84.0 65.0 14.0 91.0 62.0 41.0 1.0 78.0 26.0 100.0 66.0 37.0 8.0 96.0 51.0 83.0 94.0
100.0 28.0 77.0 7.0 3.0 72.0 82.0 97.0 3.0 25.0 64.0 66.0 25.0 75.0 54.0 29.0 15.0 36.0 25.0 57.0
7.0 63.0 64.0 46.0 41.0 13.0 6.0 34.0 55.0 32.0 99.0 87.0 12.0 100.0 9.0 31.0 48.0 22.0 7.0 52.0
45.0 43.0 9.0 17.0 39.0 98.0 20.0 51.0 80.0 91.0 52.0 18.0 40.0 39.0 5.0 97.0 74.0 5.0 9.0 99.0
67.0 53.0 24.0 67.0 32.0 17.0 35.0 98.0 37.0 68.0 18.0 31.0 71.0 95.0 59.0 100.0 97.0 60.0 22.0 6.0
97.0 73.0 69.0 93.0 13.0 87.0 29.0 69.0 10.0 62.0 35.0 26.0 62.0 51.0 30.0 15.0 25.0 62.0 1.0 61.0
59.0 53.0 22.0 74.0 5.0 52.0 13.0 29.0 14.0 18.0 9.0 99.0 70.0 3.0 53.0 100.0 39.0 20.0 73.0 30.0
66.0 72.0 92.0 26.0 77.0 67.0 74.0 60.0 55.0 68.0 71.0 61.0 32.0 72.0 100.0 49.0 80.0 74.0 8.0 1.0
53.0 4.0 82.0 72.0 86.0 82.0 92.0 22.0 52.0 8.0 31.0 30.0 6.0 2.0 14.0 73.0 34.0 60.0 84.0 96.0
I corroborated, with help of Value Tool plugin, that values extracted were correct.
Based on @Michael Stimson in above and QGIS document (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/raster.html#query-values). In below example it is possible to get pixel values without gdal.
from qgis.core import *
def pixel2coord(x, y):
xp = (pixelWidth * x) + originX + (pixelWidth/2)
yp = (pixelHeight * y) + originY + (pixelHeight /2)
return QgsPoint(xp, yp)
# Open tif file
ds = QgsRasterLayer("/Data/Sample_Polygon.tif","Study")
pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()
pntRstList = []
for i in range(0, src_cols):
for j in range(0, src_rows):
rspnt = pixel2coord(i,j)
pntRstList.append(rspnt)
for x, y in pntRstList:
ident = ds.dataProvider().identify(QgsPoint(x, y),
QgsRaster.IdentifyFormatValue)
if ident.isValid():
print ident.results()
-
1Yes, it is possible to get pixel values without gdal, only PyQGIS, but easier way is with a QgsRasterBlock object.– xunilkJul 13, 2017 at 0:47
the first solution is a really good one! I just would like to propose instead of print or instead of put it in a list like in the second solution, put the values in numpy array. If you are planning to work with large maps is better to take this approach once the numpy array is much faster.
from qgis.core import *
import numpy as np
def GetValuePixels(layer):
provider = layer.dataProvider()
extent = provider.extent()
rows = layer.height()
cols = layer.width()
block = provider.block(1, extent, cols, rows)
values = np.zeros(shape=(rows,cols)) #Creates space in the memory to avoid spend time...
for i in range(rows):
for j in range(cols):
values[i,j] = block.value(i,j)
return values
thanks!