CLI: Check input CRS consistency

isciences · Jan 21, 2025 · f312a11 · f312a11
1 parent f9926f1
commit f312a11
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 0 deletions.
diff --git a/src/exactextract.cpp b/src/exactextract.cpp
@@ -32,6 +32,8 @@
 #include "utils_cli.h"
 #include "version.h"
 
+#include <ogr_srs_api.h>
+
 using exactextract::GDALDatasetWrapper;
 using exactextract::GDALRasterWrapper;
 using exactextract::Operation;
@@ -43,6 +45,9 @@ load_dataset(const std::string& descriptor,
              const std::string& dst_id_name,
              const std::string& dst_id_type);
 
+static void
+check_crs_consistent(const GDALDatasetWrapper& features, const std::vector<std::unique_ptr<exactextract::RasterSource>>& rasters, const std::vector<std::unique_ptr<exactextract::RasterSource>>& weights);
+
 int
 main(int argc, char** argv)
 {
@@ -114,6 +119,8 @@ main(int argc, char** argv)
 
         GDALDatasetWrapper shp = load_dataset(poly_descriptor, include_cols, src_id_name, dst_id_name, dst_id_type);
 
+        check_crs_consistent(shp, rasters, weights);
+
         std::unique_ptr<exactextract::GDALWriter> gdal_writer = std::make_unique<exactextract::GDALWriter>(
           output_filename, !nested_output, shp.srs());
 
@@ -216,3 +223,21 @@ load_dataset(const std::string& descriptor,
 
     return ds;
 }
+
+static void
+check_crs_consistent(const GDALDatasetWrapper& features, const std::vector<std::unique_ptr<exactextract::RasterSource>>& rasters, const std::vector<std::unique_ptr<exactextract::RasterSource>>& weights)
+{
+    OGRSpatialReferenceH expected = features.srs();
+    for (const auto& raster : rasters) {
+        OGRSpatialReferenceH srs = static_cast<const GDALRasterWrapper*>(raster.get())->srs();
+        if (!OSRIsSame(srs, expected)) {
+            std::cerr << "Input features do not have the same CRS as raster " << raster->name() << std::endl;
+        }
+    }
+    for (const auto& raster : weights) {
+        OGRSpatialReferenceH srs = static_cast<const GDALRasterWrapper*>(raster.get())->srs();
+        if (!OSRIsSame(srs, expected)) {
+            std::cerr << "Input features do not have the same CRS as raster " << raster->name() << std::endl;
+        }
+    }
+}
diff --git a/src/gdal_raster_wrapper.cpp b/src/gdal_raster_wrapper.cpp
@@ -242,4 +242,10 @@ GDALRasterWrapper::compute_raster_grid()
     m_grid = { box, dx, dy };
 }
 
+OGRSpatialReferenceH
+GDALRasterWrapper::srs() const
+{
+    return GDALGetSpatialRef(m_rast->get());
+}
+
 }
diff --git a/src/gdal_raster_wrapper.h b/src/gdal_raster_wrapper.h
@@ -66,6 +66,8 @@ class GDALRasterWrapper : public RasterSource
 
     ~GDALRasterWrapper() override;
 
+    OGRSpatialReferenceH srs() const;
+
     bool cartesian() const;
 
   private: