Commit c1e9e6955cc7ca6443e817397bde48ba5f0ee6ad

Authored by Jordi Inglada
1 parent 7234f4d5
Exists in master

Check if the dept has rendnorme info

Showing 1 changed file with 103 additions and 29 deletions   Show diff stats
src/select.cxx
... ... @@ -15,19 +15,67 @@ struct TileIntersection{
15 15 };
16 16 using DeptToTiles = std::unordered_map<std::string, std::vector<TileIntersection>>;
17 17  
  18 +struct Field
  19 +{
  20 + std::string name;
  21 + OGRFieldType type;
  22 +};
18 23  
19   -void print_fields_in_layer(otb::ogr::Layer& layer, std::ostream& stream = std::cout)
  24 +std::vector<Field> get_fields(const otb::ogr::Layer& layer)
20 25 {
  26 + std::vector<Field> fields;
21 27 auto& layer_def = layer.GetLayerDefn();
22 28 for(auto i=0; i< layer_def.GetFieldCount(); ++i)
23   - stream << layer_def.GetFieldDefn(i)->GetNameRef() <<
24   - '\t' << layer_def.GetFieldDefn(i)->GetType() << '\n';
  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 + }
25 74 }
26 75  
27 76 void print_fields(otb::ogr::DataSource::Pointer &source, int layer_id=0,
28 77 std::ostream& stream = std::cout)
29 78 {
30   -
31 79 otb::ogr::Layer layer = source->GetLayer(layer_id);
32 80 print_fields_in_layer(layer, stream);
33 81 }
... ... @@ -114,6 +162,8 @@ void print_dept_tiles(DeptToTiles tile_list, std::ostream&amp; stream = std::cout)
114 162  
115 163  
116 164 using AreaPerTile = std::unordered_map<std::string, double>;
  165 +using AreaPerTilePerDept = std::unordered_map<std::string, AreaPerTile>;
  166 +
117 167  
118 168 AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer& polygon_source,
119 169 std::vector<TileIntersection>& tiles,
... ... @@ -123,27 +173,42 @@ AreaPerTile compute_area_per_tile(otb::ogr::DataSource::Pointer&amp; polygon_source,
123 173 AreaPerTile result;
124 174 for(const auto& parcelle : polys)
125 175 {
126   - auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
127   - if(cbutils::string::contains(code_cult, culture_patterns))
  176 + if(parcelle["LIBCULTURE"].HasBeenSet())
128 177 {
129   - auto parcelle_geom = parcelle.GetGeometry();
130   - for(const auto& t : tiles)
  178 + auto code_cult = parcelle["LIBCULTURE"].GetValue<std::string>();
  179 + if(cbutils::string::contains(code_cult, culture_patterns))
131 180 {
132   - auto tile = t.tile_feature;
133   - auto code_tile = tile["Name"].GetValue<std::string>();
134   - auto tile_geo = tile.GetGeometry();
135   - if(parcelle_geom->Within(tile_geo))
  181 + auto parcelle_geom = parcelle.GetGeometry();
  182 + for(const auto& t : tiles)
136 183 {
137   - if(result.find(code_tile) == result.end())
138   - { result[code_tile] = 0;}
139   - result[code_tile] += get_area(parcelle_geom);
  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 + }
140 193 }
141   - }
  194 + }
142 195 }
143 196 }
144 197 return result;
145 198 }
146 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 +}
  211 +
147 212 int main()
148 213 {
149 214 std::string s2_tiles_shp{"/home/inglada/stok/DATA/iota2/S22017/Features.shp"};
... ... @@ -159,27 +224,36 @@ int main()
159 224 print_dept_tiles(tile_list);
160 225  
161 226 std::string shape_dir{"/home/inglada/stok/DATA/OSO/MAA_SSP/"};
  227 + AreaPerTilePerDept aptpd{};
162 228 for(const auto& f : cbutils::file::list_files(shape_dir,
163   - ".*RPG_TERLAB_DEP88-89_2017.7z"))
  229 + ".*RPG_TERLAB_DEP0.*_2017.7z"))
164 230 {
165 231 std::cout << "Unzipping " << f << '\n';
166 232 cbutils::system::call("7z x "+f+" -aoa");
167   - for(auto& shape : cbutils::file::list_files("/home/inglada/Dev/MAA-SSP/build/",
168   - ".*shp"))
  233 + auto shapes_to_process = cbutils::file::list_files("./",
  234 + ".*shp");
  235 + for(auto& shape : shapes_to_process)
169 236 {
170 237 auto parcelles = otb::ogr::DataSource::New(shape);
171   - auto parc_geometries = otb::GeometriesSet::New(parcelles);
172   - auto parcs = reproject_geometries(parc_geometries,
173   - s2_tiles->GetLayer(0).GetProjectionRef());
174   -
175   - auto code_dept = cbutils::string::split(shape, "_")[1].substr(1,2);
176   - std::cout << "Departement " << code_dept << '\n';
177   - auto apt = compute_area_per_tile(parcs, tile_list[code_dept],
178   - {"TOURNESOL", "COLZA"});
179   - for(const auto& x : apt) std::cout << x.first << "\t" << x.second << '\n';
  238 + if(has_field(parcelles, "RENDNORME"))
  239 + {
  240 + auto parc_geometries = otb::GeometriesSet::New(parcelles);
  241 + auto parcs = reproject_geometries(parc_geometries,
  242 + s2_tiles->GetLayer(0).GetProjectionRef());
  243 +
  244 + auto code_dept = cbutils::string::split(shape, "_")[1].substr(1,2);
  245 + std::cout << "Departement " << code_dept << '\n';
  246 + auto apt = compute_area_per_tile(parcs, tile_list[code_dept],
  247 + {"TOURNESOL", "COLZA"});
  248 + aptpd[code_dept] = apt;
  249 + }
  250 + auto files_to_remove = std::string("."+cbutils::string::split(shape, ".")[1]+".*");
  251 + std::cout << "Remove " << files_to_remove << '\n';
  252 + cbutils::system::call("rm -f "+files_to_remove);
180 253 }
181 254 }
182   -
  255 + std::ofstream outfile("/tmp/stats.txt");
  256 + print_total_areas(aptpd, outfile);
183 257  
184 258 return EXIT_SUCCESS;
185 259 }
... ...