Commit ead0df824a75586f4806ee2811009b67e25d85f6

Authored by Jordi Inglada
1 parent c1e9e695
Exists in master

Move utils to header files

CMakeLists.txt
... ... @@ -8,8 +8,8 @@ message(STATUS "Found OTB: ${OTB_USE_FILE}")
8 8 find_package(GSL REQUIRED)
9 9  
10 10  
11   -include_directories(${MASSP_SOURCE_DIR}/include )
12   -include_directories(/home/inglada/Dev/cesbio-utils/include)
  11 +include_directories(${MAASSP_SOURCE_DIR}/include )
  12 +include_directories(${MAASSP_SOURCE_DIR}/../cesbio-utils/include)
13 13  
14 14 add_subdirectory( src )
15 15  
... ...
include/maassp_utils.h 0 → 100644
... ... @@ -0,0 +1,104 @@
  1 +#ifndef MAASSP_UTILS_H
  2 +#define MAASSP_UTILS_H
  3 +
  4 +#include "ogr_utils.h"
  5 +#include <unordered_map>
  6 +
  7 +struct TileIntersection{
  8 + std::string tile_name;
  9 + double intersection_area;
  10 + otb::ogr::Feature tile_feature;
  11 +};
  12 +using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;
  13 +
  14 +DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer& polygon_source,
  15 + otb::ogr::DataSource::Pointer& tile_source)
  16 +{
  17 +
  18 + DeptToTiles result{};
  19 + auto polys = get_geometries(polygon_source);
  20 + auto tiles = get_geometries(tile_source);
  21 +
  22 + for(const auto& dept : polys)
  23 + {
  24 + auto code_dept = dept["CODE_DEPT"].GetValue<std::string>();
  25 + auto dept_geo = dept.GetGeometry();
  26 + for(const auto& tile : tiles)
  27 + {
  28 + auto code_tile = tile["Name"].GetValue<std::string>();
  29 + auto tile_geo = tile.GetGeometry();
  30 + if(dept_geo->Intersects(tile_geo))
  31 + {
  32 + auto ti = TileIntersection{code_tile,
  33 + get_area(dept_geo->Intersection(tile_geo)),
  34 + tile};
  35 + result[code_dept].emplace_back(ti);
  36 + }
  37 + }
  38 + }
  39 + return result;
  40 +}
  41 +
  42 +void print_dept_tiles(DeptToTiles tile_list, std::ostream& stream = std::cout)
  43 +{
  44 + for(const auto& d : tile_list)
  45 + {
  46 + stream << "Departement " << d.first << '\n';
  47 + for(const auto& t : d.second)
  48 + {
  49 + stream << "\t(" << t.tile_name << ", " << t.intersection_area << ")\n";
  50 + }
  51 + stream << '\n';
  52 + }
  53 +}
  54 +
  55 +
  56 +using AreaPerTile = std::unordered_map<std::string, double>;
  57 +using AreaPerTilePerDept = std::unordered_map<std::string, AreaPerTile>;
  58 +
  59 +
  60 +AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer& polygon_source,
  61 + std::vector<TileIntersection>& tiles,
  62 + std::vector<std::string> culture_patterns = {"_BLE_"})
  63 +{
  64 + auto polys = get_geometries(polygon_source);
  65 + AreaPerTile result;
  66 + for(const auto& parcelle : polys)
  67 + {
  68 + if(parcelle["LIBCULTURE"].HasBeenSet())
  69 + {
  70 + auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
  71 + if(cbutils::string::contains(code_cult, culture_patterns))
  72 + {
  73 + auto parcelle_geom = parcelle.GetGeometry();
  74 + for(const auto& t : tiles)
  75 + {
  76 + auto tile = t.tile_feature;
  77 + auto code_tile = tile["Name"].GetValue<std::string>();
  78 + auto tile_geo = tile.GetGeometry();
  79 + if(parcelle_geom->Within(tile_geo))
  80 + {
  81 + if(result.find(code_tile) == result.end())
  82 + { result[code_tile] = 0;}
  83 + result[code_tile] += get_area(parcelle_geom);
  84 + }
  85 + }
  86 + }
  87 + }
  88 + }
  89 + return result;
  90 +}
  91 +
  92 +void print_total_areas(AreaPerTilePerDept aptpd, std::ostream& stream = std::cout)
  93 +{
  94 + for(const auto& dept : aptpd)
  95 + {
  96 + auto code_dept = dept.first;
  97 + auto apt = dept.second;
  98 + stream << "Departement " << code_dept << '\n';
  99 + for(const auto& tile : apt)
  100 + stream << '\t' << tile.first << ": " << tile.second << '\n';
  101 + }
  102 +}
  103 +
  104 +#endif //MAASSP_UTILS_H
