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

Question

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

---

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 using 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 shapefile 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

Answer 1

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

import arcpy, numpy

# Set input rasters
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 + rasMeanCellWidth*(0.5 + rasWidth)
    return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord

def rasCentrY(rasHeight, rasWidth):
    coordY = rasYmin + rasMeanCellHeight*(0.5 + 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.

To save out the coordinates to a csv file, run:

x = inRasterCoordX.T.flatten()
y = inRasterCoordY.T.flatten()
DF = pd.DataFrame({'x': x, 'y': y})
# save the dataframe as a csv file
DF.to_csv(r"your_output_path.csv")

Answer 2

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

Answer 3

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.

Answer 4

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
    arcpy.AddXY_management(shp_path)

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

Answer 5

A simple solution using open source python packages:

import fiona
import rasterio
from pprint import pprint


def raster_point_coords(raster, points):

    # initialize dict to hold data
    pt_data = {}

    with fiona.open(points, 'r') as src:
        for feature in src:
            # create dict entry for each feature
            pt_data[feature['id']] = feature

    with rasterio.open(raster, 'r') as src:
        # read raster into numpy array
        arr = src.read()
        # rasterio always reads into 3d array, this is 2d, so reshape
        arr = arr.reshape(arr.shape[1], arr.shape[2])
        # get affine, i.e. data needed to work between 'image' and 'raster' coords
        a = src.affine

    for key, val in pt_data.items():
        # get coordinates
        x, y = val['geometry']['coordinates'][0], val['geometry']['coordinates'][1]
        # use affine to convert to row, column
        col, row = ~a * (x, y)
        # remember numpy array is indexed array[row, column] ie. y, x
        val['raster_value'] = arr[int(row), int(col)]

    pprint(pt_data) 

if __name__ == '__main__':
    # my Landsat raster
    ras = '/data01/images/sandbox/LT05_040028_B1.tif'
    # my shapefile with two points which overlap raster area
    pts = '/data01/images/sandbox/points.shp'
    # call function
    raster_point_coords(ras, pts)

Fiona is handy as you can open a shapefile, iterate through the features, and (as I have) append them to a dict object. Indeed the Fiona feature itself is like a dict as well, so it is easy to access properties. If my points had any attributes, they would appear in this dict along with the coordinates, id, etc.

Rasterio is handy because it's easy to read in the raster as a numpy array, a light and fast data type. We also have access to a dict of raster properties including the affine, which is all the data we need to convert the raster x, y coordinates into array row, col coordinates. See @perrygeo's excellent explanation here.

We end up with a pt_data of type dict which has data for each point and the extracted raster_value. We could easily rewrite the shapefile with the extracted data as well if we wanted.

Answer 6

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

Answer 7

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.

Answer 8

Ugly but highly effective:

1. Create a new point feature with 4 points outside the corners of raster in question. Make sure in same coordinate system as raster in question.
2. Add 'xcor' and 'ycor' double fields
3. Calculate geometry to get coordinates for these fields
4. Spatial Analyst->Interpolation->Trend->Linear regression
5. Environment settings: snap raster and cell size to same as raster in question
6. Execute separately for 'xcor' and 'ycor'
7. Out comes raters with coordinates as cell values, use as input for scripts.