README

This project is not maintained. 

<a id="top"></a>
This package enables simple tiled (or untiled if desired) raster calculations 
(AKA "map algebra") from the commandline or from within your python scripts.  

* [Notes](#not)
* [Installation](#ins)
* [Requirements](#req)
* [Commandline raster calculator](#cmd)
* [Raster calculations library](#lib)

<a id="not"></a>
Notes [^](#top)
-----
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.

<a id="ins"></a>
Installation [^](#top)
------------

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.

<a id="req"></a>
Requirements [^](#top)
------------

* [Python](http://www.python.org) 2.6+
* [GDAL](http://www.gdal.org) with python bindings
* [numpy](http://www.numpy.org)

<a id="cmd"></a>
Commandline raster calculator [^](#top)
-----------------------------

    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
                

<a id="lib"></a>
Raster calculations library [^](#top)
---------------------------

    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