Skip to content

Commit

Permalink
Move station to shpfile logic into standalone function
Browse files Browse the repository at this point in the history
  • Loading branch information
Chrismarsh committed Dec 10, 2024
1 parent a28ef50 commit e8218db
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 5 deletions.
3 changes: 2 additions & 1 deletion src/core.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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>());

Expand Down
19 changes: 15 additions & 4 deletions src/metdata.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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);
}
6 changes: 6 additions & 0 deletions src/metdata.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit e8218db

Please sign in to comment.