select.cxx 5.89 KB
#include "cbutils.h"
#include <iostream>
#include "otbOGRDataSourceWrapper.h"
#include "otbOGRFeatureWrapper.h"
#include <unordered_map>
#include "otbGeometriesProjectionFilter.h"
#include "otbGeometriesSet.h"

struct TileIntersection{
  std::string tile_name;
  double intersection_area;
  otb::ogr::Feature tile_feature;
};
using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;

void print_fields_in_layer(otb::ogr::Layer& layer)
{
  auto& layer_def = layer.GetLayerDefn();
  for(auto i=0; i< layer_def.GetFieldCount(); ++i)
    std::cout << layer_def.GetFieldDefn(i)->GetNameRef() <<
      '\t' <<   layer_def.GetFieldDefn(i)->GetType() << '\n';

}

void print_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0)
{

  otb::ogr::Layer layer = source->GetLayer(layer_id);
  print_fields_in_layer(layer);
}

std::vector<otb::ogr::Feature> get_geometries(otb::ogr::DataSource::Pointer& source)
{
  std::vector<otb::ogr::Feature> result;
  otb::ogr::Layer layer = source->GetLayer(0);
  for(const auto& feature : layer)
    {
    result.emplace_back(feature);
    }
  return result;
}


otb::ogr::DataSource::Pointer reproject_geometries(otb::GeometriesSet::Pointer& source, std::string proj_ref)
{
  auto reprojVector = otb::ogr::DataSource::New();
  auto outputGeomSet = otb::GeometriesSet::New(reprojVector);
  auto projection_filter = otb::GeometriesProjectionFilter::New();
  projection_filter->SetInput(source);
  projection_filter->SetOutput(outputGeomSet);
  projection_filter->SetOutputProjectionRef(proj_ref);
  projection_filter->Update();
  return reprojVector;
}

double get_area(const OGRGeometry* geom)
{
  auto geo_name = std::string(geom->getGeometryName());
  double area{0};
  if(geo_name == "POLYGON")
    {
    area = dynamic_cast<const OGRPolygon*>(geom)->get_Area();
    }
  if(geo_name == "MULTIPOLYGON")
    {
    area = dynamic_cast<const OGRMultiPolygon*>(geom)->get_Area();
    }
  return area;
}

DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer& polygon_source,
                                   otb::ogr::DataSource::Pointer& tile_source)
{
  
  DeptToTiles result{};
  auto polys = get_geometries(polygon_source);
  auto tiles = get_geometries(tile_source);

  for(const auto& dept : polys)
    {
    auto code_dept = dept["CODE_DEPT"].GetValue<std::string>();
    auto dept_geo = dept.GetGeometry();
    for(const auto& tile : tiles)
      {
      auto code_tile = tile["Name"].GetValue<std::string>();
      auto tile_geo = tile.GetGeometry();      
      if(dept_geo->Intersects(tile_geo))
        {
        auto ti = TileIntersection{code_tile, 
                                   get_area(dept_geo->Intersection(tile_geo)),
                                   tile};
        result[code_dept].emplace_back(ti);
        }
      }
    }
  return result;
}

void print_dept_tiles(DeptToTiles tile_list)
{
  for(const auto& d : tile_list)
    {
    std::cout << "Departement " << d.first << '\n';
    for(const auto& t : d.second)
      {
      std::cout << "\t(" << t.tile_name << ", " << t.intersection_area << ")\n";
      }
    std::cout << '\n';
    }
}

bool matches_pattern(std::string s, std::vector<std::string> patterns)
{
  bool result{false};
  for(const auto& pat : patterns)
    {
    if(s.find(pat) != std::string::npos)
      {
      return true;
      }
    }
  return result;
}

using AreaPerTile = std::unordered_map<std::string, double>;

AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer& polygon_source,
                                  std::vector<TileIntersection>& tiles,
                                  std::vector<std::string> culture_patterns = {"_BLE_"})
{
  auto polys = get_geometries(polygon_source);

  AreaPerTile result;
  for(const auto& parcelle : polys)
    {
    auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
    if(matches_pattern(code_cult, culture_patterns))
      {
      auto parcelle_geom = parcelle.GetGeometry();
      for(const auto& t : tiles)
        {
        auto tile = t.tile_feature;
        auto code_tile = tile["Name"].GetValue<std::string>();
        auto tile_geo = tile.GetGeometry();      
        if(parcelle_geom->Within(tile_geo))
          {
          if(result.find(code_tile) == result.end()) 
            { result[code_tile] = 0;}
          result[code_tile] += get_area(parcelle_geom);
          }
         }
      }
    }
  return result;
}

int main()
{
  std::string s2_tiles_shp{"/home/inglada/stok/DATA/iota2/S22017/Features.shp"};
  std::string departements_shp{"/home/inglada/stok/DATA/departements_shp/"
                               "departements_lambert93.shp"};

  auto s2_tiles = otb::ogr::DataSource::New(s2_tiles_shp);
  auto departements = otb::ogr::DataSource::New(departements_shp);
  auto dept_geometries = otb::GeometriesSet::New(departements);
  auto depts = reproject_geometries(dept_geometries,
                                    s2_tiles->GetLayer(0).GetProjectionRef());
  auto tile_list = get_intersecting_tiles(depts, s2_tiles);
  print_dept_tiles(tile_list);

  std::string shaphe_file{"/home/inglada/Dev/MAA-SSP/build/SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_089_20180210.shp"};
  auto parcelles = otb::ogr::DataSource::New(shaphe_file);
  auto parc_geometries = otb::GeometriesSet::New(parcelles);
  auto parcs = reproject_geometries(parc_geometries,
                                    s2_tiles->GetLayer(0).GetProjectionRef());

  auto apt = compute_area_per_tile(parcs, tile_list["89"], {"TOURNESOL", "COLZA"});
  for(const auto& x : apt) std::cout << x.first << "\t" << x.second << '\n';

/*
  std::string shape_dir{"/home/inglada/stok/DATA/OSO/MAA_SSP/"};
  for(const auto& f : cbutils::file::list_files(shape_dir, ".*RPG_TERLAB_DEP88-89_2017.7z"))
    {
    std::cout << "Unzipping " << f << '\n';
    cbutils::system::call("7z x "+f+" -aoa");
    for(const auto& shape : cbutils::file::list_files("/home/inglada/Dev/MAA-SSP/build/",".*shp"))
      {
      }
    }
*/
  return EXIT_SUCCESS;
}