diff --git a/src/core.cpp b/src/core.cpp index b84d2a10..f9389cd5 100644 --- a/src/core.cpp +++ b/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; diff --git a/src/mesh/triangulation.cpp b/src/mesh/triangulation.cpp index e46b2bff..42715f93 100644 --- a/src/mesh/triangulation.cpp +++ b/src/mesh/triangulation.cpp @@ -1451,7 +1451,7 @@ void triangulation::reorder_faces(std::vector 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) diff --git a/src/utility/gis.cpp b/src/utility/gis.cpp index c76dfd74..18725a0e 100644 --- a/src/utility/gis.cpp +++ b/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> xy, std::string filepath, std::string proj4) + void xy2geojson(std::vector> 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); } diff --git a/src/utility/gis.hpp b/src/utility/gis.hpp index 3c39e83d..e7ec4931 100644 --- a/src/utility/gis.hpp +++ b/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> 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> xy, const std::string& filepath, const std::string& proj4); }