Ensure geojson outputs are correctly georeferenced

Author: Chrismarsh

src/core.cpp

@@ -1601,11 +1601,6 @@ void core::init(int argc, char **argv)
         }
     }
 
-    // needs to go here as both mesh needs to be loaded and the output dir needs to be known
-    boost::filesystem::create_directories(output_folder_path / "mesh_boundingbox");
-    auto mesh_boundingbox_path = output_folder_path / "mesh_boundingbox" / std::format("mesh_bbox_{}.geojson", _comm_world.rank());
-    _mesh->write_bbox_geojson(mesh_boundingbox_path.string());
-
 
     pt::json_parser::write_json((output_folder_path / "config.json" ).string(),cfg); // output a full dump of the cfg, after all modifications, to the output directory
     _cfg = cfg;

src/mesh/triangulation.cpp

@@ -1451,7 +1451,7 @@ void triangulation::reorder_faces(std::vector<size_t> permutation)
 
 void triangulation::write_bbox_geojson(const std::string& filename)
 {
-    gis::bbox2geojson(_bounding_box.x_min, _bounding_box.y_min, _bounding_box.x_max, _bounding_box.y_max, filename);
+    gis::bbox2geojson(_bounding_box.x_min, _bounding_box.y_min, _bounding_box.x_max, _bounding_box.y_max, filename, proj4());
 }
 
 void triangulation::load_partition_from_mesh(const std::string& mesh_filename)

src/utility/gis.cpp

@@ -3,7 +3,7 @@
 namespace gis
 {
 
-    void bbox2geojson(double xmin, double ymin, double xmax, double ymax, const std::string& filepath)
+    void bbox2geojson(double xmin, double ymin, double xmax, double ymax, const std::string& filepath, const std::string& proj4)
     {
         GDALAllRegister();
         const char *pszDriverName = "GeoJSON";
@@ -19,12 +19,22 @@ namespace gis
             CHM_THROW_EXCEPTION(chm_error, "Creation of output file failed");
         }
 
-        OGRLayer *poLayer = poDS->CreateLayer("bbox", nullptr, wkbPolygon, nullptr);
+        OGRSpatialReference oSRS;
+
+        if (oSRS.importFromProj4(proj4.c_str()) != OGRERR_NONE)
+        {
+            CHM_THROW_EXCEPTION(chm_error, "Failed to import PROJ.4 string");
+        }
+
+        char **papszOptions = nullptr;
+        papszOptions = CSLAddNameValue(papszOptions, "RFC7946", "TRUE");
+        OGRLayer *poLayer = poDS->CreateLayer("bbox", &oSRS, wkbPolygon, papszOptions);
         if (poLayer == nullptr)
         {
             CHM_THROW_EXCEPTION(chm_error, "Layer creation failed");
         }
 
+
         // Create a polygon from the bounding box
         OGRPolygon polygon;
         OGRLinearRing ring;
@@ -49,11 +59,13 @@ namespace gis
             CHM_THROW_EXCEPTION(chm_error, "Failed to create feature in GeoJSON");
         }
 
+        CSLDestroy(papszOptions);
         OGRFeature::DestroyFeature(poFeature);
+
         GDALClose(poDS);
     }
 
-    void xy2geojson(std::vector<std::tuple<float, float>> xy, std::string filepath, std::string proj4)
+    void xy2geojson(std::vector<std::tuple<float, float>> xy, const std::string& filepath, const std::string& proj4)
     {
         GDALAllRegister();
         const char *pszDriverName = "GeoJSON";
@@ -69,23 +81,21 @@ namespace gis
             CHM_THROW_EXCEPTION(chm_error, "Creation of output file failed");
         }
 
-        OGRLayer *poLayer = poDS->CreateLayer("layer", nullptr, wkbPoint, nullptr);
+        OGRSpatialReference oSRS;
+
+        if (oSRS.importFromProj4(proj4.c_str()) != OGRERR_NONE)
+        {
+            CHM_THROW_EXCEPTION(chm_error, "Failed to import PROJ.4 string");
+        }
+        char **papszOptions = nullptr;
+        papszOptions = CSLAddNameValue(papszOptions, "RFC7946", "TRUE");
+
+        OGRLayer *poLayer = poDS->CreateLayer("layer", &oSRS, wkbPoint, papszOptions);
         if (poLayer == nullptr)
         {
             CHM_THROW_EXCEPTION(chm_error, "Layer creation failed");
         }
 
-//        OGRSpatialReference oSRS;
-
-//        if (oSRS.importFromProj4(proj4.c_str()) != OGRERR_NONE)
-//        {
-//            CHM_THROW_EXCEPTION(chm_error, "Failed to import PROJ.4 string");
-//        }
-//        if (poDS->SetProjection(proj4.c_str()) != CE_None)
-//        {
-//            CHM_THROW_EXCEPTION(chm_error, "Failed to set spatial reference");
-//        }
-
         OGRFieldDefn oFieldX("X", OFTReal);
         oFieldX.SetWidth(32);
         poLayer->CreateField(&oFieldX);
@@ -121,6 +131,7 @@ namespace gis
 
             OGRFeature::DestroyFeature(poFeature);
         }
+        CSLDestroy(papszOptions);
         GDALClose(poDS);
     }
 

src/utility/gis.hpp

@@ -35,7 +35,7 @@
 
 namespace gis
 {
-    void bbox2geojson(double xmin, double ymin, double xmax, double ymax, const std::string& filepath);
-    void xy2geojson(std::vector<std::tuple<float, float>> xy, std::string filepath, std::string proj4);
+    void bbox2geojson(double xmin, double ymin, double xmax, double ymax, const std::string& filepath, const std::string& proj4);
+    void xy2geojson(std::vector<std::tuple<float, float>> xy, const std::string& filepath, const std::string& proj4);
 
 }