Using OTB Python API inside SNAP Operator

Introduction

This page explain how to access OTB C++ application via Python API inside a SNAP Operator Plugin.

pre-requsities

Read how to write a processor in python for a tutorial. - How to write a processor in Python. This tutorial explain all the snap related stuff in creating an java maven project to write a SNAP processor

Download and install OTB with python wrapping (see https://www.orfeo-toolbox.org/SoftwareGuide/SoftwareGuidech2.html#x16-190002).

Content

Here is a snap-engine-python-operator which use OTB's RadiometricIndices application. https://github.com/senbox-org/snap-examples/tree/master/snap-engine-otb-python-operator

This example demonstrates computing ndvi for a source product opened in SNAP through OTB. Note that RadiometricIndices application in OTB is not exclusively for computing NDVI. see its documentation here[link].

The python-operator-plugin ask for red and nir bands that will be used to compute NDVI.  The parameter options are handled by ndvi_op-in.xml

Project Setup

In order to write a Python plugin for SNAP you need the following tools.

Project structure


The recommended basic project structure looks like as shown in the following code block. We will look at the content and the meaning of each of those files later.

<PROJECT_DIR>
|   pom.xml (the project definition file for maven)
\---src
    \---main
        +---nbm
        |       manifest.mf (settings for the plugin)
        +---python
        |       my_python_op-info.xml (metadata about the operator)
        |       my_python_op.py (the actual operator implementation)
        \---resources
            |   layer.xml (integration into SNAP Desktop)
            \---META-INF
                \---services
                        org.esa.snap.python.gpf.PyOperatorSpi (make the operator known to SNAP Engine)

Operator implementation

An operator implemented in Python needs to have three methods. An initialize, a compute and a dispose method. Please have a look at the following skeleton on an operator.

skeleton
 class MyOp:

    def init(self): pass

def initialize(self, context):
        pass
  
    # You either implement this or the computeTileStack method
    def computeTile(self, context, band, tile):
        pass
  
    # You either implement this or the computeTile method
    def computeTileStack(self, context, target_tiles, target_rectangle):
        pass
  
    def dispose(self, context):
        pass

__init__ method

def __init__(self):
   self.red_band = None
   self.nir_band = None
   self.ndvi_band = None
   self.ndvi_flags_band = None

   self.otb_ndvi_app = otbApplication.Registry.CreateApplication("RadiometricIndices")

We create our otb application in the __init__ method. Variables used to store input and output band are set to None in __init__ method.

Initialize method

initialize
def initialize(self, context):

        source_product = context.getSourceProduct('source')
        print('initialize: source product location is', source_product.getFileLocation().toString())
        red_band_name = context.getParameter('redName')

        if not red_band_name:
            raise RuntimeError('Missing parameter "redName"')

        self.red_band = self._get_band(source_product, red_band_name)

		 nir_band_name = context.getParameter('nirName')
		 if not nir_band_name:
            raise RuntimeError('Missing parameter "nirName"')

        self.nir_band = self._get_band(source_product, nir_band_name)

        print('initialize: red_band =', self.red_band, ', nir_band =', self.nir_band)

        width = source_product.getSceneRasterWidth()
        height = source_product.getSceneRasterHeight()

        ndvi_product = snappy.Product('py_NDVI', 'py_NDVI', width, height)
        snappy.ProductUtils.copyGeoCoding(source_product, ndvi_product)
        self.ndvi_band = ndvi_product.addBand('ndvi', snappy.ProductData.TYPE_FLOAT32)
        self.ndvi_flags_band = ndvi_product.addBand('ndvi_flags', snappy.ProductData.TYPE_UINT8)
        context.setTargetProduct(ndvi_product)


initialize method takes input from user through UI and create source and target product. We also create band instance which are later used inside computeTileStack method. This method is called only once at start of processing.

ComputeTileStack Method

Data input to OTB application from SNAP is done using python numpy arrays.  Inside computeTileStack method, one can have access to underlying pixels which fall only under the target_rectangle. We pick the pixel data using   Tile.getSamplesFloat() and stack into an N-dimensional numpy array. The call for making numpy.ndarray or numpy.array is decided by the input image type (VectorImage, Image, ImageList) of OTB application.

Below is the snippet which create N-dimensional numpy array using two bands for the input to RadiometricIndices application.

def computeTileStack(self, context, target_tiles, target_rectangle):  
   tile_red = context.getSourceTile(self.red_band, target_rectangle)
   tile_nir = context.getSourceTile(self.nir_band, target_rectangle)

    samples_red = tile_red.getSamplesFloat()
    samples_nir = tile_nir.getSamplesFloat()

    ndvi_tile = target_tiles.get(self.ndvi_band)
    w = ndvi_tile.getWidth()
    h = ndvi_tile.getHeight()

    otb_input_array = np.dstack((np.array(samples_red, dtype=np.float32).reshape(w, h),
                                 np.array(samples_nir, dtype=np.float32).reshape(w, h)))


samples_red and sample_nir represent the underlying array of samples for red and nir bands.

w and h and width and height of single target tile respectively.

Set Input of OTB application from Numpy array

Once we have data we can set this numpy directly as the input parameter to OTB application.

self.otb_ndvi_app.SetVectorImageFromNumpyArray("in", otb_input_array)


The above code will set the input parameter to a otbVectorImage whose each Compoenent represent the values of numpy array.

Setting other parameters for OTB application

Setting other parameters of OTB application.

self.otb_ndvi_app.SetParameterStringList("list", ["Vegetation:NDVI"])
self.otb_ndvi_app.SetParameterInt("channels.red", 1)
self.otb_ndvi_app.SetParameterInt("channels.nir", 2)


The above code will tell RadiometricIndices application to compute only Vegetation:NDVI, where the 1st and 2nd band in the input image (numpy array) represent Red and NIR bands respectively.

Execute OTB application

self.otb_ndvi_app.Execute()


Store Output of OTB application into target product

ndvi = self.otb_ndvi_app.GetVectorImageAsNumpyArray("out")
ndvi_tile.setSamples(ndvi)

Note: Tile.setSamples() requires a N-dimensional numpy array even though the function set samples for one single band. This is a limitation in snap python bindings.

Our otb application returns an N-D array from GetVectorImageAsNumpyArray() but number of bands is 1 because we compute only one index NDVI.

suppose we compute two indices NDVI and WVI, then the "ndvi" will have two bands each representing one index. In that case, we have to workaround to overcome the snap python issue before calling Tile.setSamples()

See  below code snippet which create a temporary nd array and fill the first slice from ndvi

# ndvi_1 = np.ndarray(shape=(w, h, 1), dtype=np.float32)
# ndvi_1[..., 0] = ndvi[:,:,1]

Full source code

https://github.com/senbox-org/snap-examples/tree/master/snap-engine-otb-python-operator

Testing

Compile and Install OTB (instructions for linux)

you must install swig, python-devel and python-numpy package

cd $HOME/orfeotoolbox

wget https://www.orfeo-toolbox.org/packages/xdk/OTB-5.4/OTB-5.4.0-xdk-Linux64.run

chmod +x OTB-5.4.0-xdk-Linux64.run 

./OTB-5.4.0-xdk-Linux64.run

export CMAKE_PREFIX_PATH=$HOME/orfeotoolbox/OTB-5.4.0-xdk-Linux64

git clone https://git@git.orfeo-toolbox.org/git/otb.git

mkdir OTB.build && cd OTB.build

cmake -DOTB_WRAP_PYTHON=TRUE -DOTB_USE_MUPARSER=TRUE ../OTB

make && make install (default install is /usr/local)

Clone snap-examples from github fork

You can checkout the sources and import the project.

git clone https://github.com/rkanavath/snap-examples

Configure Paths


export PYTHONPATH=/usr/local/lib/otb/python
export OTB_APPLICATION_PATH=/usr/local/lib/otb/applications


Start Intelij IDE

see How to develop an extension module for configuring IDE and running the processor inside snap desktop.

Debugging and Running your Extension in an IDE  →  IntelliJ IDEA IDE → Debugging and Running with SNAP Desktop

Open snap-engine-otb-python-operator/pom.xml

References