Skip to content

NetCDF Forcings Provider: Restore ability to read older-format files #110

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
PhilMiller wants to merge 4 commits into
base: development
Choose a base branch
from
Open

NetCDF Forcings Provider: Restore ability to read older-format files #110

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3ee3c7b
Restore support for older-format NetCDF forcings files that have Time…
PhilMiller Dec 26, 2025
407ed3e
Clean up comments and formatting
PhilMiller Dec 26, 2025
e90c79c
Handle more units, with and without a specified reference epoch
PhilMiller Dec 26, 2025
3002d54
Handle different orderings of variable dimensions
PhilMiller Dec 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 86 additions & 34 deletions src/forcing/NetCDFPerFeatureDataProvider.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,6 @@ NetCDFPerFeatureDataProvider::NetCDFPerFeatureDataProvider(std::string input_pat
// correct string release
nc_free_string(num_ids,&string_buffers[0]);

// Modified code to handle units, epoch start, and reading all time values correctly - KSL

// Get the time variable - getVar collects all values at once and stores in memory
// Extremely large timespans could be problematic, but for ngen use cases, this should not be a problem
auto time_var = nc_file->getVar("Time");
Expand All @@ -147,8 +145,22 @@ NetCDFPerFeatureDataProvider::NetCDFPerFeatureDataProvider(std::string input_pat
std::vector<double> raw_time(num_times);

try {
time_var.getVar(raw_time.data());
} catch(const netCDF::exceptions::NcException& e) {
auto dim_count = time_var.getDimCount();
// Old-format files have dimensions (catchment, time), new-format
// files generated by the forcings engine have just (time)
if (dim_count == 2) {
if (time_var.getDim(0).getName() != "catchment-id" || time_var.getDim(1).getName() != "time") {
Logger::logMsgAndThrowError("In NetCDF file '" + input_path + "', 'Time' variable dimensions don't match expectations");
}
time_var.getVar({0ul, 0ul}, {1ul, num_times}, raw_time.data());
} else if (dim_count == 1) {
time_var.getVar({0ul}, {num_times}, raw_time.data());
} else {
throw std::runtime_error("Unexpected " + std::to_string(dim_count)
+ " dimensions on Time variable in NetCDF file '"
+ input_path + "'");
}
} catch(const std::exception& e) {
netcdf_ss << "Error reading time variable: " << e.what() << std::endl;
LOG(netcdf_ss.str(), LogLevel::WARNING); netcdf_ss.str("");
throw;
Expand All @@ -157,7 +169,6 @@ NetCDFPerFeatureDataProvider::NetCDFPerFeatureDataProvider(std::string input_pat
std::string time_units;
try {
time_var.getAtt("units").getValues(time_units);

} catch(const netCDF::exceptions::NcException& e) {
netcdf_ss << "Error reading time units: " << e.what() << std::endl;
LOG(netcdf_ss.str(), LogLevel::WARNING); netcdf_ss.str("");
Expand All @@ -169,27 +180,43 @@ NetCDFPerFeatureDataProvider::NetCDFPerFeatureDataProvider(std::string input_pat
double time_scale_factor = 1.0;
std::time_t epoch_start_time = 0;

//The following makes some assumptions that NetCDF output from the forcing engine will be relatively uniform
//Specifically, it assumes time values are in units since the Unix Epoch.
//If the forcings engine outputs additional unit formats, this will need to be expanded
if (time_units.find("minutes since") != std::string::npos) {
std::string time_base_unit;
auto since_index = time_units.find("since");
if (since_index != std::string::npos) {
time_base_unit = time_units.substr(0, since_index - 1);

std::string datetime_str = time_units.substr(since_index + 6);
std::tm tm = {};
std::istringstream ss(datetime_str);
ss >> std::get_time(&tm, "%Y-%m-%d %H:%M:%S"); // This may be particularly inflexible
epoch_start_time = timegm(&tm); // timegm may not be available in all environments/OSes ie: Windows
} else {
time_base_unit = time_units;
}

if (time_base_unit == "minutes") {
time_unit = TIME_MINUTES;
time_scale_factor = 60.0;
} else if (time_units.find("hours since") != std::string::npos) {
} else if (time_base_unit == "hours") {
time_unit = TIME_HOURS;
time_scale_factor = 3600.0;
} else {
} else if (time_base_unit == "seconds" || time_base_unit == "s") {
time_unit = TIME_SECONDS;
time_scale_factor = 1.0;
} else if (time_base_unit == "milliseconds" || time_base_unit == "ms") {
time_unit = TIME_MILLISECONDS;
time_scale_factor = 1.0e-3;
} else if (time_base_unit == "microseconds" || time_base_unit == "us") {
time_unit = TIME_MICROSECONDS;
time_scale_factor = 1.0e-6;
} else if (time_base_unit == "nanoseconds" || time_base_unit == "ns") {
time_unit = TIME_NANOSECONDS;
time_scale_factor = 1.0e-9;
} else {
Logger::logMsgAndThrowError("In NetCDF file '" + input_path + "', time unit '" + time_base_unit + "' could not be converted");
}
//This is also based on the NetCDF from the forcings engine, and may not be super flexible
std::string datetime_str = time_units.substr(time_units.find("since") + 6);
std::tm tm = {};
std::istringstream ss(datetime_str);
ss >> std::get_time(&tm, "%Y-%m-%d %H:%M:%S"); //This may be particularly inflexible
epoch_start_time = timegm(&tm); //timegm may not be available in all environments/OSes ie: Windows

time_vals.resize(raw_time.size());
// End modification - KSL

std::transform(raw_time.begin(), raw_time.end(), time_vals.begin(),
[&](const auto& n) {
Expand All @@ -214,13 +241,20 @@ NetCDFPerFeatureDataProvider::NetCDFPerFeatureDataProvider(std::string input_pat
#endif

netcdf_ss << "All time intervals are constant within tolerance." << std::endl;
LOG(netcdf_ss.str(), LogLevel::SEVERE); netcdf_ss.str("");
LOG(netcdf_ss.str(), LogLevel::DEBUG); netcdf_ss.str("");

// determine start_time and stop_time;
start_time = time_vals[0];
stop_time = time_vals.back() + time_stride;

sim_to_data_time_offset = sim_start_date_time_epoch - start_time;

netcdf_ss << "NetCDF Forcing from file '" << input_path << "'"
<< "Start time " << (time_t)start_time
<< ", Stop time " << (time_t)stop_time
<< ", sim_start_date_time_epoch " << sim_start_date_time_epoch
;
LOG(netcdf_ss.str(), LogLevel::DEBUG); netcdf_ss.str("");
}

NetCDFPerFeatureDataProvider::~NetCDFPerFeatureDataProvider() = default;
Expand Down Expand Up @@ -304,7 +338,8 @@ double NetCDFPerFeatureDataProvider::get_value(const CatchmentAggrDataSelector&

auto stride = idx2 - idx1;

std::vector<std::size_t> start, count;
std::vector<std::size_t> start(2), count(2);
std::vector<std::ptrdiff_t> var_index_map(2);

auto cat_pos = id_pos[selector.get_id()];

Expand All @@ -325,16 +360,29 @@ double NetCDFPerFeatureDataProvider::get_value(const CatchmentAggrDataSelector&
//TODO: Currently assuming a whole variable cache slice across all catchments for a single timestep...but some stuff here to support otherwise.
// For reference: https://stackoverflow.com/a/72030286

//Modified to work with NetCDF dimension shapes and fix errors - KSL
size_t cache_slices_t_n = (read_len + cache_slice_t_size - 1) / cache_slice_t_size; // Ceiling division to ensure remainders have a slice


//Explicitly setting dimension shapes
auto dims = ncvar.getDims();
size_t catchment_dim_size = dims[1].getSize();
size_t time_dim_size = dims[0].getSize();

//Cache slicing - modified to work with dimensions structure
int dim_time, dim_catchment;
if (dims.size() != 2) {
Logger::logMsgAndThrowError("Variable dimension count isn't 2");
}
if (dims[0].getName() == "time" && dims[1].getName() == "catchment-id") {
// Forcings Engine NetCDF output case
dim_time = 0;
dim_catchment = 1;
} else if (dims[1].getName() == "time" && dims[0].getName() == "catchment-id") {
// Classic NetCDF file case
dim_time = 1;
dim_catchment = 0;
} else {
Logger::logMsgAndThrowError("Variable dimensions aren't 'time' and 'catchment-id'");
}

size_t time_dim_size = dims[dim_time].getSize();
size_t catchment_dim_size = dims[dim_catchment].getSize();

for( size_t i = 0; i < cache_slices_t_n; i++ ) {
std::shared_ptr<std::vector<double>> cached;
size_t cache_t_idx = idx1 + i * cache_slice_t_size;
Expand All @@ -345,14 +393,18 @@ double NetCDFPerFeatureDataProvider::get_value(const CatchmentAggrDataSelector&
cached = value_cache.get(key).get();
} else {
cached = std::make_shared<std::vector<double>>(catchment_dim_size * slice_size);
start.clear();
start.push_back(cache_t_idx); // start from correct time index
start.push_back(0); // Start from the first catchment
count.clear();
count.push_back(slice_size); // Read the calculated slice size for time
count.push_back(catchment_dim_size); // Read all catchments
start[dim_time] = cache_t_idx; // start from correct time index
start[dim_catchment] = 0; // Start from the first catchment
count[dim_time] = slice_size; // Read the calculated slice size for time
count[dim_catchment] = catchment_dim_size; // Read all catchments
// Whichever order the file stores the data in, the
// resulting array should have all catchments for a given
// time step contiguous
var_index_map[dim_time] = catchment_dim_size;
var_index_map[dim_catchment] = 1;

try {
ncvar.getVar(start,count,&(*cached)[0]);
ncvar.getVar(start,count, {1l, 1l}, var_index_map, cached->data());
value_cache.insert(key, cached);
} catch (netCDF::exceptions::NcException& e) {
netcdf_ss << "NetCDF exception: " << e.what() << std::endl;
Expand All @@ -377,7 +429,7 @@ double NetCDFPerFeatureDataProvider::get_value(const CatchmentAggrDataSelector&
}
}
}
// End modification

rvalue = 0.0;

double a , b = 0.0;
Expand Down
Loading