Commit ed7b84a3d9dedb83007b5fea7cedacf15da627f5

Authored by Jordi Inglada
0 parents
Exists in master

Compute intersections and areas between depts and S2 tiles

CMakeLists.txt 0 → 100644
  1 +++ a/CMakeLists.txt
... ... @@ -0,0 +1,20 @@
  1 +cmake_minimum_required(VERSION 3.0)
  2 +project(MAASSP)
  3 +
  4 +find_package(OTB REQUIRED)
  5 +include(${OTB_USE_FILE})
  6 +message(STATUS "Found OTB: ${OTB_USE_FILE}")
  7 +
  8 +find_package(GSL REQUIRED)
  9 +
  10 +
  11 +include_directories(${MASSP_SOURCE_DIR}/include )
  12 +include_directories(/home/inglada/Dev/cesbio-utils/include)
  13 +
  14 +add_subdirectory( src )
  15 +
  16 +include(CTest)
  17 +if(BUILD_TESTING)
  18 + enable_testing()
  19 + add_subdirectory( tests )
  20 +endif()
... ...
README.org 0 → 100644
  1 +++ a/README.org
... ... @@ -0,0 +1,7 @@
  1 +* Build, install
  2 +#+BEGIN_SRC bash
  3 +mkdir build
  4 +cd build
  5 +cmake -DCMAKE_CXX_FLAGS:STRING=-std=c++17 -DOTB_DIR:PATH=/home/inglada/OTB/SuperBuild/lib/cmake/OTB-6.7 ../
  6 +#+END_SRC
  7 +
... ...
src/CMakeLists.txt 0 → 100644
  1 +++ a/src/CMakeLists.txt
... ... @@ -0,0 +1,3 @@
  1 +
  2 +add_executable(select select.cxx)
  3 +target_link_libraries(select stdc++fs ${OTB_LIBRARIES})
... ...
src/select.cxx 0 → 100644
  1 +++ a/src/select.cxx
