Chapter 7
Reading and Writing Auxiliary Data

As we have seen in the previous chapter, OTB has a great capability to read and process images. However, images are not the only type of data we will need to manipulate. Images are characterized by a regular sampling grid. For some data, such as Digital Elevation Models (DEM) or Lidar, this is too restrictive and we need other representations.

Vector data are also used to represent cartographic objects, segmentation results, etc: basically, everything which can be seen as points, lines or polygons. OTB provides functionnalities for accessing this kind of data.

7.1 Reading DEM Files

The source code for this example can be found in the file
Examples/IO/DEMToImageGenerator.cxx.

The following example illustrates the use of the otb::DEMToImageGenerator class. The aim of this class is to generate an image from the srtm data (precising the start extraction latitude and longitude point). Each pixel is a geographic point and its intensity is the altitude of the point. If srtm doesn’t have altitude information for a point, the altitude value is set at -32768 (value of the srtm norm).

Let’s look at the minimal code required to use this algorithm. First, the following header defining the otb::DEMToImageGenerator class must be included.

  #include "otbDEMToImageGenerator.h"

The image type is now defined using pixel type and dimension. The output image is defined as an otb::Image .

    const unsigned int Dimension = 2;
    typedef otb::Image<double, Dimension> ImageType;

The DEMToImageGenerator is defined using the image pixel type as a template parameter. After that, the object can be instancied.

    typedef otb::DEMToImageGenerator<ImageType> DEMToImageGeneratorType;
  
    DEMToImageGeneratorType::Pointer object = DEMToImageGeneratorType::New();

Input parameter types are defined to set the value in the otb::DEMToImageGenerator .

    typedef DEMToImageGeneratorType::SizeType       SizeType;
    typedef DEMToImageGeneratorType::SpacingType    SpacingType;
    typedef DEMToImageGeneratorType::PointType      PointType;

The path to the DEM folder is given to the otb::DEMHandler .

    otb::DEMHandler::Instance()->OpenDEMDirectory(folderPath);

The origin (Longitude/Latitude) of the output image in the DEM is given to the filter.

    PointType origin;
    origin[0] = ::atof(argv[3]);
    origin[1] = ::atof(argv[4]);
  
    object->SetOutputOrigin(origin);

The size (in Pixel) of the output image is given to the filter.

    SizeType size;
    size[0] = ::atoi(argv[5]);
    size[1] = ::atoi(argv[6]);
  
    object->SetOutputSize(size);

The spacing (step between to consecutive pixel) is given to the filter. By default, this spacing is set at 0.001.

    SpacingType spacing;
    spacing[0] = ::atof(argv[7]);
    spacing[1] = ::atof(argv[8]);
  
    object->SetOutputSpacing(spacing);

The output image name is given to the writer and the filter output is linked to the writer input.

    writer->SetFileName(outputName);
  
    writer->SetInput(object->GetOutput());

The invocation of the Update() method on the writer triggers the execution of the pipeline. It is recommended to place update calls in a try/catch block in case errors occur and exceptions are thrown.

    try
      {
      writer->Update();
      }
  
    catch (itk::ExceptionObject& err)
      {
      std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
      std::cout << err << std::endl;
      return EXIT_FAILURE;
      }

Let’s now run this example using as input the SRTM data contained in DEM_srtm folder. Figure 7.1 shows the obtained DEM. Invalid data values – hidden areas due to SAR shadowing – are set to zero.


PIC

Figure 7.1: DEMToImageGenerator image.


7.2 Elevation management with OTB

The source code for this example can be found in the file
Examples/IO/DEMHandlerExample.cxx.

OTB relies on OSSIM for elevation handling. Since release 3.16, there is a single configuration class otb::DEMHandler to manage elevation (in image projections or localization functions for example). This configuration is managed by the a proper instantiation and parameters setting of this class. These instantiations must be done before any call to geometric filters or functionalities. Ossim internal accesses to elevation are also configured by this class and this will ensure consistency throughout the library.

