recording pixel value with Python script in ArcGIS 10.x

[DeanPodolsky]

A scripting newbie question...

I'd like to write a Python script for ArcGIS 10.x that would record the pixel value for all pixels in a raster. Furthermore, I'd like to only record a given value once in the results. For example, if there are 23 values of 180 overall, only record one instance of that value.

The raster in question is 4-band so the results would have to include values for all pixels in each band, recorded separately by band. The Pixel Inspector tool (see image below) will show you the value per pixel (and surrounding pixels), but I'd like to record the values as described.

Any help and examples would be greatly appreciated. Thanks in advance.


[ATTACH=CONFIG]30411[/ATTACH]

[NeilAyres]

There is a GetCellValue function that can be used in python to query values in a raster, but it is quite slow.
And there maybe really fast ways of doing this using numpy arrays.
But why not use RasterToPoint, then process the resulting point attribute table. Can't remember if R2P will export all 4 bands or if you would have to do each band separately.
If the raster is very large, it will of course give you a very large number of points to process.
Cheers and good luck,
Neil

[DeanPodolsky]

Thanks for your reply Neil.

I found this on another forum that goes toward what I'm trying to do:

import arcpy
raster = arcpy.Raster("image.tif")
array = arcpy.RasterToNumPyArray(raster) (height, width)=array.shape 
    for row in range(0,height): 
        for col in range(0,width): 
            print str(row)+","+str(col)+":"+str(array.item(row,col))

It only does a single band and prints the result to the screen but I'm trying to address both of those at the moment. The raster is 6017 x 7715 pixels so I'd like to remove any intermediate steps if possible.

Any thoughts or comments?

[PhilMorefield]

It appears that your rasters hold integer values. If so, you could just build and export the raster attribute tables. That will always give you a list of unique values.

[DeanPodolsky]

Thanks for your reply Phil.

Any thoughts on how best to do that? As I mentioned I'm new to scripting and have been researching a lot of different options. Lately I've been looking at rec2csv and numpy.savetxt as well as NumPyArrayToTable.

I'm getting closer to the goal but it's slow progress!

[ChrisSnyder]

Do you want a listing of the unique pixel values or to actually store each pixel's value?

[DeanPodolsky]

Ideally, to lessen processing, I'd like the unique values per row but I'd be happy with each pixel's value per row (and can sort it out in Excel).

[ChrisSnyder]

Not sure I understand your purpose here, but the most efficient way to store a pixel value from a raster is.... the original raster.

What information/analysis are you wanting to do? Knowing that will dictate the best way to proceed.

That said, I noticed that v10.2+ now supports multiband rasters in the RasterToNumpy tool: http://resources.arcgis.com/en/help/main/10.2/index.html#//018v00000023000000

Otherwise you would have to cobble several 2D arrays together yourself, which I believe there are methods to do this. Numpy is probably your best bet here I think, but I have little/no expertise in it's use.

[DeanPodolsky]

What I'm trying to do is capture the value of each pixel in a raster, but I'm only interested in each unique value per row.

In the 10 x 10 example below, the first row pixel values that would be listed are:
191,192,188,190,181,203,183,192,203,195

where the second row would have 8 unique values because 181 and 172 are repeated:
181,187,172,186,182,189,161,208

[ATTACH=CONFIG]30468[/ATTACH]

I suppose ideally the output would be like this:

row:values

0:191,192,188,190,181,203,183,192,203,195
1:181,187,172,186,182,189,161,208

I'm stuck with 10.1 at the moment so I'll have to work within those parameters.

[DouglasSands]

What I'm trying to do is capture the value of each pixel in a raster, but I'm only interested in each unique value per row.

In the 10 x 10 example below, the first row pixel values that would be listed are:
191,192,188,190,181,203,183,192,203,195

where the second row would have 8 unique values because 181 and 172 are repeated:
181,187,172,186,182,189,161,208

[ATTACH=CONFIG]30468[/ATTACH]

I suppose ideally the output would be like this:

row:values

0:191,192,188,190,181,203,183,192,203,195
1:181,187,172,186,182,189,161,208

I'm stuck with 10.1 at the moment so I'll have to work within those parameters.

There are two straightforward ways you could do this, I have to do something kind of similar. The fastest way would be to convert the raster into a NumPy array using arcpy.RasterToNumpyArray, and then iterate over the rows in the array and use set(). Also using nparray.tolist() tends to speed things up in my experience. Code might look something like this:

import arcpy
import numpy

raster = 'path\to\raster'
output = 'the\output.file'

data = arcpy.RasterToNumPyArray(raster).tolist
result = open(output, 'w')

rowNum = 0
for row in data:
    unique = set(row)
    result.write('%s:'%rowNum)
    line = ''
    for val in unique:
        line += '%s,'%val
    line = line[:-1] + '\n' #eliminates a trailing comma
    result.write(line)
    rowNum += 1
    
result.close()

I haven't tested that, but something like that should be what you need. However this will fail to work if your raster is to large to put the whole thing into memory at the same time. In this case, you could export the file to an ASCII raster first, and then do the following:

import arcpy

raster = 'path\to\ASCIIraster.asc'
output = 'the\output.file'

ASCIIData = open(raster, 'r')
OutData = open(output, 'w')

rowNum = 0
for line in ASCIIData:
    values = line.split()
    writeLine = '%s:'%rowNum
    try:
        test = int(values[0]) #This will fail for the headers
        unique = set(values) #All strings already
        for val in unique:
            writeLine += val + ','
        writeLine = writeLine[:-1] + '\n' #eliminates a trailing comma
        OutData.write(writeLine)
        rowNum += 1 #Has to be here because we don't increment for the header
    except ValueError:
        pass #Pass because we are skipping the text header

OutData.close()

This will be slower, but for small datasets the difference will be minimal and it will work for datasets that will fail with NumPy.

Hope this helps.
-Doug