Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/runoff/Runoff.cc

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#include <alpine3d/runoff/Runoff.h>
19
20const double Runoff::DISTANCE_ABSOLUTE_PRECISION = 1e-3; //< [m]
21const double Runoff::MIN_CELL_SIZE = 5.0; //< [m]
22
23/**
24 * @brief Constructor of Runoff instance.
25 * @param in_cfg reference to the config object used by SnowpackInterfaceWorker in order to build our own IOManager in Runoff
26 * @param in_dem reference to the DEM used by Alpine3D
27 * @param in_thresh_rain rain/snow temperature threshold. This will be used to correctly split the output between melt/precipitation
28 */
29Runoff::Runoff(const mio::Config& in_cfg, const mio::DEMObject& in_dem, const double& in_thresh_rain) :
30io(NULL), snowpack(NULL), thresh_rain(in_thresh_rain),
31tz_out(in_cfg.get("TIME_ZONE", "Output")), output_grids(false),
32output_sums(false), catchment_out_path(), resampling_cell_size(),
33grid_size_factor(), n_grid_cells(in_dem.getNx()*in_dem.getNy()),
34catchment_numbering(getCatchmentNumberingScheme(in_cfg)),
35timer(), slope_correction(getSlopeCorrection(in_dem)),
36is_glacier_mask_dynamic(getIsGlacierDynamic(in_cfg)),
37is_glacier_mask_set(false), total_runoff(), glacier_mask(),
38extra_meteo_variables(getExtraMeteoVariables(in_cfg)),
39n_extra_meteo_variables(extra_meteo_variables.size()), catchment_masks(),
40use_external_iomanager_for_grids(false)
41{
42in_cfg.getValue("WRITE_RUNOFF_GRIDS", "OUTPUT", output_grids, mio::IOUtils::nothrow);
43
44std::string catchmentInFile;
45in_cfg.getValue("CATCHMENT", "INPUT", catchmentInFile, mio::IOUtils::nothrow);
46output_sums = !catchmentInFile.empty();
47
48if (output_grids || output_sums) {
49io = getIOManager(in_cfg, output_grids);
50total_runoff.set(in_dem, mio::IOUtils::nodata);
51}
52
53// Set the value of use_external_io_grids if plugin is netcdf and output files
54// are the same
55bool keep_output_files_open = false;
56in_cfg.getValue("NC_KEEP_FILES_OPEN", "Output", keep_output_files_open, mio::IOUtils::nothrow);
57if(keep_output_files_open){
58std::string runoff_grid2d("ARC");
59in_cfg.getValue("RUNOFF_GRID2D", "Output", runoff_grid2d, mio::IOUtils::nothrow);
60mio::IOUtils::toUpper(runoff_grid2d);
61std::string grid2d("ARC");
62in_cfg.getValue("GRID2D", "Output", grid2d, mio::IOUtils::nothrow);
63mio::IOUtils::toUpper(grid2d);
64
65if(runoff_grid2d == "NETCDF" && grid2d == "NETCDF") {
66
67std::string runoff_grid2d_file("NULL");
68in_cfg.getValue("RUNOFF_GRID2DFILE", "Output", runoff_grid2d_file, mio::IOUtils::nothrow);
69mio::IOUtils::toUpper(runoff_grid2d_file);
70
71std::string grid2d_file("NULL");
72in_cfg.getValue("GRID2DFILE", "Output", grid2d_file, mio::IOUtils::nothrow);
73mio::IOUtils::toUpper(grid2d_file);
74
75std::string runoff_grid2d_path("NULL");
76in_cfg.getValue("RUNOFF_GRID2DPATH", "Output", runoff_grid2d_path, mio::IOUtils::nothrow);
77mio::IOUtils::toUpper(runoff_grid2d_path);
78
79std::string grid2d_path("NULL");
80in_cfg.getValue("GRID2DPATH", "Output", grid2d_path, mio::IOUtils::nothrow);
81mio::IOUtils::toUpper(grid2d_path);
82
83if(runoff_grid2d_file == grid2d_file && runoff_grid2d_path == grid2d_path) {
84use_external_iomanager_for_grids=true;
85}
86}
87}
88
89if (output_sums) {
90in_cfg.getValue("CATCHMENTS_PATH", "OUTPUT", catchment_out_path);
91mio::Grid2DObject catchmentGrid;
92io->read2DGrid(catchmentGrid, catchmentInFile);
93resampling_cell_size = getResamplingCellSize(in_dem, catchmentGrid);
94grid_size_factor = in_dem.cellsize/resampling_cell_size;
95constructCatchmentMasks(catchmentGrid);
96initializeOutputFiles(in_dem);
97glacier_mask.set(in_dem, mio::IOUtils::nodata);
98}
99if (MPIControl::instance().master()) {
100std::cout << "[i] Runoff initialised";
101if (output_sums) std::cout << " - Number of catchments: " << catchment_masks.size() << "\n";
102else std::cout << "\n";
103}
104}
105
106/*
107 * @brief Copy constructor
108 * @param copy Runoff object to copy
109 */
110Runoff::Runoff(const Runoff& copy) :
111io(new mio::IOManager(*(copy.io))), snowpack(copy.snowpack),
112thresh_rain(copy.thresh_rain), tz_out(copy.tz_out),
113output_grids(copy.output_grids), output_sums(copy.output_sums),
114catchment_out_path(copy.catchment_out_path),
115resampling_cell_size(copy.resampling_cell_size),
116grid_size_factor(copy.grid_size_factor), n_grid_cells(copy.n_grid_cells),
117catchment_numbering(copy.catchment_numbering),
118timer(copy.timer), slope_correction(copy.slope_correction),
119is_glacier_mask_dynamic(copy.is_glacier_mask_dynamic),
120is_glacier_mask_set(copy.is_glacier_mask_set),
121total_runoff(copy.total_runoff), glacier_mask(copy.glacier_mask),
122extra_meteo_variables(copy.extra_meteo_variables),
123n_extra_meteo_variables(copy.n_extra_meteo_variables),
124catchment_masks(copy.catchment_masks)
125{}
126
127std::string Runoff::getGridsRequirements() const
128{
129return "PSUM TA MS_SOIL_RUNOFF GLACIER MS_SNOWPACK_RUNOFF";
130}
131
132
133/**
134 * @brief Sets the internal reference to SnowpackInterface object
135 * @param sn_interface Reference to the SnowpackInterface object
136 */
137void Runoff::setSnowPack(SnowpackInterface &sn_interface) {
138snowpack = &sn_interface;
139}
140
141
142/**
143 * @brief Writes the results for a specific day
144 * @param[in] i_date the date at which the results should be written
145 * @param[in] psum grid of precipitation in mm/h
146 * @param[in] ta grid of air temperature
147 * @param io_in External IOManager that can be used to write all outputs into the same file (if GRID2DFILE=RUNOFF_GRID2DFILE)
148 */
149void Runoff::output(const mio::Date& i_date, const mio::Grid2DObject& psum,
150const mio::Grid2DObject& ta, mio::IOManager& io_in)
151{
152if (!output_grids && !output_sums) return;
153
154timer.restart();
155updateTotalRunoffGrid();
156
157if(MPIControl::instance().master() && output_grids){
158if(use_external_iomanager_for_grids){
159io_in.write2DGrid(total_runoff, mio::MeteoGrids::ROT, i_date);
160}
161else{
162io->write2DGrid(total_runoff, mio::MeteoGrids::ROT, i_date);
163}
164}
165
166if (output_sums) {
167if(!is_glacier_mask_set || is_glacier_mask_dynamic) {
168updateGlacierMask();
169is_glacier_mask_set = true;
170}
171
172//Compute precip, glacier and melt runoffs
173mio::Grid2DObject precipRunoff(computePrecipRunoff(psum, ta));
174mio::Grid2DObject meltRunoff(total_runoff - precipRunoff);
175mio::Grid2DObject glacierRunoff(meltRunoff*glacier_mask);
176mio::Grid2DObject snowRunoff(meltRunoff - glacierRunoff);
177mio::Grid2DObject totalRunoff(total_runoff);
178
179//Get the grids of the additional meteo variables
180std::vector<mio::Grid2DObject> extraGrids;
181getExtraMeteoGrids(extraGrids);
182
183
184if(MPIControl::instance().master()) {
185//Resample runoff grids to match the cell size of the catchment masks
186if (fabs(grid_size_factor - 1.0) > 1e-5) {
187totalRunoff = mio::LibResampling2D::Nearest(totalRunoff, grid_size_factor);
188precipRunoff = mio::LibResampling2D::Nearest(precipRunoff, grid_size_factor);
189snowRunoff = mio::LibResampling2D::Nearest(snowRunoff, grid_size_factor);
190glacierRunoff = mio::LibResampling2D::Nearest(glacierRunoff, grid_size_factor);
191for (size_t iVar(0); iVar < n_extra_meteo_variables; ++iVar)
192extraGrids[iVar] = mio::LibResampling2D::Nearest(extraGrids[iVar], grid_size_factor);
193}
194
195//Sum the runoffs over each mask and write them in the output files
196double currTotalRunoff, currPrecipRunoff, currSnowRunoff, currGlacierRunoff;
197std::vector<double> currMeteoVars(n_extra_meteo_variables);
198std::map<size_t, mio::Grid2DObject>::const_iterator itMask;
199for (itMask = catchment_masks.begin(); itMask != catchment_masks.end(); ++itMask) {
200currTotalRunoff = sumOverMask(totalRunoff, itMask->second);
201currPrecipRunoff = sumOverMask(precipRunoff, itMask->second);
202currSnowRunoff = sumOverMask(snowRunoff, itMask->second);
203currGlacierRunoff = sumOverMask(glacierRunoff, itMask->second);
204for (size_t iVar(0); iVar < n_extra_meteo_variables; ++iVar)
205currMeteoVars[iVar] = averageOverMask(extraGrids[iVar], itMask->second);
206
207updateOutputFile(itMask->first, i_date, currTotalRunoff,
208currPrecipRunoff, currSnowRunoff, currGlacierRunoff,
209currMeteoVars);
210}
211}
212}
213timer.stop();
214}
215
216
217/**
218 * @brief Destructor of class Runoff
219 */
220Runoff::~Runoff()
221{
222delete io;
223}
224
225
226/**
227 * @brief Initializes private attribute catchment_masks
228 * @param catchmentGrid grid defining the catchments. The catchment numbering
229 * scheme must be specified in the ini file using the key CATCHMENT_NUMBERING
230 * in section INPUT. This scheme can be either ALPINE3D_OLD (for catchments
231 * numbered with powers of 2), or TAUDEM (for standard numbering).
232 */
233void Runoff::constructCatchmentMasks(mio::Grid2DObject catchmentGrid)
234{
235const double factor = catchmentGrid.cellsize/resampling_cell_size;
236if (fabs(factor - 1.0) > 1e-5) {
237catchmentGrid = mio::LibResampling2D::Nearest(catchmentGrid, factor);
238}
239mio::Grid2DObject emptyGrid(catchmentGrid, mio::IOUtils::nodata);
240
241std::vector<size_t> currIndices;
242for (size_t iCell = 0; iCell < catchmentGrid.size(); ++iCell) {
243if (catchmentGrid(iCell) == mio::IOUtils::nodata) continue;
244const longuint currValue = static_cast<longuint>( round(catchmentGrid(iCell)) );
245
246if (catchment_numbering == Alpine3DOld) {
247currIndices = factorizeCatchmentNumber(currValue);
248} else {
249if (currValue > std::numeric_limits<size_t>::max()) {
250std::ostringstream os;
251os << "The ID numbers of some of the subwatersheds defined in "
252 << "the catchment file exceed the maximum index value ("
253 << std::fixed << std::numeric_limits<size_t>::max() << ")."
254 << " Did you forget to add \"CATCHMENT_NUMBERING = ALPINE3D_OLD\""
255 << " in section [INPUT] of your configuration file?";
256throw mio::IndexOutOfBoundsException(os.str(), AT);
257}
258currIndices.assign(1, static_cast<size_t>( currValue ));
259}
260
261for (std::vector<size_t>::const_iterator it = currIndices.begin(); it != currIndices.end(); ++it) {
262if (catchment_masks.find(*it) == catchment_masks.end())
263catchment_masks[*it] = emptyGrid;
264catchment_masks[*it](iCell) = 1.0;
265}
266}
267
268for (std::map<size_t, mio::Grid2DObject>::iterator it = catchment_masks.begin(); it != catchment_masks.end(); ++it) {
269cropMask(it->second);
270}
271}
272
273
274/**
275 * @brief Updates protected attribute total_runoff.
276 */
277void Runoff::updateTotalRunoffGrid()
278{
279total_runoff = snowpack->getGrid(SnGrids::MS_SOIL_RUNOFF);
280
281if ( (total_runoff.getNx() != slope_correction.getNx()) ||
282 (total_runoff.getNy() != slope_correction.getNy()) )
283throw mio::InvalidArgumentException("DEM and soil runoff grid have "
284"incompatible dimensions!", AT);
285
286total_runoff *= slope_correction;
287}
288
289
290/**
291 * @brief Updates protected attribute glacier_mask.
292 */
293void Runoff::updateGlacierMask() {
294glacier_mask = snowpack->getGrid(SnGrids::GLACIER);
295
296if ( (glacier_mask.getNx() != slope_correction.getNx()) ||
297 (glacier_mask.getNy() != slope_correction.getNy()) )
298throw mio::InvalidArgumentException("DEM and glacier mask grid have "
299"incompatible dimensions!", AT);
300
301for (size_t iCell = 0; iCell < n_grid_cells; ++iCell)
302glacier_mask(iCell) = static_cast<double>(glacier_mask(iCell) == mio::IOUtils::nodata);
303}
304
305
306/**
307 * @brief Computes the grid storing the runoff which originates from liquid
308 * precipitation
309 * @param psum grid of precipitation in mm/h
310 * @param ta grid of air temperature
311 * @return Precipitation runoff grid
312 */
313mio::Grid2DObject Runoff::computePrecipRunoff(const mio::Grid2DObject& psum, const mio::Grid2DObject& ta) const
314{
315const mio::Grid2DObject surfRunoff( snowpack->getGrid(SnGrids::MS_SNOWPACK_RUNOFF) );
316if ( (psum.getNx() != slope_correction.getNx()) ||
317 (ta.getNx() != slope_correction.getNx()) ||
318 (surfRunoff.getNx() != slope_correction.getNx()) ||
319 (psum.getNy() != slope_correction.getNy()) ||
320 (ta.getNy() != slope_correction.getNy()) ||
321 (surfRunoff.getNy() != slope_correction.getNy()) )
322throw mio::InvalidArgumentException("DEM and input grids have incompatible dimensions!", AT);
323
324mio::Grid2DObject precipRunoff(total_runoff, mio::IOUtils::nodata);
325for (size_t iCell = 0; iCell < n_grid_cells; ++iCell) {
326const double& precip = psum(iCell);
327const double& runoff = total_runoff(iCell);
328if (runoff == mio::IOUtils::nodata || precip == mio::IOUtils::nodata ||
329ta(iCell) == mio::IOUtils::nodata || surfRunoff(iCell) == mio::IOUtils::nodata)
330precipRunoff(iCell) = mio::IOUtils::nodata;
331else if (ta(iCell) > thresh_rain)
332precipRunoff(iCell) = std::min(precip, runoff);
333else
334precipRunoff(iCell) = 0.0;
335}
336
337return precipRunoff;
338}
339
340
341/**
342 * @brief Returns the grids corresponding to the additional meteo variables
343 * which have to be averaged over the subcatchments and written in the output
344 * files.
345 * @param[out] grids vector containing the meteo grids
346 */
347void Runoff::getExtraMeteoGrids(std::vector<mio::Grid2DObject>& grids) const {
348grids.clear();
349grids.reserve(n_extra_meteo_variables);
350
351for(size_t iVar(0); iVar < n_extra_meteo_variables; ++iVar) {
352const SnGrids::Parameters& currParam(extra_meteo_variables.at(iVar));
353const mio::Grid2DObject tmp( snowpack->getGrid(currParam) );
354
355if(tmp.grid2D.getCount() == 0)
356throw mio::InvalidArgumentException("Cannot average parameter " +
357SnGrids::getParameterName(currParam) + " over the "
358"subcatchments: cannot retrieve the parameter value over "
359"the entire watershed", AT);
360else if ( (tmp.getNx() != slope_correction.getNx()) ||
361 (tmp.getNy() != slope_correction.getNy()) )
362throw mio::InvalidArgumentException("DEM and " +
363SnGrids::getParameterName(currParam) + " grid have "
364"incompatible dimensions!", AT);
365
366grids.push_back(tmp);
367}
368}
369
370
371/**
372 * @brief Initializes the SMET files in which the per-catchment-aggregated
373 * runoff values will be written
374 * @param dem Reference to the DEM object
375 */
376void Runoff::initializeOutputFiles(const mio::Grid2DObject& dem) const
377{
378std::stringstream ss;
379const double cellArea = pow(catchment_masks.begin()->second.cellsize, 2); // in m^2
380double catchArea;
381
382for (std::map<size_t, mio::Grid2DObject>::const_iterator itMask = catchment_masks.begin();
383itMask != catchment_masks.end(); ++itMask, ss.str(""), ss.clear()) {
384catchArea = double(itMask->second.grid2D.getCount())*cellArea*1e-6; //< in km^2
385
386ss << "catch" << std::setfill('0') << std::setw(2) << itMask->first;
387const std::string id = ss.str();
388
389const std::string filename = catchment_out_path + "/" + id + ".smet";
390std::ofstream smet_out; //Output file streams
391smet_out.open(filename.c_str());
392if (smet_out.fail()) throw mio::AccessException(filename.c_str(), AT);
393
394smet_out << "SMET 1.1 ASCII\n";
395smet_out << "[HEADER]\n";
396smet_out << "station_id = " << id << "\n";
397smet_out << "station_name = " << id << "\n";
398smet_out << std::right;
399smet_out << std::fixed;
400
401smet_out << "latitude = " << std::setw(11) << std::setprecision(8) << dem.llcorner.getLat() << "\n";
402smet_out << "longitude = " << std::setw(11) << std::setprecision(8) << dem.llcorner.getLon() << "\n";
403smet_out << "altitude = " << std::setw(7) << std::setprecision(1) << dem.llcorner.getAltitude() << "\n";
404smet_out << "catchment_surface = " << std::setw(7) << std::setprecision(3) << catchArea << "\n";
405smet_out << "comment = surface is in km^2, runoff (Ro) in m^3/h\n";
406smet_out << "nodata = " << std::setw(7) << std::setprecision(0) << mio::IOUtils::nodata << "\n";
407smet_out << "tz = " << std::setw(7) << std::setprecision(0) << tz_out << "\n";
408
409smet_out << "fields = timestamp precip_Ro snowmelt_Ro glacier_melt_Ro total_Ro";
410for (size_t iVar(0); iVar < n_extra_meteo_variables; ++iVar)
411smet_out << " " << SnGrids::getParameterName(extra_meteo_variables.at(iVar));
412smet_out << "\n[DATA]\n";
413
414smet_out.close();
415}
416}
417
418
419/**
420 * @brief Writes the runoff values aggregated over a given catchment in the
421 * corresponding SMET file
422 * @param catchId catchment id number
423 * @param currTime time corresponding to the runoff values
424 * @param totalRunoff total runoff (corrected for slope) over the catchment,
425 * in mm/h
426 * @param precipRunoff part of the total runoff attributable to liquid
427 * precipitation, in mm/h
428 * @param snowRunoff part of the total runoff attributable to snow melt, in
429 * mm/h
430 * @param glacierRunoff part of the total runoff attributable to glacier melt,
431 * in mm/h
432 * @param meteoVars additional meteorological variables averaged over the
433 * catchment area. The units of each variable correspond to those used in
434 * Alpine3D or Snowpack.
435 */
436void Runoff::updateOutputFile(const size_t& catchId, const mio::Date& currTime,
437const double& totalRunoff, const double& precipRunoff,
438const double& snowRunoff, const double& glacierRunoff,
439const std::vector<double>& meteoVars) const
440{
441const double cellArea = resampling_cell_size*resampling_cell_size;
442
443std::stringstream ss;
444ss << "catch" << std::setfill('0') << std::setw(2) << catchId;
445const std::string filename = catchment_out_path + "/" + ss.str() + ".smet";
446std::ofstream smet_out; //Output file streams
447smet_out.open(filename.c_str(), std::ofstream::out | std::ofstream::app);
448if (smet_out.fail()) throw mio::AccessException(filename.c_str(), AT);
449
450smet_out.fill(' ');
451smet_out << std::right;
452smet_out << std::fixed;
453smet_out << currTime.toString(mio::Date::ISO) << " ";
454smet_out << std::setw(10) << std::setprecision(3) << cellArea*precipRunoff*1e-3 << " ";
455smet_out << std::setw(10) << std::setprecision(3) << cellArea*snowRunoff*1e-3 << " ";
456smet_out << std::setw(10) << std::setprecision(3) << cellArea*glacierRunoff*1e-3 << " ";
457smet_out << std::setw(10) << std::setprecision(3) << cellArea*totalRunoff*1e-3;
458for (size_t iVar(0); iVar < n_extra_meteo_variables; ++iVar)
459smet_out << " " << std::setw(10) << std::setprecision(2) << meteoVars.at(iVar);
460
461smet_out << "\n";
462smet_out.close();
463}
464
465
466/**
467 * @brief Returns a IOManager object which will write the runoff grids in the
468 * correct output folder
469 * @param cfg Config object used in Alpine3D
470 * @param outputGrids should the runoff grids be written out?
471 */
472mio::IOManager* Runoff::getIOManager(mio::Config cfg, const bool& outputGrids)
473{
474if (outputGrids) {
475std::string runoff_grid2d("ARC");
476cfg.getValue("RUNOFF_GRID2D", "Output", runoff_grid2d, mio::IOUtils::nothrow);
477mio::IOUtils::toUpper(runoff_grid2d);
478
479cfg.addKey("GRID2D", "Output", runoff_grid2d);
480
481if(runoff_grid2d == "NETCDF") {
482const std::string runoff_grid2d_file = cfg.get("RUNOFF_GRID2DFILE", "Output");
483cfg.addKey("GRID2DFILE", "Output", runoff_grid2d_file);
484} else {
485const std::string runoff_grid2d_path = cfg.get("RUNOFF_GRID2DPATH", "Output");
486cfg.addKey("GRID2DPATH", "Output", runoff_grid2d_path);
487}
488}
489
490return new mio::IOManager(cfg);
491}
492
493
494/**
495 * @brief Computes the grid used to correct the runoff values for slope
496 * (Alpine3D assumes the cells to be flat, the predicted runoffs therefore
497 * have to be corrected for terrain inclination)
498 * @param dem Reference to the DEM object used by Alpine3D
499 */
500mio::Grid2DObject Runoff::getSlopeCorrection(const mio::DEMObject& dem)
501{
502const double to_rad = M_PI/180.;
503
504mio::Grid2DObject slope_corr(dem, mio::IOUtils::nodata);
505for (size_t iCell = 0; iCell < dem.getNx()*dem.getNy(); ++iCell) {
506if (dem(iCell) != mio::IOUtils::nodata)
507slope_corr(iCell) = 1.0/cos(dem.slope(iCell)*to_rad);
508}
509
510return slope_corr;
511}
512
513
514/**
515 * @brief Returns the numbering scheme used in the grid defining the catchments
516 * @param in_cfg Reference to the Config object used in Alpine3D
517 * @return Catchment numbering scheme
518 */
519Runoff::NumberingType Runoff::getCatchmentNumberingScheme(const mio::Config& in_cfg)
520{
521NumberingType type;
522std::string numbering("Alpine3D_old");
523in_cfg.getValue("CATCHMENT_NUMBERING", "INPUT", numbering, mio::IOUtils::nothrow);
524mio::IOUtils::toUpper(numbering);
525if (numbering == "TAUDEM")
526type = TauDEM;
527else if (numbering == "ALPINE3D_OLD")
528type = Alpine3DOld;
529else
530throw mio::IOException("Key CATCHMENT_NUMBERING can only be either "
531"'TauDEM' or 'Alpine3D_old'", AT);
532return type;
533}
534
535
536/**
537 * @brief Returns whether the glaciers are dynamic, i.e. whether the glacier
538 * mask should be expected to change from one time step to the other.
539 * @param cfg Reference to the Config object holding the simulation parameters
540 * @return Is glacier mask dynamic?
541 */
542bool Runoff::getIsGlacierDynamic(const mio::Config& cfg) {
543bool isDynamic(false);
544cfg.getValue("MASK_DYNAMIC", "Output", isDynamic, mio::IOUtils::nothrow);
545return isDynamic;
546}
547
548
549/**
550 * @brief This method returns the meteo variables which have to be averaged over
551 * the subwatershed areas and written in the output files on top of runoff.
552 * @param cfg Reference to the Config object holding the simulation parameters
553 * @return Variables which have to be averaged over the subwatershed areas
554 */
555std::vector<SnGrids::Parameters> Runoff::getExtraMeteoVariables(const mio::Config& cfg) {
556std::vector<std::string> extraDataNames;
557cfg.getValue("RUNOFF_FILES_EXTRA_DATA", "Output", extraDataNames, mio::IOUtils::nothrow);
558const size_t nData(extraDataNames.size());
559
560std::vector<SnGrids::Parameters> extraData;
561extraData.reserve(nData);
562for(size_t iData(0); iData < nData; ++iData) {
563const size_t iParam(SnGrids::getParameterIndex(extraDataNames[iData]));
564extraData.push_back((SnGrids::Parameters)iParam);
565}
566
567return extraData;
568}
569
570
571/**
572 * @brief Returns the cell size which must be used to resample both the
573 * catchment masks and the runoff grids
574 * @param in_dem reference to the Config object used in Alpine3D
575 * @param catchmentGrid grid defining the catchments
576 * @return Resampling cell size in m
577 */
578double Runoff::getResamplingCellSize(const mio::DEMObject& in_dem,
579const mio::Grid2DObject& catchmentGrid)
580{
581double resampCellSize;
582const double minSize = std::min(catchmentGrid.cellsize, in_dem.cellsize);
583const double maxSize = std::max(catchmentGrid.cellsize, in_dem.cellsize);
584if (isMultiple(maxSize, minSize)) {
585const double llxOffset = fabs(catchmentGrid.llcorner.getEasting() -
586in_dem.llcorner.getEasting());
587const double llyOffset = fabs(catchmentGrid.llcorner.getNorthing() -
588in_dem.llcorner.getNorthing());
589const double xCellSize(estimateResamplingCellSize(llxOffset, minSize));
590const double yCellSize(estimateResamplingCellSize(llyOffset, minSize));
591if (isMultiple(xCellSize, yCellSize) || isMultiple(yCellSize, xCellSize))
592resampCellSize = std::min(xCellSize, yCellSize);
593else
594resampCellSize = MIN_CELL_SIZE;
595} else {
596resampCellSize = MIN_CELL_SIZE;
597}
598if (resampCellSize <= MIN_CELL_SIZE + DISTANCE_ABSOLUTE_PRECISION)
599std::cout << "[w] The resolution of the grid defining the sub-catchments "
600 << "matches poorly with the resolution of the DEM. The runoff "
601 << "module will re-interpolate all runoff grids using a new cell "
602 << "size of " << MIN_CELL_SIZE << " meters, which might "
603 << "significantly impact simulation time in case of a large DEM. "
604 << "It is recommended that you define your sub-catchments using the "
605 << "same DEM as the one used by Alpine3D to solve this issue.\n";
606
607return resampCellSize;
608}
609
610
611/**
612 * @brief Returns whether the first input is a multiple of the second one
613 */
614bool Runoff::isMultiple(const double& a, const double& b) {
615const double res = fmod(a, b);
616return fabs(res - b/2.0) > b/2.0 - DISTANCE_ABSOLUTE_PRECISION;
617}
618
619
620double Runoff::estimateResamplingCellSize(const double& llOffset,
621const double& currSizeEstimate) {
622double cellSize;
623if (isMultiple(llOffset, currSizeEstimate)) {
624cellSize = std::max(currSizeEstimate, MIN_CELL_SIZE);
625} else {
626double xShift = fmod(llOffset, currSizeEstimate);
627xShift = std::min(xShift, currSizeEstimate - xShift);
628cellSize = currSizeEstimate/round(currSizeEstimate/xShift);
629cellSize = std::max(cellSize, MIN_CELL_SIZE);
630}
631return cellSize;
632}
633
634
635/**
636 * @brief This function splits up a given unsigned int value into a sum of powers of 2
637 * @param value The value that shall be split up into a sum of powers of two
638 * @return A vector that will hold all exponents of the sum of powers
639 */
640std::vector<size_t> Runoff::factorizeCatchmentNumber(longuint value)
641{
642std::vector<size_t> result;
643while (value != 0) {
644const double lnresult = log(double(value))/log(2);
645const size_t exp = (unsigned int)floor(lnresult);
646const longuint pot = longuint(floor(pow(double(2), double(exp)) + 0.00001)); //numerically stable
647result.push_back(exp);
648value -= pot;
649}
650return result;
651}
652
653
654/**
655 * @brief Crops the mask so as to remove as many nodata cells as possible
656 * (note: this method possibly changes the georeferencing of the mask by
657 * changing its lower-left corner as well as its number of rows or columns)
658 */
659void Runoff::cropMask(mio::Grid2DObject& mask)
660{
661size_t min_ix(mask.getNx()-1), max_ix(0), min_iy(mask.getNy()-1), max_iy(0);
662for (size_t iy = 0; iy < mask.getNy(); ++iy) {
663for (size_t ix = 0; ix < mask.getNx(); ++ix) {
664if (mask(ix,iy) != mio::IOUtils::nodata) {
665if (ix < min_ix) min_ix = ix;
666if (ix > max_ix) max_ix = ix;
667if (iy < min_iy) min_iy = iy;
668if (iy > max_iy) max_iy = iy;
669}
670}
671}
672mio::Grid2DObject tmp(mask, min_ix, min_iy, max_ix-min_ix+1, max_iy-min_iy+1);
673mask = tmp;
674}
675
676
677/**
678 * @brief Sums the values of the grid cells which are located over the mask
679 * @param grid grid whose cell values have to be summed
680 * @param mask mask over which the grid cell values have to be summed
681 * @return Sum of the masked grid cell values
682 */
683double Runoff::sumOverMask(const mio::Grid2DObject& grid, const mio::Grid2DObject& mask)
684{
685mio::Coords llCorner(mask.llcorner);
686grid.gridify(llCorner);
687mio::Grid2DObject subgrid;
688try {
689mio::Grid2DObject tmp(grid, llCorner.getGridI(), llCorner.getGridJ(),
690 mask.getNx(), mask.getNy());
691subgrid = tmp;
692} catch (...) {
693throw mio::InvalidFormatException("Catchment mask extends beyond the DEM boundaries", AT);
694}
695
696double sum(0.0);
697for (size_t iCell = 0; iCell < mask.getNx()*mask.getNy(); ++iCell) {
698if (mask(iCell) != mio::IOUtils::nodata && subgrid(iCell) != mio::IOUtils::nodata)
699sum += subgrid(iCell);
700}
701return sum;
702}
703
704
705/**
706 * @brief Averages the values of the grid cells which are located over the mask
707 * @param grid grid whose cell values have to be averaged
708 * @param mask mask over which the grid cell values have to be averaged
709 * @return Average of the masked grid cell values
710 */
711double Runoff::averageOverMask(const mio::Grid2DObject& grid, const mio::Grid2DObject& mask)
712{
713return sumOverMask(grid, mask)/double(mask.grid2D.getCount());
714}
715
716double Runoff::getTiming() const
717{
718return timer.getElapsed();
719}

Archive Download this file

Revision: HEAD