This class is a singleton, the New() method is deprecated and will be removed in future release. We need to use the Instance() method instead.

    otb::DEMHandler::Pointer demHandler = otb::DEMHandler::Instance();

It allows configuring a directory containing DEM tiles (DTED or SRTM supported) using the OpenDEMDirectory() method. The OpenGeoidFile() method allows inputting a geoid file as well. Last, a default height above ellipsoid can be set using the SetDefaultHeightAboveEllipsoid() method.

    demHandler->SetDefaultHeightAboveEllipsoid(defaultHeight);
  
    if(!demHandler->IsValidDEMDirectory(demdir.c_str()))
      {
      std::cerr<<"IsValidDEMDirectory("<<demdir<<") = false"<<std::endl;
      fail = true;
      }
  
    demHandler->OpenDEMDirectory(demdir);
    demHandler->OpenGeoidFile(geoid);

We can now retrieve height above ellipsoid or height above Mean Sea Level (MSL) using the methods GetHeightAboveEllipsoid() and GetHeightAboveMSL(). Outputs of these methods depend on the configuration of the class otb::DEMHandler and the different cases are:

For GetHeightAboveEllipsoid():

For GetHeightAboveMSL():

    otb::DEMHandler::PointType point;
    point[0] = longitude;
    point[1] = latitude;
  
    double height = -32768;
  
    height = demHandler->GetHeightAboveMSL(point);
    std::cout<<"height above MSL ("<<longitude<<","
             <<latitude<<") = "<<height<<" meters"<<std::endl;
  
    height = demHandler->GetHeightAboveEllipsoid(point);
    std::cout<<"height above ellipsoid ("<<longitude
             <<", "<<latitude<<") = "<<height<<" meters"<<std::endl;

Note that OSSIM internal calls for sensor modelling use the height above ellipsoid, and follow the same logic as the GetHeightAboveEllipsoid() method.

More examples about representing DEM are presented in section 23.1.4.

7.3 Reading and Writing Shapefiles and KML

The source code for this example can be found in the file
Examples/IO/VectorDataIOExample.cxx.

Unfortunately, many vector data formats do not share the models for the data they represent. However, in some cases, when simple data is stored, it can be decomposed in simple objects as for instance polylines, polygons and points. This is the case for the Shapefile and the KML (Keyhole Markup Language) formats, for instance.

Even though specific reader/writer for Shapefile and the Google KML are available in OTB, we designed a generic approach for the IO of this kind of data.

The reader/writer for VectorData in OTB is able to access a variety of vector file formats (all OGR supported formats)

In section 11.4, you will find more information on how projections work for the vector data and how you can export the results obtained with OTB to the real world.

This example illustrates the use of OTB’s vector data IO framework.

We will start by including the header files for the classes describing the vector data and the corresponding reader and writer.

  #include "otbVectorData.h"
  #include "otbVectorDataFileReader.h"
  #include "otbVectorDataFileWriter.h"

We will also need to include the header files for the classes which model the individual objects that we get from the vector data structure.

  #include "itkPreOrderTreeIterator.h"
  #include "otbObjectList.h"
  #include "otbPolygon.h"

We define the types for the vector data structure and the corresponding file reader.

    typedef otb::VectorData<PixelType, 2> VectorDataType;
  
    typedef otb::VectorDataFileReader<VectorDataType>
    VectorDataFileReaderType;

We can now instantiate the reader and read the data.

    VectorDataFileReaderType::Pointer reader = VectorDataFileReaderType::New();
    reader->SetFileName(argv[1]);
    reader->Update();

The vector data obtained from the reader will provide a tree of nodes containing the actual objects of the scene. This tree will be accessed using an itk::PreOrderTreeIterator .

    typedef VectorDataType::DataTreeType            DataTreeType;
    typedef itk::PreOrderTreeIterator<DataTreeType> TreeIteratorType;

In this example we will only read polygon objects from the input file before writing them to the output file. We define the type for the polygon object as well as an iterator to the vertices. The polygons obtained will be stored in an otb::ObjectList .

    typedef otb::Polygon<double>                     PolygonType;
    typedef otb::ObjectList<PolygonType>             PolygonListType;
  
    PolygonListType::Pointer polygonList = PolygonListType::New();

