Skip to content

Instantly share code, notes, and snippets.

@scw
Created June 17, 2016 21:50
Show Gist options
  • Save scw/106677c4ba0c60fec9a806619814b04d to your computer and use it in GitHub Desktop.
Save scw/106677c4ba0c60fec9a806619814b04d to your computer and use it in GitHub Desktop.
# Second example from RasterToNumPyArray (arcpy)
# http://resources.arcgis.com/en/help/main/10.2/index.html#//03q300000029000000
"""
A Python NumPy array is designed to deal with large arrays. There are many exist
ing Python functions that have been created to process NumPy arrays, the most no
ted being contained in the SciPy scientific computing package for Python. You ma
y want to convert an ArcGIS raster to a NumPy array to
1. Implement one of the many existing Python functions that can be applied to a
NumPy array (for example, run filters on the data, perform multidimensional
analysis, or utilize optimization routines).
2. Develop a custom function by accessing the individual cells within the NumPy
array (for example, to implement neighborhood notation, change individual
cell values, or run accumulative operators on an entire raster).
If the array def'n (the lower left corner and the number of rows and columns)
exceeds the extent of the in_raster, the array values will be assigned NoData.
If the lower_left_corner does not coincide with the corner of a cell, it will au
tomatically be snapped to the lower left of the nearest cell corner applying the
same rules as the Snap Raster environment setting. This snapping action within
the RasterToNumPy function is not to be confused with the Snap Raster environmen
t setting; the function only uses the same interaction.
For more information, see the following:
- Learn more about how the Snap Raster environment works
- RasterToNumPyArray supports the direct conversion of multiband rasters to an
N-dimensional array (ndarray).
If the input Raster instance is based on a multiband raster, it returns an ndarr
y, where the length of the first dimension represents the number of bands. The n
darray will have the dimensions (bands, rows, columns).
If the input Raster instance is based on a single raster or a specific band from
a multiband raster, it returns a two-dimensional array with the dimensions (row
s, columns).
RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})
...
The other helper functions:
ExtendTable
Join the contents of a NumPy structured array to a table based on a common attribute field.
FeatureClassToNumPyArray
Convert a feature class to NumPy structured array.
NumPyArrayToFeatureClass
Convert a NumPy structured array to a feature class.
NumPyArrayToTable
Convert a NumPy structured array to a table.
TableToNumPyArray
Convert a table to NumPy structured array.
Data types:
Single numpy.float32
Double numpy.float64
SmallInteger numpy.int32
Integer numpy.int32
OID numpy.int32
GUID <U64
String <u1, <u10, and so on
# XXX look at this as a potential bug: numpy.int16 should be mapped to SmallInteger, right? it's 2 bytes (16bit) same as the
# SmallInteger in ArcGIS world
Do we have an Esri internal representation for 64-bit integers? [64 bit floats map to doubles]. Nothing looks to be exposed at the ArcPy level at least, but would take research. Apparently, GDAL also doesn't directly support 64-bit integers, as [seen here](https://trac.osgeo.org/gdal/wiki/rfc31_ogr_64). An ArcGIS FAQ on [the problem as documented with SDE](http://support.esri.com/es/knowledgebase/techarticles/detail/35865).
"""
# Note that, if the input raster is multiband, the data blocks will also be
# multiband, having dimensions (bands, rows, columns). Otherwise, they will
# have dimensions (rows, columns).
import arcpy
import numpy
import os
# Input raster
filein = os.path.join(os.getcwd(),r"input\input.tif")
# Output raster (after processing)
fileout = os.path.join(os.getcwd(),r"output\blockprocessingrdb22.tif")
# Size of processing data block
# where memorysize = datatypeinbytes*nobands*blocksize^2
blocksize = 512
# ----------------------------------------------------------------------------
# Create raster object from file
myRaster = arcpy.Raster(filein)
# Set environmental variables for output
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = filein
arcpy.env.cellSize = filein
# Loop over data blocks
filelist = []
blockno = 0
for x in range(0, myRaster.width, blocksize):
for y in range(0, myRaster.height, blocksize):
# Lower left coordinate of block (in map units)
mx = myRaster.extent.XMin + x * myRaster.meanCellWidth
my = myRaster.extent.YMin + y * myRaster.meanCellHeight
# Upper right coordinate of block (in cells)
lx = min([x + blocksize, myRaster.width])
ly = min([y + blocksize, myRaster.height])
# noting that (x, y) is the lower left coordinate (in cells)
# Extract data block
myData = arcpy.RasterToNumPyArray(myRaster, arcpy.Point(mx, my),
lx-x, ly-y)
# PROCESS DATA BLOCK -----------------------------
# e.g. Calculate mean of each cell of all bands.
myData -= numpy.mean(myData, axis=0, keepdims=True)
# ------------------------------------------------
# Convert data block back to raster
myRasterBlock = arcpy.NumPyArrayToRaster(myData, arcpy.Point(mx, my),
myRaster.meanCellWidth,
myRaster.meanCellHeight)
# Save on disk temporarily as 'filename_#.ext'
filetemp = ('_%i.' % blockno).join(fileout.rsplit('.',1))
myRasterBlock.save(filetemp)
# Maintain a list of saved temporary files
filelist.append(filetemp)
blockno += 1
# Mosaic temporary files
arcpy.Mosaic_management(';'.join(filelist[1:]), filelist[0])
if arcpy.Exists(fileout):
arcpy.Delete_management(fileout)
arcpy.Rename_management(filelist[0], fileout)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment