[GIS] How to extract X Y coordinates and RGB from pixels using python in QGIS

coordinatespythonqgisrasterrgb

I'm trying to white a Python plugin in QGIS (2.2) to get, for each pixel of a raster image, its X Y coordinates, as well as, its RGB (3 bands) values.

At first I opened the raster file using: rlayer = QgsRasterLayer(rfileName, rbaseName)

Now I don't know how to get, for example, for pixel (1,1) its coordinates (X,Y) and its RGB color values.

Anybody could help me?

I know that I'll need to implement a DO WHILE, but I don't know the commands to extract this information of each pixel.

Thanks in advance,
Mateus

Best Answer

As a geologist, I make geological cross section using the elevations values from a DEM and the colors from a geological map, look at the pure Python solution with osgeo.gdal in Python Script for getting elevation difference between two points

But now, since PyQGIS 2.x, it is easier with PyQGIS and the QgsRaster.IdentifyFormatValuefunction

example with a DEM (one band):

myDEM = qgis.utils.iface.activeLayer()  
print myDEM.bandCount()  
1   
print myDEM.dataProvider().identify(QgsPoint(229774,111171), QgsRaster.IdentifyFormatValue)  
{1: 221.0}

example with a classic raster (three bands -> R,G,B values):

 print myraster.dataProvider().identify(QgsPoint(229774,111171), QgsRaster.IdentifyFormatValue).results()
 {1: 180.0, 2: 142.0, 3: 125.0}

The result is a dictionary (key= band number) and you can create a simple function that returns the values:

def val_raster(point,raster):
    return raster.dataProvider().identify(point, QgsRaster.IdentifyFormatValue).results().values()

elevation = val_raster(QgsPoint(229774,111171),myDEM)
print elevation
221.0
R,G,B = val_raster(QgsPoint(229774,111171),myraster)
print R,G,B
180.0, 142.0, 125.0

I presented a complete solution with PyQGIS in French PyQGIS (QGIS 2): geological cross-sections (colorization of a topographic profile build from a DEM with the colors of a raster and placement of the intersection points of the geological layers boundaries )

enter image description here

Related Solutions

[GIS] How to get X Y coordinates and cell value of each pixel in a raster using Python

Following @Dango's idea I created and tested (on small rasters with the same extent and cell size) the following code:

import arcpy, numpy

inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"

##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)

#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin

# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width

##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
    coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
    return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord

def rasCentrY(rasHeight, rasWidth):
    coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
    return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord

#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))


##Raster conversion to NumPy Array
#create NumPy array from input rasters 
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)

#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)


# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))

Based on @hmfly code, you can have access to desired values:

(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
    for col in range(0,width):
        #now you have access to single array of values for one cell location

Unfortunately there's one 'but' - the code is right for NumPy arrays which can be handled by system memory. For my system (8GB), the largest array was about 9000,9000.

As my experience doesn't let me provide more help, you can consider some suggestions about dealing wiht large arrays: https://stackoverflow.com/questions/1053928/python-numpy-very-large-matrices

arcpy.RasterToNumPyArray method allows to specify the subset of raster converted to NumPy array (ArcGIS10 help page) what can be useful when chunking large dataset into submatrices.

[GIS] How to get feature values from multiple polygons of shapefile while using extract function on raster

You have to unlist() the output of extract, and maintain that list-level grouping to know which object the value is from - this is a pain, and can be tricky for non-matches, so I put the workflow into package 'tabularaster' https://CRAN.R-project.org/package=tabularaster

You might find this is enough, where r is a raster and poly is a polygondataframe:

library(tabularaster)
cells <- cellnumbers(r, poly)
cells$value <- extract(r, cells$cell_)

So value is the pixel value, and cell_ is the pixel index, and object_ is the row-number of poly.

Then use object_ to get a value from the query layer:

## replace "polyID" with whatever column you want
cells$polyID <- poly$polyID[cells$object_]

For a case like this, making the cell_ explicit is overkill, but it works well for general extraction (i.e. time series coming in the door).

Best Answer

Related Solutions

[GIS] How to get X Y coordinates and cell value of each pixel in a raster using Python

[GIS] How to get feature values from multiple polygons of shapefile while using extract function on raster

Related Question