... ...
include/ogr_utils.h 0 → 100644
... ... @@ -0,0 +1,118 @@
  1 +#ifndef OGR_UTILS_H
  2 +#define OGR_UTILS_H
  3 +
  4 +#include "cbutils.h"
  5 +#include <iostream>
  6 +#include "otbOGRDataSourceWrapper.h"
  7 +#include "otbOGRFeatureWrapper.h"
  8 +#include "otbGeometriesProjectionFilter.h"
  9 +#include "otbGeometriesSet.h"
  10 +
  11 +using namespace cbutils::operators;
  12 +
  13 +struct Field
  14 +{
  15 + std::string name;
  16 + OGRFieldType type;
  17 +};
  18 +
  19 +std::vector<Field> get_fields(const otb::ogr::Layer& layer)
  20 +{
  21 + std::vector<Field> fields;
  22 + auto& layer_def = layer.GetLayerDefn();
  23 + for(auto i=0; i< layer_def.GetFieldCount(); ++i)
  24 + {
  25 + fields.push_back({layer_def.GetFieldDefn(i)->GetNameRef(),
  26 + layer_def.GetFieldDefn(i)->GetType()});
  27 + }
  28 + return fields;
  29 +}
  30 +
  31 +std::vector<Field> get_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0)
  32 +{
  33 + return get_fields(source->GetLayer(layer_id));
  34 +}
  35 +
  36 +bool has_field(const otb::ogr::Layer& layer, std::string field_name)
  37 +{
  38 + auto fields = get_fields(layer);
  39 + return find_if(cbegin(fields), cend(fields),
  40 + [&field_name](auto f){
  41 + return field_name == f.name;
  42 + })!=cend(fields);
  43 +}
  44 +
  45 +bool has_field(otb::ogr::DataSource::Pointer &source, std::string field_name,
  46 + int layer_id=0)
  47 +{
  48 + return has_field(source->GetLayer(layer_id), field_name);
  49 +}
  50 +
  51 +std::vector<std::string> get_field_names(const otb::ogr::Layer& layer)
  52 +{
  53 + std::vector<std::string> field_names;
  54 + for(const auto& f : get_fields(layer)) field_names.emplace_back(f.name);
  55 + return field_names;
  56 +}
  57 +
  58 +std::vector<std::string> get_field_names(otb::ogr::DataSource::Pointer &source, int layer_id=0)
  59 +{
  60 + return get_field_names(source->GetLayer(layer_id));
  61 +}
  62 +
  63 +void print_fields_in_layer(otb::ogr::Layer& layer, std::ostream& stream = std::cout)
  64 +{
  65 + for(const auto& f : get_fields(layer))
  66 + {
  67 + stream << f.name << '\t' << f.type << '\n';
  68 + }
  69 +}
  70 +
  71 +void print_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0,
  72 + std::ostream& stream = std::cout)
  73 +{
  74 + otb::ogr::Layer layer = source->GetLayer(layer_id);
  75 + print_fields_in_layer(layer, stream);
  76 +}
  77 +
  78 +std::vector<otb::ogr::Feature> get_geometries(otb::ogr::DataSource::Pointer& source)
  79 +{
  80 + std::vector<otb::ogr::Feature> result;
  81 + otb::ogr::Layer layer = source->GetLayer(0);
  82 + for(const auto& feature : layer)
  83 + {
  84 + result.emplace_back(feature);
  85 + }
  86 + return result;
  87 +}
  88 +
  89 +
  90 +otb::ogr::DataSource::Pointer reproject_geometries(otb::GeometriesSet::Pointer& source, std::string proj_ref)
  91 +{
  92 + auto reprojVector = otb::ogr::DataSource::New();
  93 + auto outputGeomSet = otb::GeometriesSet::New(reprojVector);
  94 + auto projection_filter = otb::GeometriesProjectionFilter::New();
  95 + projection_filter->SetInput(source);
  96 + projection_filter->SetOutput(outputGeomSet);
  97 + projection_filter->SetOutputProjectionRef(proj_ref);
  98 + projection_filter->Update();
  99 + return reprojVector;
  100 +}
  101 +
  102 +double get_area(const OGRGeometry* geom)
  103 +{
  104 + auto geo_name = std::string(geom->getGeometryName());
  105 + double area{0};
  106 + if(geo_name == "POLYGON")
  107 + {
  108 + area = dynamic_cast<const OGRPolygon*>(geom)->get_Area();
  109 + }
  110 + if(geo_name == "MULTIPOLYGON")
  111 + {
  112 + area = dynamic_cast<const OGRMultiPolygon*>(geom)->get_Area();
  113 + }
  114 + return area;
  115 +}
  116 +
  117 +
  118 +#endif //OGR_UTILS_H
