Running map algebra in python

Running Map Algebra in Python

MrGeo includes python bindings that make it simple to execute map algebra functionality within a python script. The python bindings provide functions corresponding to all of the functionality described in the map algebra reference with the same parameters.

Install pymrgeo

We are working on a pip installer. Coming soon...

Configuration

Running against MrGeo installed on the same machine

Make sure that pyspark is accessible by python. The easiest way to do that is to set the PYTHONPATH environment variable to include "$SPARK_HOME/python:$SPARK_HOME/python/lib" so that python can find pyspark.

Also make sure that MrGeo is installed and configured properly.

Running against a remote MrGeo instance

Coming soon...

MrGeo Python Syntax

The following is a boilerplate python script for running MrGeo map algebra:

#!/usr/bin/python

from pymrgeo import MrGeo

# Initialize MrGeo context
print("Initializing")
mrgeo = MrGeo()

mrgeo.start()

try:
    # Map Algebra code goes here...

    elev = mrgeo.load_image("aster-30m")
    result = elev.slope(elev)
    result.save("aster-slope-30m")

finally:
    mrgeo.stop()

The map algebra syntax in python is object-oriented as shown in the above sample. An image is initially loaded by calling the MrGeo.load_image function which returns a MrGeo raster object. That object can then be used to invoke any of the raster map algebra functions (like "slope" in the above example). The raster functions also return raster objects which can be used to call other raster functions.

Also of note in using python is that the arithmetic operations (+, -, /, *, and ^ for power) can all be invoked as infix operations on a raster variable - for example:

myImage = mrgeo.load_image("my-image")
myNewImage = myImage + 1

Also be sure to run the "save" method for any results you wish to save.

The MrGeo Python Class

As shown in the above sample, the MrGeo class must be instantiated to provide the context within which MrGeo executes. That class also provides the following functions:

• create_points(lon, lat, lon, lat, ...) - where the arguments consist of pairs of longitude and latitude world coordinate pairs. The return value is a vector map op that can be passed to other map algebra functions that take vector map op arguments.
• ingest_image(name, zoom, categorical) - where name is the image source path, zoom is the optional zoom level at which to ingest the image, and categorical is an optional boolean which should be True if the source image contains categorical data. The return value is a raster map op that can be passed to other map algebra functions that accept image arguments.
• load_image(imageName) - where imageName is the name of an ingested MrGeo image. The return value is a raster map op that can be passed to other map algebra functions that accept image arguments.
• load_vector(vectorName) - where vectorName is the name of a vector source (e.g. a shapefile, tsv file, csv file, or GeoWave data source). The return value is a vector map op that can be passed to other map algebra functions that take vector map op arguments.

Examples

In the mrgeo-python sub-directory within the repo, there are several scripts that show how to make use of pymrgeo. tester.py shows a simple example of how to initialize MrGeo, load an image, run a slope operation, and cleanup MrGeo. For more complex examples, take a look at landsat.py and landsat-ndvi.py.

landsat.py shows how to dynamically ingest the red, green and blue bands from a landat8 scene, convert each to reflectance units, and build and export an RGB output image.

landsat-ndvi.py shows how to loop through multiple landsat8 scenes, and convert bands 4 (red) and 5 (NIR) to reflectance units and store the results to temp files using GDAL and numpy, then dynamically ingest all of those files together as a single MrGeo image, and finally compute the normalized difference vegetation index (NDVI) across that entire image using MrGeo map algebra.

Landsat8 imagery is available in Amazon S3 as a public data set. More information here.