src/terralib/attributefill/Utils.cpp

Go to the documentation of this file.

/*  Copyright (C) 2008 National Institute For Space Research (INPE) - Brazil.

    This file is part of the TerraLib - a Framework for building GIS enabled applications.

    TerraLib is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License,
    or (at your option) any later version.

    TerraLib is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with TerraLib. See COPYING. If not, write to
    TerraLib Team at <terralib-team@terralib.org>.
 */

/*!
  \file terralib\attributefill\core\Utils.cpp
   
  \brief Utility functions for Attribute Fill.
*/

// TerraLib
#include "../core/translator/Translator.h"
#include "../geometry/Envelope.h"
#include "../geometry/MultiPolygon.h"
#include "../geometry/Polygon.h"
#include "../geometry/Utils.h"
#include "../raster/Band.h"
#include "../raster/BandProperty.h"

#include "Utils.h"

// STL
#include <tuple>

std::string te::attributefill::GetOperationFullName(const int& e)
{
  switch(e)
  {
    case te::attributefill::VALUE:
      return "Value";
    case te::attributefill::MIN_VALUE:
      return "Minimum value";
    case te::attributefill::MAX_VALUE:
      return "Maximum value";
    case te::attributefill::MEAN:
      return "Mean";
    case te::attributefill::SUM:
      return "Sum of values";
    case te::attributefill::COUNT:
      return "Total number of values";
    case te::attributefill::VALID_COUNT:
      return "Total not null values";
    case te::attributefill::COUNT_DISTINCT:
      return "Total number of distinct values";
    case te::attributefill::STANDARD_DEVIATION:
      return "Standard deviation";
    case te::attributefill::VARIANCE:
      return "Variance";
    case te::attributefill::SKEWNESS:
      return "Skewness";
    case te::attributefill::KURTOSIS:
      return "Kurtosis";
    case te::attributefill::AMPLITUDE:
      return "Amplitude";
    case te::attributefill::MEDIAN:
      return "Median";
    case te::attributefill::VAR_COEFF:
      return "Coefficient variation";
    case te::attributefill::MODE:
      return "Mode";
    case te::attributefill::HIGHEST_OCCURRENCE:
      return "Class with highest occurrence";
    case te::attributefill::HIGHEST_INTERSECTION:
      return "Class with highest intersection area";
    case te::attributefill::PERCENT_CLASS:
      return "Percentage per Class";
    case te::attributefill::MIN_DISTANCE:
      return "Minimum Distance";
    case te::attributefill::MIN_DISTANCE_CENTROID:
      return "Minimum Distance From Centroid";
    case te::attributefill::PRESENCE:
      return "Presence";
    case te::attributefill::WEIGHTED:
      return "Weighted by Area";
    case te::attributefill::WEIGHTED_SUM:
      return "Weighted Sum by Area";
    case te::attributefill::PERCENT_EACH_CLASS:
      return "Percentage of each Class by Area";
    case te::attributefill::PERCENT_TOTAL_AREA:
      return "Percentage of Total Area";
    default:
      return "";
  }
}

std::pair<uint32_t, uint32_t> te::attributefill::GeoToGrid(const te::gm::Coord2D& coord, const te::rst::Grid& grid)
{
  double colD=-1, rowD=-1;
  grid.geoToGrid(coord.getX(), coord.getY(), colD, rowD);
  int col = static_cast<int>(std::round(colD));
  int row = static_cast<int>(std::round(rowD));

  if(col < 0 || col >= (int)grid.getNumberOfColumns())
    col = std::numeric_limits<uint32_t>::max();
  if(row < 0 || row >= (int)grid.getNumberOfRows())
    row = std::numeric_limits<uint32_t>::max();

  return std::make_pair(col, row);
}

void te::attributefill::GetMinMaxLineAndColumn(const te::rst::Raster &raster,
                                               const te::gm::Polygon &polygon,
                                               uint32_t &minimumRow, uint32_t &minimumColumn,
                                               uint32_t &maximumRow, uint32_t &maximumColumn)
{
  std::unique_ptr<te::gm::Polygon> polygonTemp(
      dynamic_cast<te::gm::Polygon*>(polygon.clone()));

  if (polygonTemp->getSRID() != raster.getSRID())
    polygonTemp->transform(raster.getSRID());

  auto polygonEnvelope = polygonTemp->getMBR();

  auto lowerLeftCoord2D = polygonEnvelope->getLowerLeft();
  auto upperRightCoord2D = polygonEnvelope->getUpperRight();

  auto grid = raster.getGrid();

  std::tie(minimumColumn, maximumRow) = GeoToGrid(lowerLeftCoord2D, *grid);
  if(minimumColumn == std::numeric_limits<uint32_t>::max())
    minimumColumn = 0;

  if(maximumRow == std::numeric_limits<uint32_t>::max())
    maximumRow = grid->getNumberOfRows()-1;

  std::tie(maximumColumn, minimumRow) = GeoToGrid(upperRightCoord2D, *grid);
  if(minimumRow == std::numeric_limits<uint32_t>::max())
    minimumRow = 0;

  if(maximumColumn == std::numeric_limits<uint32_t>::max())
    maximumColumn = grid->getNumberOfColumns()-1;
}

void te::attributefill::GetValuesFromBand(const te::rst::Raster &raster,
                                          const unsigned int &band,
                                          const te::gm::Polygon &polygon,
                                          const uint32_t &minimumRow,
                                          const uint32_t &minimumColumn,
                                          const uint32_t &maximumRow,
                                          const uint32_t &maximumColumn,
                                          std::map<double, int> &values)
{
  assert(band < raster.getNumberOfBands());

  auto noDataValue = raster.getBand(band)->getProperty()->m_noDataValue;

  std::unique_ptr<te::gm::Geometry> polygonTemp(
      dynamic_cast<te::gm::Geometry*>(polygon.clone()));

  if (polygonTemp->getSRID() != raster.getSRID())
    polygonTemp->transform(raster.getSRID());

  auto grid = raster.getGrid();

  auto xResolution = grid->getResolutionX();
  auto yResolution = grid->getResolutionY();

  for(uint32_t row = minimumRow; row <= maximumRow; ++row)
  {
    for(uint32_t col = minimumColumn; col <= maximumColumn; ++col)
    {
      auto pixelCenter = grid->gridToGeo(col, row);

      std::unique_ptr<te::gm::Envelope> pixelEnvelope( new te::gm::Envelope(pixelCenter.getX()-(xResolution/2.),
                                                                            pixelCenter.getY()-(yResolution/2.),
                                                                            pixelCenter.getX()+(xResolution/2.),
                                                                            pixelCenter.getY()+(yResolution/2.)));

      std::unique_ptr<te::gm::Geometry> pixelEnvelopeAsGeometry(
            te::gm::GetGeomFromEnvelope(pixelEnvelope.get(), raster.getSRID()));

      if(pixelEnvelopeAsGeometry->intersects(polygonTemp.get()) == false)
        continue;

      double value;
      raster.getValue(col, row, value, band);

      if(value == noDataValue)
        continue;

      ++values[value];
    }
  }
}

void te::attributefill::GetPercentOfEachClassByArea(const te::rst::Raster &raster,
                                                    const unsigned int &band,
                                                    const te::gm::Polygon &polygon,
                                                    const uint32_t &minimumRow,
                                                    const uint32_t &minimumColumn,
                                                    const uint32_t &maximumRow,
                                                    const uint32_t &maximumColumn,
                                                    std::map<double, double> &percentOfEachClassByArea)
{
  assert(band < raster.getNumberOfBands());

  double polygonArea = polygon.getArea();

  std::unique_ptr<te::gm::Geometry> polygonTemp(
      dynamic_cast<te::gm::Geometry*>(polygon.clone()));

  if (polygonTemp->getSRID() != raster.getSRID())
    polygonTemp->transform(raster.getSRID());

  auto grid = raster.getGrid();

  auto xResolution = grid->getResolutionX();
  auto yResolution = grid->getResolutionY();

  for(uint32_t row = minimumRow; row <= maximumRow; ++row)
  {
    for(uint32_t col = minimumColumn; col <= maximumColumn; ++col)
    {
      auto pixelCenter = grid->gridToGeo(col, row);

      std::unique_ptr<te::gm::Envelope> pixelEnvelope( new te::gm::Envelope(pixelCenter.getX()-(xResolution/2.),
                                                                            pixelCenter.getY()-(yResolution/2.),
                                                                            pixelCenter.getX()+(xResolution/2.),
                                                                            pixelCenter.getY()+(yResolution/2.)));

      std::unique_ptr<te::gm::Geometry> pixelEnvelopeAsGeometry(
            te::gm::GetGeomFromEnvelope(pixelEnvelope.get(), raster.getSRID()));

      std::unique_ptr<te::gm::Geometry> pixelIntersectedGeometry(
            pixelEnvelopeAsGeometry->intersection(polygonTemp.get()));

      if(pixelIntersectedGeometry == nullptr)
        continue;

      double currentPixelArea;

      switch (pixelIntersectedGeometry->getGeomTypeId()) {
      case te::gm::PolygonType:
      {
        te::gm::Polygon polygon = *dynamic_cast<te::gm::Polygon*>(pixelIntersectedGeometry.get());
        currentPixelArea = polygon.getArea()/polygonArea;
      }
        break;
      case te::gm::MultiPolygonType:
      {
        te::gm::MultiPolygon multiPolygon = *dynamic_cast<te::gm::MultiPolygon*>(pixelIntersectedGeometry.get());
        currentPixelArea = multiPolygon.getArea()/polygonArea;
      }
        break;
      default:
        currentPixelArea = 0;
        break;
      }

      double value;
      raster.getValue(col, row, value, band);

      std::map<double, double>::iterator it = percentOfEachClassByArea.find(value);
      if(it != percentOfEachClassByArea.end())
        percentOfEachClassByArea[it->first] += currentPixelArea;
      else
        percentOfEachClassByArea.insert(std::pair<double, double>(value, currentPixelArea));
    }
  }
}