OTB Mad Lab

It seems that C++ programming can be cumbersome for some users. Well, OTB allows you to work also in other programming languages as for instance Python. The Pylab Project which gathers numeric and scientific Python packages such as NumPy, SciPy and Matplotlib is considered as having the same features as some commercial scientific environments. Some people consider that Pylab is even more powerful than these commercial tools, since you get for free the power of a full featured and general purpose programming language as Python.

Recently I needed to quickly plot the temporal profiles of several pixels coming from temporal image series in order to easily assess the impact of a correction algorithm. I have to admit that my first idea was to generate simple ASCII files with the pixel values and then feed these files to Gnuplot. However, I found a more elegant (at least from my point of view) way to do this.

I use the OTB Python bindings in order to read my images (2 GeoTiff files with 49 bands each). Using the “GetPixel” method, I access the pixels I want to plot. Finally I use pylab to plot 2 temporal profiles (before and after correction). Here goes the code.

import sys

if( len(sys.argv) != 5):
    print("Usage: "+sys.argv[0]+" series1 series2 x0 y0")

import itk
import otbImages
import otbIO
import pylab

# Types

ImageType = otbImages.VectorImage[itk.D, 2]
ReaderType = otbIO.ImageFileReader[ImageType]

# Pixel coordinates

index = itk.Index[2]()
index[0] = int(sys.argv[3])
index[1] = int(sys.argv[4])

# Import first series
readerBefore = ReaderType.New()

pixBefore = readerBefore.GetOutput().GetPixel( index )

# Import second series
readerAfter = ReaderType.New()

pixAfter = readerAfter.GetOutput().GetPixel( index )

# Chek that both series have the same length

if( pixBefore.GetSize() != pixAfter.GetSize() ):
    print("Series have different lengths")

# Plot

vectorBefore = 
 [pixBefore.GetElement(i) for i in range(pixBefore.GetSize())]
vectorAfter = 
 [pixAfter.GetElement(i) for i in range(pixAfter.GetSize())]

pylab.plot(vectorBefore, label="Raw")
pylab.plot(vectorAfter, label="Clean")
pylab.legend(loc='upper left', shadow=True)

And here you can see a result example.

Temporal profiles

Of course, you can imagine to make this script a little bit more interactive in order to choose other pixel coordinates, etc., but I hope you grasp the power of this approach. By using the OTB Python bindings and other freely available open source Python goodies (NumPy, SciPy, Matplotlib, rpy) you get the user friendly approach of Octave, Scilab and their commercial counterparts plus the power of OTB (which in turn makes available ITK, OSSIM, etc.).

Isn’t it great?

4 thoughts on “OTB Mad Lab

  1. Welcome to my favourite language admittedly I do the same with gdal python bindings … Try orange classifiers and learning stuff and hook up otb to it.

    Yay for big snakes

  2. Great example

    I would also like to peform OTB processing tasks using Python as I am no software engineer nor have any experience in C++. I do however have experience using python for geo-processing tasks.

    I am forced to use windows and this is the OS for most of our geoprocessing software

    I have downloaded the OTB wrapping package windows installer and am using python 2.7. When I run the OTB python “hello world” script there is no recognition of the otb modules. “ImportError: No module named itk”

    Not sure how to get python to recognise the OTB modules.

    Any help would be greatly appeciated


  3. Is it possible to load numpy array (image) and use it with other filters?

    For example, I have image of instance python.pillow.Image. I can get numpy array version easily. Now if i want to use this python object as an image source just like a gdalIO/JpegIO. I am asking because itk gives numpy array and not sure if possible with OTB?

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.