Commit a41f6fdebdc6589c9fb9d94ee3022db462a82464

Authored by Jordi Inglada
1 parent ed7b84a3
Exists in master

Compute area per cult per dept per tile

Showing 1 changed file with 78 additions and 20 deletions   Show diff stats
src/select.cxx
... ... @@ -9,6 +9,7 @@
9 9 struct TileIntersection{
10 10 std::string tile_name;
11 11 double intersection_area;
  12 + otb::ogr::Feature tile_feature;
12 13 };
13 14 using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;
14 15  
... ... @@ -86,7 +87,8 @@ DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer&amp; polygon_source
86 87 if(dept_geo->Intersects(tile_geo))
87 88 {
88 89 auto ti = TileIntersection{code_tile,
89   - get_area(dept_geo->Intersection(tile_geo))};
  90 + get_area(dept_geo->Intersection(tile_geo)),
  91 + tile};
90 92 result[code_dept].emplace_back(ti);
91 93 }
92 94 }
... ... @@ -94,42 +96,98 @@ DeptToTiles get_intersecting_tiles(otb::ogr::DataSource::Pointer&amp; polygon_source
94 96 return result;
95 97 }
96 98  
  99 +void print_dept_tiles(DeptToTiles tile_list)
  100 +{
  101 + for(const auto& d : tile_list)
  102 + {
  103 + std::cout << "Departement " << d.first << '\n';
  104 + for(const auto& t : d.second)
  105 + {
  106 + std::cout << "\t(" << t.tile_name << ", " << t.intersection_area << ")\n";
  107 + }
  108 + std::cout << '\n';
  109 + }
  110 +}
  111 +
  112 +bool matches_pattern(std::string s, std::vector<std::string> patterns)
  113 +{
  114 + bool result{false};
  115 + for(const auto& pat : patterns)
  116 + {
  117 + if(s.find(pat) != std::string::npos)
  118 + {
  119 + return true;
  120 + }
  121 + }
  122 + return result;
  123 +}
  124 +
  125 +using AreaPerTile = std::unordered_map<std::string, double>;
  126 +
  127 +AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer& polygon_source,
  128 + std::vector<TileIntersection>& tiles,
  129 + std::vector<std::string> culture_patterns = {"_BLE_"})
  130 +{
  131 + auto polys = get_geometries(polygon_source);
  132 +
  133 + AreaPerTile result;
  134 + for(const auto& parcelle : polys)
  135 + {
  136 + auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
  137 + if(matches_pattern(code_cult, culture_patterns))
  138 + {
  139 + auto parcelle_geom = parcelle.GetGeometry();
  140 + for(const auto& t : tiles)
  141 + {
  142 + auto tile = t.tile_feature;
  143 + auto code_tile = tile["Name"].GetValue<std::string>();
  144 + auto tile_geo = tile.GetGeometry();
  145 + if(parcelle_geom->Within(tile_geo))
  146 + {
  147 + if(result.find(code_tile) == result.end())
  148 + { result[code_tile] = 0;}
  149 + result[code_tile] += get_area(parcelle_geom);
  150 + }
  151 + }
  152 + }
  153 + }
  154 + return result;
  155 +}
  156 +
97 157 int main()
98 158 {
99 159 std::string s2_tiles_shp{"/home/inglada/stok/DATA/iota2/S22017/Features.shp"};
100 160 std::string departements_shp{"/home/inglada/stok/DATA/departements_shp/"
101 161 "departements_lambert93.shp"};
102 162  
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);
  163 + auto s2_tiles = otb::ogr::DataSource::New(s2_tiles_shp);
  164 + auto departements = otb::ogr::DataSource::New(departements_shp);
108 165 auto dept_geometries = otb::GeometriesSet::New(departements);
109 166 auto depts = reproject_geometries(dept_geometries,
110 167 s2_tiles->GetLayer(0).GetProjectionRef());
111 168 auto tile_list = get_intersecting_tiles(depts, s2_tiles);
  169 + print_dept_tiles(tile_list);
112 170  
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/"};
  171 + std::string shaphe_file{"/home/inglada/Dev/MAA-SSP/build/SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_089_20180210.shp"};
  172 + auto parcelles = otb::ogr::DataSource::New(shaphe_file);
  173 + auto parc_geometries = otb::GeometriesSet::New(parcelles);
  174 + auto parcs = reproject_geometries(parc_geometries,
  175 + s2_tiles->GetLayer(0).GetProjectionRef());
  176 +
  177 + auto apt = compute_area_per_tile(parcs, tile_list["89"], {"TOURNESOL", "COLZA"});
  178 + for(const auto& x : apt) std::cout << x.first << "\t" << x.second << '\n';
  179 +
  180 +/*
  181 + std::string shape_dir{"/home/inglada/stok/DATA/OSO/MAA_SSP/"};
123 182 for(const auto& f : cbutils::file::list_files(shape_dir, ".*RPG_TERLAB_DEP88-89_2017.7z"))
124 183 {
125 184 std::cout << "Unzipping " << f << '\n';
126 185 cbutils::system::call("7z x "+f+" -aoa");
127   - for(const auto& shape : cbutils::file::list_files("./",".*shp"))
  186 + for(const auto& shape : cbutils::file::list_files("/home/inglada/Dev/MAA-SSP/build/",".*shp"))
128 187 {
129   -
130 188 }
131   - }*/
132   -
  189 + }
  190 +*/
133 191 return EXIT_SUCCESS;
134 192 }
135 193  
... ...