... ...
src/select.cxx
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   -using namespace cbutils::operators;
10   -
11   -struct TileIntersection{
12   - std::string tile_name;
13   - double intersection_area;
14   - otb::ogr::Feature tile_feature;
15   -};
16   -using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;
17   -
18   -struct Field
19   -{
20   - std::string name;
21   - OGRFieldType type;
22   -};
23   -
24   -std::vector<Field> get_fields(const otb::ogr::Layer& layer)
25   -{
26   - std::vector<Field> fields;
27   - auto& layer_def = layer.GetLayerDefn();
28   - for(auto i=0; i< layer_def.GetFieldCount(); ++i)
29   - {
30   - fields.push_back({layer_def.GetFieldDefn(i)->GetNameRef(),
31   - layer_def.GetFieldDefn(i)->GetType()});
32   - }
33   - return fields;
34   -}
35   -
36   -std::vector<Field> get_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0)
37   -{
38   - return get_fields(source->GetLayer(layer_id));
39   -}
40   -
41   -bool has_field(const otb::ogr::Layer& layer, std::string field_name)
42   -{
43   - auto fields = get_fields(layer);
44   - return find_if(cbegin(fields), cend(fields),
45   - [&field_name](auto f){
46   - return field_name == f.name;
47   - })!=cend(fields);
48   -}
49   -
50   -bool has_field(otb::ogr::DataSource::Pointer &source, std::string field_name,
51   - int layer_id=0)
52   -{
53   - return has_field(source->GetLayer(layer_id), field_name);
54   -}
55   -
56   -std::vector<std::string> get_field_names(const otb::ogr::Layer& layer)
57   -{
58   - std::vector<std::string> field_names;
59   - for(const auto& f : get_fields(layer)) field_names.emplace_back(f.name);
60   - return field_names;
61   -}
62   -
63   -std::vector<std::string> get_field_names(otb::ogr::DataSource::Pointer &source, int layer_id=0)
64   -{
65   - return get_field_names(source->GetLayer(layer_id));
66   -}
67   -
68   -void print_fields_in_layer(otb::ogr::Layer& layer, std::ostream& stream = std::cout)
69   -{
70   - for(const auto& f : get_fields(layer))
71   - {
72   - stream << f.name << '\t' << f.type << '\n';
73   - }
74   -}
75   -
76   -void print_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0,
77   - std::ostream& stream = std::cout)
78   -{
79   - otb::ogr::Layer layer = source->GetLayer(layer_id);
80   - print_fields_in_layer(layer, stream);
81   -}
82   -
83   -std::vector<otb::ogr::Feature> get_geometries(otb::ogr::DataSource::Pointer& source)
84   -{
85   - std::vector<otb::ogr::Feature> result;
86   - otb::ogr::Layer layer = source->GetLayer(0);
87   - for(const auto& feature : layer)
88   - {
89   - result.emplace_back(feature);
90   - }
91   - return result;
92   -}
93   -
94   -
95   -otb::ogr::DataSource::Pointer reproject_geometries(otb::GeometriesSet::Pointer& source, std::string proj_ref)
96   -{
97   - auto reprojVector = otb::ogr::DataSource::New();
98   - auto outputGeomSet = otb::GeometriesSet::New(reprojVector);
99   - auto projection_filter = otb::GeometriesProjectionFilter::New();
100   - projection_filter->SetInput(source);
101   - projection_filter->SetOutput(outputGeomSet);
102   - projection_filter->SetOutputProjectionRef(proj_ref);
103   - projection_filter->Update();
104   - return reprojVector;
105   -}
106   -
107   -double get_area(const OGRGeometry* geom)
108   -{
109   - auto geo_name = std::string(geom->getGeometryName());
110   - double area{0};
111   - if(geo_name == "POLYGON")
112   - {
113   - area = dynamic_cast<const OGRPolygon*>(geom)->get_Area();
114   - }
115   - if(geo_name == "MULTIPOLYGON")
116   - {
117   - area = dynamic_cast<const OGRMultiPolygon*>(geom)->get_Area();
118   - }
119   - return area;
120   -}
121   -
122   -DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer& polygon_source,
123   - otb::ogr::DataSource::Pointer& tile_source)
124   -{
125   -
126   - DeptToTiles result{};
127   - auto polys = get_geometries(polygon_source);
128   - auto tiles = get_geometries(tile_source);
129   -
130   - for(const auto& dept : polys)
131   - {
132   - auto code_dept = dept["CODE_DEPT"].GetValue<std::string>();
133   - auto dept_geo = dept.GetGeometry();
134   - for(const auto& tile : tiles)
135   - {
136   - auto code_tile = tile["Name"].GetValue<std::string>();
137   - auto tile_geo = tile.GetGeometry();
138   - if(dept_geo->Intersects(tile_geo))
139   - {
140   - auto ti = TileIntersection{code_tile,
141   - get_area(dept_geo->Intersection(tile_geo)),
142   - tile};
143   - result[code_dept].emplace_back(ti);
144   - }
145   - }
146   - }
147   - return result;
148   -}
149   -
150   -void print_dept_tiles(DeptToTiles tile_list, std::ostream& stream = std::cout)
151   -{
152   - for(const auto& d : tile_list)
153   - {
154   - stream << "Departement " << d.first << '\n';
155   - for(const auto& t : d.second)
156   - {
157   - stream << "\t(" << t.tile_name << ", " << t.intersection_area << ")\n";
158   - }
159   - stream << '\n';
160   - }
161   -}
162   -
163   -
164   -using AreaPerTile = std::unordered_map<std::string, double>;
165   -using AreaPerTilePerDept = std::unordered_map<std::string, AreaPerTile>;
166   -
167   -
168   -AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer& polygon_source,
169   - std::vector<TileIntersection>& tiles,
170   - std::vector<std::string> culture_patterns = {"_BLE_"})
171   -{
172   - auto polys = get_geometries(polygon_source);
173   - AreaPerTile result;
174   - for(const auto& parcelle : polys)
175   - {
176   - if(parcelle["LIBCULTURE"].HasBeenSet())
177   - {
178   - auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
179   - if(cbutils::string::contains(code_cult, culture_patterns))
180   - {
181   - auto parcelle_geom = parcelle.GetGeometry();
182   - for(const auto& t : tiles)
183   - {
184   - auto tile = t.tile_feature;
185   - auto code_tile = tile["Name"].GetValue<std::string>();
186   - auto tile_geo = tile.GetGeometry();
187   - if(parcelle_geom->Within(tile_geo))
188   - {
189   - if(result.find(code_tile) == result.end())
190   - { result[code_tile] = 0;}
191   - result[code_tile] += get_area(parcelle_geom);
192   - }
193   - }
194   - }
195   - }
196   - }
197   - return result;
198   -}
199   -
200   -void print_total_areas(AreaPerTilePerDept aptpd, std::ostream& stream = std::cout)
201   -{
202   - for(const auto& dept : aptpd)
203   - {
204   - auto code_dept = dept.first;
205   - auto apt = dept.second;
206   - stream << "Departement " << code_dept << '\n';
207   - for(const auto& tile : apt)
208   - stream << '\t' << tile.first << ": " << tile.second << '\n';
209   - }
210   -}
  1 +#include "maassp_utils.h"
211 2  
212 3 int main()
213 4 {
... ...