Finding coordinates of given raster value? [closed]

Question

I would like to get the coordinates of a raster cell with a given value. I suppose many cells can take the same value, in which case several coordinates should be returned. In my first problem I'm interested in the location of the maximum value(s), but also seeking other values in planned future analysis.

I will work with very small rasters, which I cut based on polygons before processing, so speed is not so important in this case. Although fast solutions are welcome, also for wider applicability for others.

I do batch processing using python and QGIS, and can of course use other tools such as SAGA, GRASS, numpy, scipy to do this. The problem seems simple, but I have been completely stuck at this for a while now. I'm quite new to using python and QGIS.

Answer 1

Okay, here's a way to do it using Numpy and Rasterio. I'll show you an excerpt from Rasterio's rio-insp prompt using a test raster (rio insp tests/data/RGB.byte.tif).

>>> r, g, b = src.read()
>>> rc = np.transpose(np.nonzero(r==42))
>>> xy = [src.ul(row, col) for row, col in rc]
>>> print xy[0:5]
[(308411.0935524652, 2797210.8635097495), (331213.9759797724, 2793610.3621169915), (303910.5246523388, 2793010.278551532), (330913.93805309734, 2791510.069637883), (331213.9759797724, 2791510.069637883)]

I've read a RGB GeoTIFF raster into three arrays. The expression r==42 gives you a new array which is True where the value of r is 42 and False elsewhere. The transpose of the array returned by np.nonzero(r==42) contains a 1D array of (row, col) pairs. Those are the pixels of r where the value is 42. To get their georeferenced coordinates, map the .ul() method of the Rasterio dataset over that array. In the case above you get easting, northing pairs (upper left corner of your selected pixels) in EPSG:32618.

Answer 2

With GRASS, you can use a combination of r.mapcalc and r.out.xyz.

First, you use r.mapcalc to select only those cells with the desired value: something like intermediate_map=if(source_map==value) should create a map that contains only the cells with the desired value.

Second, you export this intermediate map using r.out.xyz input=intermediate_map to create a list of the centerpoints of the cells. A more detailed output can be achieved with r.stats.

Answer 3

I tried to use rasterio and it didn't work with sgillies'example. I don't know why. For this reason, I used next code (gdal/struct), with a very simple raster (test.tif: 20x20, values between 1 and 50), to verify correct operation for value 15.

from osgeo import gdal
import struct

nameraster = "test.tif"

dataset = gdal.Open(nameraster)
geotransform = dataset.GetGeoTransform()
band = dataset.GetRasterBand(1)

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 
            'Int32':'i', 'Float32':'f', 'Float64':'d'}

print "rows = %d columns = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)

print "Data type = ", BandType

print "Executing with %s" % nameraster

print "test_value = 15"

X = geotransform[0] #x coordinate
Y = geotransform[3] #y coordinate

for y in range(band.YSize):

    scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)

    for value in values:

        if(value == 15):        
            print "%.4f %.4f %.2f" % (X, Y, value)
        X += geotransform[1] #x pixel size
    X = geotransform[0]
    Y += geotransform[5] #y pixel size

dataset = None

The result obtained in bash console (GNU/Linux Debian), visually verified in QGIS, is:

rows = 20 columns = 20
Data type =  Float64
Executing with test.tif
test_value = 15
609891.6076 995435.4677 15.00
609769.1654 995221.1939 15.00
609677.3338 995159.9728 15.00
609891.6076 995159.9728 15.00
609952.8287 995159.9728 15.00
610014.0498 995068.1412 15.00
609616.1127 995037.5306 15.00