We get the data tree and instantiate an iterator to walk through it.

    TreeIteratorType it(reader->GetOutput()->GetDataTree());
  
    it.GoToBegin();

We check that the current object is a polygon using the IsPolygonFeature() method and get its exterior ring in order to store it into the list.

    while (!it.IsAtEnd())
      {
      if (it.Get()->IsPolygonFeature())
        {
        polygonList->PushBack(it.Get()->GetPolygonExteriorRing());
        }
      ++it;
      }

Before writing the polygons to the output file, we have to build the vector data structure. This structure will be built up of nodes. We define the types needed for that.

    VectorDataType::Pointer outVectorData = VectorDataType::New();
  
    typedef VectorDataType::DataNodeType DataNodeType;

We fill the data structure with the nodes. The root node is a document which is composed of folders. A list of polygons can be seen as a multi polygon object.

    DataNodeType::Pointer document = DataNodeType::New();
    document->SetNodeType(otb::DOCUMENT);
    document->SetNodeId("polygon");
    DataNodeType::Pointer folder = DataNodeType::New();
    folder->SetNodeType(otb::FOLDER);
    DataNodeType::Pointer multiPolygon = DataNodeType::New();
    multiPolygon->SetNodeType(otb::FEATURE_MULTIPOLYGON);

We assign these objects to the data tree stored by the vector data object.

    DataTreeType::Pointer tree = outVectorData->GetDataTree();
    DataNodeType::Pointer root = tree->GetRoot()->Get();
  
    tree->Add(document, root);
    tree->Add(folder, document);
    tree->Add(multiPolygon, folder);

We can now iterate through the polygon list and fill the vector data structure.

    for (PolygonListType::Iterator pit = polygonList->Begin();
         pit != polygonList->End(); ++pit)
      {
      DataNodeType::Pointer newPolygon = DataNodeType::New();
      newPolygon->SetPolygonExteriorRing(pit.Get());
      tree->Add(newPolygon, multiPolygon);
      }

And finally we write the vector data to a file using a generic otb::VectorDataFileWriter .

    typedef otb::VectorDataFileWriter<VectorDataType> WriterType;
  
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput(outVectorData);
    writer->SetFileName(argv[2]);
    writer->Update();

This example can convert an ESRI Shapefile to a MapInfo File but you can also access with the same OTB source code to a PostgreSQL datasource, using a connection string as : PG:”dbname=’databasename’ host=’addr’ port=’5432’ user=’x’ password=’y’” Starting with GDAL 1.6.0, the set of tables to be scanned can be overridden by specifying tables=schema.table.

7.4 Handling large vector data through OGR

The source code for this example can be found in the file
Examples/IO/OGRWrappersExample.cxx.

Starting with the version 3.14.0 of the OTB library, a wrapper around OGR API is provided. The purposes of the wrapper are:

As OGR already provides a rich set of geometric related data, as well as the algorithms to manipulate and serialize them, we’ve decided to wrap it into a new exception-safe interface.

This example illustrates the use of OTB’s OGR wrapper framework. This program takes a source of polygons (a shape file for instance) as input and produces a datasource of multi-polygons as output.

We will start by including the header files for the OGR wrapper classes, plus other header files that are out of scope here.

  #include "otbOGRDataSourceWrapper.h"

The following declarations will permit to merge the otb::ogr::Field s from each otb::ogr::Feature into list-fields. We’ll get back to this point later.

  #include <string>
  #include <vector>
  #include <boost/variant.hpp>
  #include "otbJoinContainer.h"
  
  typedef std::vector<int>                                    IntList_t;
  typedef std::vector<std::string>                            StringList_t;
  typedef std::vector<double>                                 RealList_t;
  // TODO: handle non recognized fields
  typedef boost::variant<IntList_t, StringList_t, RealList_t> AnyListField_t;
  typedef std::vector<AnyListField_t>                         AnyListFieldList_t;
  
  AnyListFieldList_t prepareNewFields(
    OGRFeatureDefn /⋆const⋆/& defn, otb::ogr::Layer & destLayer);
  void printField(
    otb::ogr::Field const& field, AnyListField_t const& newListField);
  void assignField(
    otb::ogr::Field field, AnyListField_t const& newListFieldValue);
  void  pushFieldsToFieldLists(
    otb::ogr::Feature const& inputFeature, AnyListFieldList_t & field);

We caw now instantiate first the input otb::ogr::DataSource .

      otb::ogr::DataSource::Pointer source = otb::ogr::DataSource::New(
        argv[1], otb::ogr::DataSource::Modes::Read);

And then, we can instantiate the output otb::ogr::DataSource and its unique otb::ogr::Layer multi-polygons.

      otb::ogr::DataSource::Pointer destination = otb::ogr::DataSource::New(
        argv[2], otb::ogr::DataSource::Modes::Update_LayerCreateOnly);
      otb::ogr::Layer destLayer = destination->CreateLayer(
        argv[2], ITK_NULLPTR, wkbMultiPolygon);

The data obtained from the reader mimics the interface of OGRDataSource. To access the geometric objects stored, we need first to iterate on the otb::ogr::Layer s from the otb::ogr::DataSource , then on the otb::ogr::Feature from each layer.

      // for (auto const& inputLayer : ⋆source)
      for (otb::ogr::DataSource::const_iterator lb=source->begin(), le=source->end()
  ; lb != le
  ; ++lb)
        {
        otb::ogr::Layer const& inputLayer = lb;

In this example we will only read polygon objects from the input file before writing them to the output file. As all features from a layer share the same geometric type, we can filter on the layer geometric type.

        if (inputLayer.GetGeomType() != wkbPolygon)
          {
          std::cout << "Warning: Ignoring layer: ";
          inputLayer.PrintSelf(std::cout, 2);
          continue; // skip to next layer
          }

In order to prepare the fields for the new layer, we first need to extract the fields definition from the input layer in order to deduce the new fields of the result layer.

        OGRFeatureDefn & sourceFeatureDefn = inputLayer.GetLayerDefn();
        AnyListFieldList_t fields = prepareNewFields(sourceFeatureDefn, destLayer);

The result layer will contain only one feature, per input layer, that stores a multi-polygon shape. All geometric shapes are plain OGRGeometry objects.

        OGRMultiPolygon destGeometry; // todo: use UniqueGeometryPtr

The transformation algorithm is as simple as aggregating all the polygons from the features from the input layer into the destination multi-polygon geometric object.

Note that otb::ogr::Feature ::GetGeometry() provides a direct access to a non-mutable OGRGeometry pointer and that OGRGeometryCollection::addGeometry() copies the received pointer. As a consequence, the following code is optimal regarding the geometric objects manipulated.

This is also at this point that we fetch the field values from the input features to accumulate them into the fields list.

        // for (auto const& inputFeature : inputLayer)
        for (otb::ogr::Layer::const_iterator fb=inputLayer.begin(), fe=inputLayer.end()
  ; fb != fe
  ; ++fb)
          {
          otb::ogr::Feature const& inputFeature = fb;
          destGeometry.addGeometry(inputFeature.GetGeometry());
          pushFieldsToFieldLists(inputFeature,fields);
          } // for each feature

Then the new geometric object can be added to a new feature, that will be eventually added to the destination layer.

        otb::ogr::Feature newFeature(destLayer.GetLayerDefn());
        newFeature.SetGeometry(&destGeometry); // SetGeom -> copies

We set here the fields of the new feature with the ones accumulated over the features from the input layer.

        for (size_t i=0, N=sourceFeatureDefn.GetFieldCount(); i!=N; ++i)
          {
          printField(newFeature[i], fields[i]);
          assignField(newFeature[i], fields[i]);
          }

Finally we add (with otb::ogr::Layer::CreateFeature() the new feature to the destination layer, and we can process the next layer from the input datasource.

        destLayer.CreateFeature(newFeature); // adds feature to the layer
        } // for each layer

In order to simplify the manipulation of otb::ogr::Field s and to avoid copy-paste for each possible case of field-type, this example relies on boost.Variant.

As such, we have defined AnyListField_t as a variant type on all possible types of field. Then, the manipulation of the variant field values is done through the templatized functions otb::ogr::Field::SetValue<>() and otb::ogr::Field::GetValue<>(), from the various variant-visitors.

Before using the visitors, we need to operate a switch on the exact type of each field from the input layers. An empty field-values container is first added to the set of fields containers. Finally, the destination layer is completed with a new field of the right deduced type.

  AnyListFieldList_t prepareNewFields(
    OGRFeatureDefn /⋆const⋆/& defn, otb::ogr::Layer & destLayer)
  {
    AnyListFieldList_t fields;
    for (int i=0, N=defn.GetFieldCount(); i!=N; ++i)
      {
      const char name = defn.GetFieldDefn(i)->GetNameRef();
      OGRFieldType type  = static_cast<OGRFieldType>(-1);
      switch (defn.GetFieldDefn(i)->GetType())
        {
      case OFTInteger:
        fields.push_back (IntList_t());
        type = OFTIntegerList;
        break;
      case OFTString:
        fields.push_back (StringList_t());
        type = OFTStringList;
        break;
      case OFTReal:
        fields.push_back (RealList_t());
        type = OFTRealList;
        break;
      default:
        std::cerr << "Unsupported field type: " <<
          OGRFieldDefn::GetFieldTypeName(defn.GetFieldDefn(i)->GetType())
          << " for " << name << "\n";
        break;
        }
      OGRFieldDefn newFieldDefn(name, type); // name is duplicated here => no dangling pointer
      destLayer.CreateField(newFieldDefn, false);
      }
    return fields;
  }

The first visitor, PushVisitor(), takes the value from one field and pushes it into a container of list-variant. The type of the field to fetch is deduced from the type of the values stored in the container. It is called by pushFieldsToFieldLists(), that for each field of the input feature applies the visitor on the container.

  struct PushVisitor : boost::static_visitor<>
  {
    PushVisitor(otb::ogr::Field const& f) : m_f(f) {}
  
    template <typename T> void operator()(T & container) const
      {
      typedef typename T::value_type value_type;
      value_type const value = m_f.GetValue<value_type>();
      container.push_back(value);
      }
  private:
    otb::ogr::Field const& m_f;
  };
  
  void  pushFieldsToFieldLists(
    otb::ogr::Feature const& inputFeature, AnyListFieldList_t & fields)
  {
    // For each field
    for (size_t i=0, N=inputFeature.GetSize(); i!=N; ++i)
      {
      otb::ogr::Field field = inputFeature[i];
      boost::apply_visitor(PushVisitor(field), fields[i]);
      }
  }

A second simple visitor, PrintVisitor, is defined to trace the values of each field (which contains a list of typed data (integers, strings or reals).

  struct PrintVisitor : boost::static_visitor<>
  {
    template <typename T> void operator()(T const& container) const
      {
      otb::Join(std::cout, container, ", ");
      }
  };
  
  void printField(otb::ogr::Field const& field, AnyListField_t const& newListField)
  {
    std::cout << field.GetName() << " -> ";
    boost::apply_visitor(PrintVisitor(), newListField);
  }

The third visitor, SetFieldVisitor, sets the field of the destination features, which have been accumulated in the list of typed values.

  struct SetFieldVisitor : boost::static_visitor<>
  {
    SetFieldVisitor(otb::ogr::Field f) : m_f(f) {}
    // operator() from visitors are expected to be const
    template <typename T> void operator()(T const& container) const
      {
      m_f.SetValue(container);
      }
  private:
    otb::ogr::Field mutable m_f; // this is a proxy -> a reference in a sort
  };
  
  void assignField(otb::ogr::Field field, AnyListField_t const& newListFieldValue)
  {
    boost::apply_visitor(SetFieldVisitor(field), newListFieldValue);
  }

Note that this example does not handle the case when the input layers don’t share a same fields-definition.