"Built-in" rasters in python map algebra

5555
10
12-22-2009 02:31 PM
curtvprice
MVP Esteemed Contributor

Any progress in implementing something like Map Algebra's

$$ROWMAP: Row address of the processing cell
$$COLMAP: Column address of the processing cell
$$XMAP: x-map coordinate for the center of the processing cell
$$YMAP: y-map coordinate for the center of the processing cell

0 Kudos
10 Replies
EricRice
Esri Regular Contributor
Hi Curtis,

This functionality will not be provided at  Beta 2 using Map Algebra in the Python Window.  You can however, script it just as you did in versions prior to 9.4.  All pre-9.4 scripts/models will still execute just fine.
0 Kudos
curtvprice
MVP Esteemed Contributor
I had a nice discussion with a product team member (who shall remain anonymous) about this a couple of weeks ago.

He suggested getting conversant with the numpy raster object -- arcpy will have hooks to easily go back and forth, which may be a really handy way to create function-based rasters that don't lend themselves easily to map algebra, or are slower to manipulate within the map algebra environment.

Is this "hook" to numpy rasters going to be available in beta 2? If not what is a workaround to try in the meantime?
0 Kudos
EricRice
Esri Regular Contributor
Hi Curtis,

I'm looking into this for you and hopefully will be able to post later this afternoon.

Thanks,
0 Kudos
RyanDeBruyn
Esri Contributor
RasterToNumpyArray and NumpyArrayToRaster functions are are available in arcpy in 9.4 beta2.  Be aware there are some known issues but in general they are working.

Here is a basic syntax

RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})

NumPyArrayToRaster (in_array, {lower_left_corner}, {x_cell_size}, {y_cell_size}, {value_to_nodata})



link to beta web help:
http://helpdev.esri.com/EN/ArcGISDesktop/9.4/Help/index.html#/NumPyArrayToRaster/000V00000130000000/
http://helpdev.esri.com/EN/ArcGISDesktop/9.4/Help/index.html#/RasterToNumPyArray/000V0000012Z000000/



-Ryan
0 Kudos
curtvprice
MVP Esteemed Contributor
I don't seem to have access to helpdev.esri.com.

Is there another site that is accessible to beta testers outside ESRI?
0 Kudos
KevinHibma
Esri Regular Contributor
Ooops!

You can just take the "dev" out of the link and it should line up the same.

I believe these are the external version of the links Ryan was posting.

http://help.esri.com/EN/ArcGISDesktop/9.4/Help/index.html#/NumPyArrayToRaster/000V00000130000000/
http://help.esri.com/EN/ArcGISDesktop/9.4/Help/index.html#/RasterToNumPyArray/000V0000012Z000000/
0 Kudos
curtvprice
MVP Esteemed Contributor
That was it! Those links work.  Thanks!
0 Kudos
curtvprice
MVP Esteemed Contributor
A mole from the geoprocessing team sent me this neat example.

This method could be used to generate $$XMAP and $$YMAP using the environment raster extent and cell size.

Sure would be nice to include a script tool in the Arc 10 toolbox to both serve as a simple example of how to use NumpyArrayToRaster, and keep the functionality in the standard product. Otherwise the handfull of geeks that use this functionality will have to do it ourselves. Here's the (now trivial, with the example below) tool I will probably write if you don't release one:

CreateLocationRaster <Out_raster> {ROWMAP | COLMAP | XMAP | YMAP} {Template_raster}

## NumPyArrayToRaster (in_array, {lower_left_corner},
##       {x_cell_size}, {y_cell_size}, {value_to_nodata})

import os, sys, arcpy 
import numpy as NUM
 
# Set path
sPath = sys.path[0]
outPath = os.path.join(os.path.dirname(sPath),"output")
 
# Set environment
arcpy.env.workspace = outPath
arcpy.env.overwriteOutput = 1
 
#------------------------------------------
# Create and populate an array
myArray = NUM.zeros((10,16),dtype='int16')
 
# Fill row with row number
for i in range(0,10):
  myArray = i + 1
print myArray
 
# Create a raster
rowRas = arcpy.NumPyArrayToRaster(myArray)
rowRas.save("rowRas")
print "...saved rowras"
 
#------------------------------------------
# Fill col with col number
for j in range(0,16):
  myArray[:,j] = j+1
print myArray

# Create a raster
colRas = arcpy.NumPyArrayToRaster(myArray)
colRas.save("colRas")
print "...saved colras"


And here's another (certifiably cool) way to do it, sent by a USGS colleague... using the NumPy "indices" method:

    import numpy as np
    import arcpy
    # Setup some rasters of rows and columns
    # (like old $$ROWMAP and $$COLMAP)
    nprows = np.indices((10,10))[0]
    npcols = np.indices((10,10))[1]

    # Convert the numpy arrays to ESRI rasters (ie Raster objects)
    # called 'rows' and 'cols'
    rows = arcpy.NumPyArrayToRaster(nprows)
    cols = arcpy.NumPyArrayToRaster(npcols)
0 Kudos
MarkCarroll
New Contributor
Has there been any movement on this in terms of including a script or built-in tool in Arc 10?  If not, can you tell me is the script above intended to be run in the raster calculator or as a script in a model or toolbox?
0 Kudos