From e8218db0ab53f634144f5ff015884e1d9b338606 Mon Sep 17 00:00:00 2001
From: Chris Marsh <chrismarsh.c2@gmail.com>
Date: Tue, 10 Dec 2024 11:15:17 -0600
Subject: [PATCH] Move station to shpfile logic into standalone function

---
 src/core.cpp    |  3 ++-
 src/metdata.cpp | 19 +++++++++++++++----
 src/metdata.hpp |  6 ++++++
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/src/core.cpp b/src/core.cpp
index 60b3478a..93635fea 100644
--- a/src/core.cpp
+++ b/src/core.cpp
@@ -609,8 +609,9 @@ void core::config_forcing(pt::ptree &value)
     }
 
 
-    auto f = output_folder_path / "stations.vtp";
+    auto f = output_folder_path / std::format("stations_{}.vtp", _comm_world.rank());
     _metdata->write_stations_to_ptv(f.string());
+    _metdata->write_stations_to_shp((output_folder_path / std::format("stations_{}.shp",_comm_world.rank())).string());
 
     SPDLOG_DEBUG("Finished reading stations. Took {} s", c.toc<s>());
 
diff --git a/src/metdata.cpp b/src/metdata.cpp
index 9a379dbe..6483963e 100644
--- a/src/metdata.cpp
+++ b/src/metdata.cpp
@@ -241,7 +241,7 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou
 
         SPDLOG_DEBUG("Initializing datastructure");
 
-        std::vector<std::tuple<float, float>> xy;
+
         auto e = _nc->get_z();
 
         // #pragma omp parallel for
@@ -321,12 +321,12 @@ void metdata::load_from_netcdf(const std::string& path, const triangulation::bou
 
                 _dD_tree.insert( boost::make_tuple(Kernel::Point_2(s->x(),s->y()),s) );
 
-                xy.emplace_back(longitude, latitude);
+
             }
         }
         SPDLOG_DEBUG("Done initializing datastructure");
-        gis::xy2shp(xy, "forcing_points.shp", _mesh_proj4);
-        SPDLOG_DEBUG("This rank is using # grid cells = {}", xy.size());
+
+        SPDLOG_DEBUG("This rank is using # grid cells = {}", _stations.size());
         if( skipped == _nstations)
         {
             CHM_THROW_EXCEPTION(forcing_error,
@@ -825,4 +825,15 @@ std::vector< std::shared_ptr<station>>& metdata::stations()
 bool metdata::is_multipart_nc()
 {
     return _is_multipart_nc;
+}
+
+void metdata::write_stations_to_shp(const std::string& fname)
+{
+    std::vector<std::tuple<float, float>> xy;
+    for(auto itr: _stations)
+    {
+        if(itr)
+            xy.emplace_back(itr->x(), itr->y());
+    }
+    gis::xy2shp(xy, fname, _mesh_proj4);
 }
\ No newline at end of file
diff --git a/src/metdata.hpp b/src/metdata.hpp
index c4baeb35..ae83401f 100644
--- a/src/metdata.hpp
+++ b/src/metdata.hpp
@@ -182,6 +182,12 @@ class metdata
      */
     bool is_netcdf();
 
+    /**
+     * Writes the lat and lon of the subsetted forcing poitns to shape file
+     */
+    void write_stations_to_shp(const std::string& fname);
+
+
     /// Subsets all timeseries to begin at [start, end]. For ascii, the underlying timeseries is modified.
     /// For nc, internal offsets are computed to start, end.
     /// This updates the internal start and end times, as well as resets the current time to be = start