... ... @@ -0,0 +1,135 @@
  1 +#include "cbutils.h"
  2 +#include <iostream>
  3 +#include "otbOGRDataSourceWrapper.h"
  4 +#include "otbOGRFeatureWrapper.h"
  5 +#include <unordered_map>
  6 +#include "otbGeometriesProjectionFilter.h"
  7 +#include "otbGeometriesSet.h"
  8 +
  9 +struct TileIntersection{
  10 + std::string tile_name;
  11 + double intersection_area;
  12 +};
  13 +using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;
  14 +
  15 +void print_fields_in_layer(otb::ogr::Layer& layer)
  16 +{
  17 + auto& layer_def = layer.GetLayerDefn();
  18 + for(auto i=0; i< layer_def.GetFieldCount(); ++i)
  19 + std::cout << layer_def.GetFieldDefn(i)->GetNameRef() <<
  20 + '\t' << layer_def.GetFieldDefn(i)->GetType() << '\n';
  21 +
  22 +}
  23 +
  24 +void print_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0)
  25 +{
  26 +
  27 + otb::ogr::Layer layer = source->GetLayer(layer_id);
  28 + print_fields_in_layer(layer);
  29 +}
  30 +
  31 +std::vector<otb::ogr::Feature> get_geometries(otb::ogr::DataSource::Pointer& source)
  32 +{
  33 + std::vector<otb::ogr::Feature> result;
  34 + otb::ogr::Layer layer = source->GetLayer(0);
  35 + for(const auto& feature : layer)
  36 + {
  37 + result.emplace_back(feature);
  38 + }
  39 + return result;
  40 +}
  41 +
  42 +
  43 +otb::ogr::DataSource::Pointer reproject_geometries(otb::GeometriesSet::Pointer& source, std::string proj_ref)
  44 +{
  45 + auto reprojVector = otb::ogr::DataSource::New();
  46 + auto outputGeomSet = otb::GeometriesSet::New(reprojVector);
  47 + auto projection_filter = otb::GeometriesProjectionFilter::New();
  48 + projection_filter->SetInput(source);
  49 + projection_filter->SetOutput(outputGeomSet);
  50 + projection_filter->SetOutputProjectionRef(proj_ref);
  51 + projection_filter->Update();
  52 + return reprojVector;
  53 +}
  54 +
  55 +double get_area(const OGRGeometry* geom)
  56 +{
  57 + auto geo_name = std::string(geom->getGeometryName());
  58 + double area{0};
  59 + if(geo_name == "POLYGON")
  60 + {
  61 + area = dynamic_cast<const OGRPolygon*>(geom)->get_Area();
  62 + }
  63 + if(geo_name == "MULTIPOLYGON")
  64 + {
  65 + area = dynamic_cast<const OGRMultiPolygon*>(geom)->get_Area();
  66 + }
  67 + return area;
  68 +}
  69 +
  70 +DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer& polygon_source,
  71 + otb::ogr::DataSource::Pointer& tile_source)
  72 +{
  73 +
  74 + DeptToTiles result{};
  75 + auto polys = get_geometries(polygon_source);
  76 + auto tiles = get_geometries(tile_source);
  77 +
  78 + for(const auto& dept : polys)
  79 + {
  80 + auto code_dept = dept["CODE_DEPT"].GetValue<std::string>();
  81 + auto dept_geo = dept.GetGeometry();
  82 + for(const auto& tile : tiles)
  83 + {
  84 + auto code_tile = tile["Name"].GetValue<std::string>();
  85 + auto tile_geo = tile.GetGeometry();
  86 + if(dept_geo->Intersects(tile_geo))
  87 + {
  88 + auto ti = TileIntersection{code_tile,
  89 + get_area(dept_geo->Intersection(tile_geo))};
  90 + result[code_dept].emplace_back(ti);
  91 + }
  92 + }
  93 + }
  94 + return result;
  95 +}
  96 +
  97 +int main()
  98 +{
  99 + std::string s2_tiles_shp{"/home/inglada/stok/DATA/iota2/S22017/Features.shp"};
  100 + std::string departements_shp{"/home/inglada/stok/DATA/departements_shp/"
  101 + "departements_lambert93.shp"};
  102 +
  103 + otb::ogr::DataSource::Pointer s2_tiles = otb::ogr::DataSource::New(s2_tiles_shp);
  104 +
  105 + otb::ogr::DataSource::Pointer departements = otb::ogr::DataSource::New(departements_shp);
  106 +
  107 + print_fields(s2_tiles);
  108 + auto dept_geometries = otb::GeometriesSet::New(departements);
  109 + auto depts = reproject_geometries(dept_geometries,
  110 + s2_tiles->GetLayer(0).GetProjectionRef());
  111 + auto tile_list = get_intersecting_tiles(depts, s2_tiles);
  112 +
  113 + for(const auto& d : tile_list)
  114 + {
  115 + std::cout << "Departement " << d.first << '\n';
  116 + for(const auto& t : d.second)
  117 + {
  118 + std::cout << "\t(" << t.tile_name << ", " << t.intersection_area << ")\n";
  119 + }
  120 + std::cout << '\n';
  121 + }
  122 + /* std::string shape_dir{"/home/inglada/stok/DATA/OSO/MAA_SSP/"};
  123 + for(const auto& f : cbutils::file::list_files(shape_dir, ".*RPG_TERLAB_DEP88-89_2017.7z"))
  124 + {
  125 + std::cout << "Unzipping " << f << '\n';
  126 + cbutils::system::call("7z x "+f+" -aoa");
  127 + for(const auto& shape : cbutils::file::list_files("./",".*shp"))
  128 + {
  129 +
  130 + }
  131 + }*/
  132 +
  133 + return EXIT_SUCCESS;
  134 +}
  135 +
... ...