arcpy-numpy-blocked.py

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