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.
- Maven >= v3.0
Maven is the build tool used to build and package the final plugin file. - Python >= 2.7 (recommend is >=3.3)
- OTB 5.6 with python wrapping
- Python Numpy
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.
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
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
- http://otbcb.readthedocs.io/en/latest/Applications/app_RadiometricIndices.html
- https://www.orfeo-toolbox.org/SoftwareGuide/SoftwareGuidech2.html#x16-190002