From f312a11ac96f3e28274c2de9086fe0d7310c7ac5 Mon Sep 17 00:00:00 2001
From: Dan Baston <dbaston@gmail.com>
Date: Tue, 21 Jan 2025 11:01:55 -0500
Subject: [PATCH] CLI: Check input CRS consistency

---
 src/exactextract.cpp        | 25 +++++++++++++++++++++++++
 src/gdal_raster_wrapper.cpp |  6 ++++++
 src/gdal_raster_wrapper.h   |  2 ++
 3 files changed, 33 insertions(+)

diff --git a/src/exactextract.cpp b/src/exactextract.cpp
index 81b0ec1..f3299c1 100644
--- 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
index 55add6d..b1a6944 100644
--- 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
index d6d6297..eef50b8 100644
--- 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: