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

I’m really new to Python and I would like to know whether there is a quick method to get cell values of a raster pixel by pixel and the coordinates (map X Y coordinate of centre of each pixel) using Python in ArcGIS 10? To describe this further, I need to get the map X, map Y and cell value of the first pixel and assign those three values in to three variables and repeat this step for rest of the other pixels (loop through the entire raster). Thanks.

I think I need to describe my question more. Problem is, I need to get the X Y location of a pixel of the first raster and get cell values of several other rasters corresponding to that X Y location. This process should be loop through every pixel of the first raster without creating any intermediate point shapefile as it is going to be really really time consuming since I have to handle a raster with nearly 8 billion pixels. Also, I need to do this usin Python in ArcGIS 10.

@JamesS: Thank you very much for your suggestion. Yes this would work for one raster but I need to collect the cell values for several other rasters as well. The problem is, after getting the X and Y coordinate of the first pixel of the first raster, I need to get the cell value of the second raster corresponding to that X, Y location of the first raster, then third raster and so on. So, I think when looping through the first raster, getting X and Y location of a pixel and getting cell values of the other raster corresponding to that location should be done simultaneously but I’m not sure. This can be done by converting the first raster in to a point shapefile and performing Extract multivalues to point function in ArcGIS 10 but I’m unable to follow that method because I don’t want to create a shapfile as it is going to be really slow and I HAVE TO FIND A SOLUTION USING PYTHON, not with any existing tool in ArcGIS.

@hmfly: Thanks, Yes this method (RastertoNumpyarray) will work if I can get the coordinate of a known row and column value of the array.

@whuber: I don’t want to perform any calculations, all I need to do is write X Y coordinates and cell values in to a text file and that’s all

-
Perhaps you just want to do some math on the whole raster? Raster calculators work pixel by pixel. –  BWill Feb 24 '12 at 10:50
please describe your purpose in more detail. –  BWill Feb 24 '12 at 10:50
Normally, efficient and reliable solutions are obtained by using Map Algebra operations rather than looping over points. Limitations in Spatial Analyst's map algebra implementation preclude this approach from working in every case, but in a surprisingly large number of situations you don't have to code a loop. What calculation do you need to perform, exactly? –  whuber Feb 26 '12 at 20:03
Re your edit: of course that's a legitimate purpose. The format may be imposed on you by the needs of software further down the pipeline. But considering that writing 8 billion (X,Y, value1, ..., value3) tuples will require between 224 billion bytes (in binary) and perhaps 400 billion bytes (in ASCII), either of which is a rather large dataset, it may be worth finding alternative approaches to whatever it is you're ultimately trying to accomplish! –  whuber Feb 27 '12 at 2:12

If you just want to get the pixel values through (row,column), you can write an arcpy script like this:

import arcpy
raster = arcpy.Raster("yourfilepath")
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))


But, if you want to get the pixel's coordinate, the NumPyArray can't help you. You could convert the raster to point by RasterToPoint Tool, and then you can get the coordinate by Shape filed.

-

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))

rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))


(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: http://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.

-
Marcin's code is super! thanks, but it doesn't write the X,Y of the raster with the same resolution of the raster I mean the x and y grow 1 m and not , for example) 100 meters.... Do you have a suggestion to fix that Thanks –  user6517 Mar 23 '12 at 9:02

One way to do this would be to use the Raster_To_Point tool followed by the Add_XY_Coordinates tool. You'll end up with a shapefile where each row in the attribute table represents a pixel from your raster with columns for X_Coord, Y_Coord and Cell_Value. You can then loop over this table using a cursor (or export it to something like Excel if you prefer).

If you only have one raster to process, it's probably not worth scripting - just use the tools from ArcToolbox. If you need to do this for many rasters, you could try something like this:

[Note: I don't have ArcGIS 10 and am not familiar with ArcPy, so this is just a very rough outline. It's untested and will almost certainly need tweaking to get it to work.]

import arcpy, os
from arcpy import env

# User input
ras_fold = r'path/to/my/data'           # The folder containing the rasters
out_fold = r'path/to/output/shapefiles' # The folder in which to create the shapefiles

# Set the workspace
env.workspace = ras_fold

# Get a list of raster datasets in the raster folder
raster_list = arcpy.ListRasters("*", "All")

# Loop over the rasters
for raster in raster_list:
# Get the name of the raster dataset without the file extension
dataset_name = os.path.splitext(raster)[0]

# Build a path for the output shapefile
shp_path = os.path.join(out_fold, '%s.shp' % dataset_name)

# Convert the raster to a point shapefile
arcpy.RasterToPoint_conversion(raster, shp_path, "VALUE")

# Add columns to the shapefile containing the X and Y co-ordinates


You can then loop over the shapefile attribute tables using a Search Cursor or (possibly simpler) using dbfpy. This will allow you to read the data from your raster (now stored in a shapefile .dbf table) into python variables.

from dbfpy import dbf

# Path to shapefile .dbf
dbf_path = r'path\to\my\dbf_file.dbf'

# Open the dbf file
db = dbf.Dbf(dbf_path)

# Loop over the records
for rec in db:
cell_no = rec['POINTID'] # Numbered from top left, running left to right along each row
cell_x = rec['POINT_X']
cell_y = rec['POINT_Y']
cell_val = rec['GRID_CODE']

# Print values
print cell_no, cell_x, cell_y, cell_val

-

The simplest method to output coordinates and cell values to a text file in ArcGIS 10 is the sample function, no need for code and especially no need to loop over each cell. In ArcGIS<=9.3x raster calculator it used to be as simple as outfile.csv = sample(someraster) which would output a text file of all (non null) cell values and coordinates (in z,x,y format). In ArcGIS 10, it looks like the "in_location_data" argument is now mandatory so you need to use the syntax Sample(someraster, someraster, outcsvfile).

Edit: You can also specify multiple rasters: Sample([someraster, anotherraster, etc], someraster, outcsvfile). Whether this would work on 8 billion cells, I have nooo idea...

Edit: Note, I haven't tested this in ArcGIS 10, but have used the sample function for years in <=9.3 (and Workstation).

Edit: I have now tested in ArcGIS 10 and it won't output to a text file. The tool changes the file extension to ".dbf" automatically. However... the following python code works as SOMA and MOMA map algebra statements are still supported in ArcGIS 10:

import arcgisscripting
gp=arcgisscripting.create()
gp.multioutputmapalgebra(r'%s=sample(%s)' % (outputcsv,inputraster))

-
Very nice. Thanks for pointing this out - I hadn't noticed this tool before. Certainly much neater and simpler than my solution! –  JamesS Mar 26 '12 at 9:06

Maybe you could create a world file for the raster, covert the raster to a numpy array. then if you loop over the array you will get the cell values and if you increamentaly update the x,y from the world file you will also have the coordinates for each cell value. hope that is usefull.

-
If you're not interested in the Raster to Point tool method suggested by JamesS, I'd say this is the way to go. –  nmpeterson Feb 24 '12 at 21:11

Marcin's code worked fine except a problem in the rasCentrX and rasCentrY functions was causing the ouput coordinates to appear at a different resolution (as Grazia observed). My fix was to change

coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)


to

coordX = rasXmin + ((0.5 + rasWidth) * rasMeanCellWidth)


and

  coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)


to

  coordY = rasYmin + ((0.5 + rasHeight) * rasMeanCellHeight)


I used the code to convert an ESRI Grid to a CSV file. This was achieved by removing the reference to inRaster2, then using a csv.writer to output the coordinates and values:

out = csv.writer(open(outputCSV,"wb"), delimiter=',', quoting=csv.QUOTE_NONNUMERIC)
out.writerow(['X','Y','Value'])
(height, width, dim )=inRasterFullArray.shape
for row in range(0,height):
for col in range(0,width):
out.writerow(inRasterFullArray[row,col])


I also didn't find the transpose was needed in

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))


so converted that to

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray))

-