README

This project is not maintained.

This package enables simple tiled (or untiled if desired) raster calculations (AKA "map algebra") from the commandline or from within your python scripts.

• Notes
• Installation
• Requirements
• Commandline raster calculator
• Raster calculations library

Notes ^

Both the raster calculator and the calculations library can handle rasters with different extents,cellsizes and coordinate systems as long as they overlap.

If extents/cellsizes/coordinate systems differ, the output extent/cellsize will be the MINOF/MAXOF of input datasets, while the output coordinate system will be . that of the leftmost Dataset in the expression unless Env.extent/Env.cellsize/Env.srs are specified.

gdal.Dataset and gdal.RasterBand and numpy.ndarray method and attribute calls are passed down to the underlying gdal.Dataset, gdal.RasterBand and ndarray objects.

If numexpr is installed, it can be used to evaluate your expressions as it is much faster. However, numexpr expressions are very limited: tiled processing, on-the-fly reprojection, extent clipping/snapping, method/function calls and subscripting are not supported.

Installation ^

Windows: Download and run the 32bit or 64bit installer as appropriate. All: Download and extract the gdal-calculations-{version}.tar.gz and run python setup.py install with administrative or root privileges.

Requirements ^

• Python 2.6+
• GDAL with python bindings
• numpy

Commandline raster calculator ^

Name: gdal_calculate
Purpose: Perform simple tiled raster calculations (AKA "map algebra")
         from the commandline

Required parameters:
     --calc     : calculation in numpy syntax, rasters specified as using
                  any legal python variable name syntax, band numbers are
                  specified using square brackets (zero based indexing)
     --outfile  : output filepath
     --{*}      : filepaths for raster variables used in --calc
                  e.g. --calc='(someraster[0]+2)*c' --someraster='foo.tif' --c='bar.tif'

Optional parameters:
     --of            : GDAL format for output file (default "GTiff")')
     --co            : Creation option to the output format driver.
                       Multiple options may be listed.
     --cellsize      : one of DEFAULT|MINOF|MAXOF|"xres yres"|xyres
                       (Default=DEFAULT, leftmost dataset in expression)
     --extent        : one of MINOF|INTERSECT|MAXOF|UNION|"xmin ymin xmax ymax"
                       (Default=MINOF)
     --nodata        : handle nodata using masked arrays (Default=False)
                       uses numpy.ma.MaskedArray to handle NoData values
                       MaskedArrays can be much slower...
     --notile        : don't use tiled processing, faster but uses more memory (Default=False)
     --numexpr       : Enable numexpr evaluation (Default=False)
     --overwrite     : overwrite if required (Default=False)
     -q --quiet      : Don't display progress (Default=False)
     --reproject     : reproject if required (Default=False)
                       datasets are projected to the SRS of the first input
                       dataset in an expression
     --resampling    : one of "AVERAGE"|"BILINEAR"|"CUBIC"|"CUBICSPLINE"|
                       "LANCZOS"|"MODE"|"NEAREST"|gdal.GRA_*)
                       (Default="NEAREST")
    --snap           : filepath of a raster to snap extent coordinates to.
    --srs            : the output spatial reference system
                       one of osgeo.osr.SpatialReference (object)|WKT (string)|EPSG code (integer)
                       (Default = None)
    --tempdir        : filepath to temporary working directory (can also use /vsimem for in memory tempdir)
    --tempoptions    : list of GTIFF creation options to use when creating temp rasters
                       (Default = ['BIGTIFF=IF_SAFER'])

Example:
       gdal_calculate --outfile=../testdata/ndvi.tif       \
            --calc="((nir[3]-red[2].astype(numpy.float32))/(nir[3]+red[2].astype(numpy.float32)))" \
            --red=../testdata/landsat_utm50.tif  \
            --nir=../testdata/landsat_geo.tif    \
            --overwrite --reproject --extent=MAXOF
            

Raster calculations library ^

Name: gdal_calculations
Purpose: GDAL Dataset and Band abstraction for simple tiled raster calculations
         (AKA "map algebra")

Author: Luke Pinner
Contributors: Matt Gregory

Classes/Objects:
    Dataset(filepath_or_dataset ,*args)
        - Base Dataset class.
        - Instantiate by passing a path or gdal.Dataset object.
        - Supports gdal.Dataset and numpy.ndarray method and attribute calls.
        - Supports arithmetic operations (i.e ds1 + ds2)
    Band
        - Returned from Dataset[i] (zero based) or Dataset.GetRasterBand(j) (1 based)
          methods, not instantiated directly.
        - Supports gdal.RasterBand and numpy.ndarray method and attribute calls.
        - Supports arithmetic operations (i.e ds1[0] + ds2.GetRasterBand(1))
    ClippedDataset(dataset_or_band, extent)
        - Subclass of Dataset.
        - Uses VRT functionality to modify extent.
    ConvertedDataset(dataset_or_band, datatype)
        - Subclass of Dataset.
        - Uses VRT functionality to modify datatype.
        - Returned by the type conversion functions, not instantiated directly.
    TemporaryDataset(cols,rows,bands,datatype,srs='',gt=[],nodata=[])
        - Subclass of Dataset.
        - A temporary raster that only persists until it goes out of scope.
        - Stored in the Env.tempdir directory, which may be on disk or in
          memory (/vsimem).
        - Can be made permanent with the `save` method.
    WarpedDataset(dataset_or_band, wkt_srs, snap_ds=None, snap_cellsize=None)
        - Subclass of Dataset.
        - Uses VRT functionality to warp Dataset.
    ArrayDataset(array,extent=[],srs='',gt=[],nodata=[], prototype_ds=None)
        - Subclass of TemporaryDataset.
        - Instantiate by passing a numpy ndarray and georeferencing information
          or a prototype Dataset.
    DatasetStack(filepaths, band=0)
        - Stack of bands from multiple datasets
        - Similar to gdalbuildvrt -separate etc... functionality, except the class
          can handle rasters with different extents,cellsizes and coordinate systems
          as long as they overlap.
    Env - Object for setting various environment properties.
        - This is instantiated on import.
        - The following properties are supported:
            cellsize
              - one of 'DEFAULT','MINOF','MAXOF', [xres,yres], xyres
              - Default = "DEFAULT"
            enable_numexpr
              - this can break core numpy methods, such as numpy.sum([Dataset(foo),Dataset(bar)]
              - Default = False
            extent
              - one of "MINOF", "INTERSECT", "MAXOF", "UNION", [xmin,ymin,xmax,ymax]
              - Default = "MINOF"
            nodata
              - handle nodata using masked arrays - True/False
              - Default = False
            ntiles
              - number of tiles to process at a time
              - Default = 1
            overwrite
              - overwrite if required - True/False
              - Default = False
            reproject
              - reproject if required - True/False
              - datasets are projected to the SRS of the first input dataset in an expression
              - Default = False
            resampling
              - one of "AVERAGE"|"BILINEAR"|"CUBIC"|"CUBICSPLINE"|"LANCZOS"|"MODE"|"NEAREST"|gdal.GRA_*)
              - Default = "NEAREST"
            snap
              - a gdal_calculations.Dataset/Band object
              - Default = None
            srs
              - the output spatial reference system
              - one of osgeo.osr.SpatialReference (object)|WKT (string)|EPSG code (integer)
              - Default = None
            tempdir
              - temporary working directory
              - Default = tempfile.tempdir
            tempoptions
              - list of GTIFF creation options to use when creating temp rasters
              - Default = ['BIGTIFF=IF_SAFER']
            tiled
              - use tiled processing - True/False
              - Default = True
    Byte, UInt16, Int16, UInt32, Int32, Float32, Float64
        - Type conversions functions
        - Returns a ConvertedDataset object

Examples:
    from gdal_calculations import *
    
    Env.extent = [xmin, ymin, xmax, ymax] # Other acceptable values:
                                          #  'INTERSECT' or 'MINOF' (default)
                                          #  'UNION' or 'MAXOF'
    Env.resampling = 'CUBIC'              # Other acceptable values:
                                          #  'NEAREST' (default)
                                          #  'AVERAGE'|'BILINEAR'|'CUBIC'
                                          #  'CUBICSPLINE'|'LANCZOS'|'MODE'
                                          #   gdal.GRA_* constant

    Env.reproject=True  #reproject on the fly if required

    Env.nodata=True  #Use a numpy.ma.MaskedArray to handle NoData values
                      #Note MaskedArrays are much slower...

    Env.overwrite=True

    gdal.UseExceptions()

    ds1=Dataset('../testdata/landsat_utm50.tif')#Projected coordinate system
    ds2=Dataset('../testdata/landsat_geo.tif')  #Geographic coordinate system

    #red=ds1[2].astype(np.float32) #You can use numpy type conversion (is slower)
    red=Float32(ds1[2]) #or use one of the provided type conversion functions (quicker as they use VRT's)
    nir=ds2[3]

    ndvi=(nir-red)/(nir+red)

    #Or in one go
    #ndvi=(ds2[3]-Float32(ds1[2])/(ds2[3]+Float32(ds1[2]))
    ndvi=ndvi.save(r'../testdata/ndvi1.tif')

    #If you want to speed things up... use numexpr!
    #but there are a few limitations...
    import numexpr as ne

    #Must not be tiled for numexpr
    Env.tiled=False

    #No subscripting or methods in the expression
    #red=ds1[2].astype(np.float32)
    red=Float32(ds1[2])
    nir=ds2[3] #Some Int*/UInt* datasets cause segfaults, workaround is cast to Float32

    #Must be same coordinate systems and dimensions
    #The check_extent method will reproject and clip if required
    #This is done using virtual rasters (VRT) so is very quick
    nir,red=nir.check_extent(red)

    expr='(nir-red)/(nir+red)'
    ndvi=ne.evaluate(expr)

    #evaluate returns an ndarray not a Dataset
    #So need to write to a Temporary ArrayDataset
    ndvi=ArrayDataset(ndvi,prototype_ds=nir)
    ndvi=ndvi.save(r'../testdata/ndvi2.tif',options=['compress=LZW','TILED=YES'])

    #Get the raw numpy array data
    for block in red.ReadBlocksAsArray():
        print block.x_off,block.y_off,block.data.shape

    rawdata=red.ReadAsArray()
    print rawdata.shape