Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/TechSnowA3D.cc

1/***********************************************************************************/
2/* Copyright 2018-2018 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/TechSnowA3D.h>
19#include <alpine3d/MPIControl.h>
20#include <alpine3d/AlpineMain.h>
21
22#include <errno.h>
23
24#include <algorithm> // for std::find
25#include <iterator>
26
27using namespace std;
28using namespace mio;
29
30/**
31 * @brief Reading the slope file and slope conditions.
32 * @details It reads in the slope file and the defined slope conditions for three different priorities
33 * @param[in] cfg is used to get 'TechSnow' conditions for snow production and reading the slope conditions file, defined in a separate file
34 * @param[in] dem file to check whether the geolocalization of the slope file is the same
35 */
36TechSnowA3D::TechSnowA3D(const mio::Config& cfg, const mio::DEMObject& dem)
37 : skiRunsMap(), psum_tech(dem, IOUtils::nodata), slope_conditions(),
38 start_season(), end_season(), earliest_production(), priorities(),
39 max_snowgun_water(), mean_elevation(dem.grid2D.getMean()), slope_open(), slope_closed()
40{
41const double TZ = cfg.get("TIME_ZONE", "Input");
42
43const std::string mapFile = cfg.get("SKIRUNS_FILE", "TechSnow");
44const bool isMaster = MPIControl::instance().master();
45if (isMaster) {
46mio::IOManager io( cfg );
47io.read2DGrid(skiRunsMap, "/" + mapFile);
48if (!skiRunsMap.isSameGeolocalization( dem ))
49throw InvalidArgumentException("The ski runs map does not has the same geolocalization as the DEM", AT);
50}
51MPIControl::instance().broadcast( skiRunsMap );
52
53// Read slope conditions
54setSlopeConditions(cfg.get("SLOPE_CONDITIONS", "TechSnow"));
55cfg.getValue("SEASON_OPENING", "TechSnow", start_season, TZ);
56cfg.getValue("SEASON_CLOSING", "TechSnow", end_season, TZ);
57cfg.getValue("SLOPE_OPEN", "TechSnow", slope_open);
58cfg.getValue("SLOPE_CLOSED", "TechSnow", slope_closed);
59cfg.getValue("MAX_SNOWGUN_WATER", "TechSnow", max_snowgun_water);
60
61for (unsigned int ii=1; ii<10000; ii++) {
62const std::string prio_key( "PRIO"+mio::IOUtils::toString(ii)+"::START_PROD" );
63if (!cfg.keyExists(prio_key, "TechSnow") ) break; //stop looking for other priorities
64
65priorities.push_back( setSnowStrategy(cfg, TZ, ii) );
66if (earliest_production.isUndef()) earliest_production = priorities.back().startProd;
67if (earliest_production > priorities.back().startProd) earliest_production = priorities.back().startProd;
68}
69}
70
71TechSnowA3D::snowStrategy TechSnowA3D::setSnowStrategy(const mio::Config& cfg, const double& TZ, const unsigned int& nr)
72{
73const std::string root_key( "PRIO"+mio::IOUtils::toString(nr)+"::" );
74TechSnowA3D::snowStrategy ppt;
75
76cfg.getValue(root_key+"start_prod", "TechSnow", ppt.startProd, TZ);
77cfg.getValue(root_key+"end_prod", "TechSnow", ppt.endProd, TZ);
78cfg.getValue(root_key+"start_aim", "TechSnow", ppt.startAim);
79cfg.getValue(root_key+"end_aim", "TechSnow", ppt.endAim);
80cfg.getValue(root_key+"gun_operation", "TechSnow", ppt.gunOperation);
81return ppt;
82}
83
84/**
85 * @brief Creating a vector to store all conditions for the defined slope sections.
86 * @details Creating a vector to store all conditions (slope number, slope area, number of snowguns, priority, minimum snow height, wetbulb treshold)
87 * for each slope section.
88 * @param[in] filename of the file containing the slope conditions for the defined slope sections.
89*/
90void TechSnowA3D::setSlopeConditions(const std::string& filename)
91{
92//first, compute the surface of each slope
93for (size_t ii=0; ii<skiRunsMap.size(); ii++) {
94if (skiRunsMap(ii)==mio::IOUtils::nodata) continue;
95slope_conditions[ getSlopeNumber(skiRunsMap(ii)) ].slope_area++;
96}
97
98std::ifstream fin( filename.c_str() );
99if (fin.fail()) {
100std::ostringstream ss;
101ss << "error opening slope conditions file \"" << filename << "\", possible reason: " << std::strerror(errno);
102throw AccessException(ss.str(), AT);
103}
104
105const char eoln = mio::FileUtils::getEoln(fin); //get the end of line character for the file
106size_t lcount=0;
107try {
108do {
109std::string line;
110getline(fin, line, eoln); //read complete line
111lcount++;
112mio::IOUtils::stripComments(line);
113mio::IOUtils::trim(line);
114if (line.empty()) continue;
115
116std::istringstream iss( line );
117iss.setf(std::ios::fixed);
118iss.precision(std::numeric_limits<double>::digits10);
119
120size_t slope_number;
121unsigned int number_snowguns, priority;
122double min_height, wet_bulb_thresh;
123iss >> std::skipws >> slope_number;
124iss >> std::skipws >> number_snowguns;
125iss >> std::skipws >> priority;
126iss >> std::skipws >> min_height;
127iss >> std::skipws >> wet_bulb_thresh;
128
129if ( !iss || !iss.eof()) {
130std::ostringstream ss;
131ss << "TechSnow: invalid line in file " << filename << " at line " << lcount;
132throw InvalidArgumentException(ss.str(), AT);
133}
134
135slope_conditions[ slope_number ].setUserProperties(number_snowguns, priority, min_height, wet_bulb_thresh);
136} while (!fin.eof());
137
138fin.close();
139} catch (const std::exception&){
140//close fin if open
141if (fin.is_open()) fin.close();
142throw;
143}
144
145const double pixel_area = mio::Optim::pow2( skiRunsMap.cellsize );
146for (std::map<size_t, condition>::iterator it = slope_conditions.begin(); it != slope_conditions.end(); ++it) {
147if (it->second.slope_area==0.) {
148std::cerr << "[W] ski slope " << it->first << " declared in the SLOPE_CONDITIONS file but not found in the SKIRUNS_FILE grid\n";
149slope_conditions.erase( it );
150continue;
151}
152if (it->second.priority==0.)
153throw mio::InvalidArgumentException("Could not find slope " + mio::IOUtils::toString(it->first) + " in the SLOPE_CONDITIONS file", AT);
154
155it->second.slope_area *= pixel_area;
156
157 std::cout << "Found slope " << it->first << " -> " << it->second.toString() << "\n";
158}
159}
160
161/**
162 * @brief Get the slope number for each pixel.
163 * @details Getting the slope number for each pixel. The pixel code is 1xxxx. This function gets rid of the 1 before xxxx.
164 * @param[in] dbl_code of the slope map.
165 * @return value of the slope pixel defined in the SLOPE_CONDITIONS file.
166*/
167size_t TechSnowA3D::getSlopeNumber(const double& dbl_code)
168{
169static const double epsilon = 0.0001;
170return static_cast<size_t>( (dbl_code - 10000. + epsilon) );
171}
172
173/**
174 * @brief Setting priority of the technical snow production
175 * @details The priority of the technical snow production is set of the defined slope conditions in the [TechSnow] in the io.ini file
176 * @param[in] date actual date
177 * @param[in] ppt definition of the snow production strategy
178 * @param[in] snow_height actual snow height
179 * @param[in] slope properties of the slope section
180 * @param[in] date_hour time
181 * @return psum_technical_snow sum of the technical snow production for each pixel.
182*/
183double TechSnowA3D::setPriority(const mio::Date& date,
184 const TechSnowA3D::snowStrategy &ppt, const double& snow_height,
185 const TechSnowA3D::condition& slope,
186 const int date_hour) const
187{
188double psum_technical_snow = IOUtils::nodata;
189
190// Start snow production of winter season; ski resort still closed
191// Technical snow production the whole day
192// Check min snow height is reached before opening
193if (date >= ppt.startProd && date < start_season && date <= ppt.endProd) {
194if (snow_height < slope.min_height * ppt.startAim)
195psum_technical_snow = max_snowgun_water * slope.number_snowguns / (slope.slope_area) * dt_main * (ppt.gunOperation);
196} else if (date >= ppt.startProd && date >= start_season && date <= ppt.endProd) {
197// Continue snow production after season opening; ski resort is open
198// Technical snow production only when ski resort is closed
199if (date_hour > slope_closed || date_hour < slope_open) {
200if (snow_height < slope.min_height * ppt.endAim)
201psum_technical_snow = max_snowgun_water * slope.number_snowguns / (slope.slope_area) * dt_main * (ppt.gunOperation);
202}
203}
204
205// Technical snow production parameters [mm]
206return psum_technical_snow;
207}
208
209/**
210 * @brief Get the grooming and amount of technical snow production map
211 * @details Setting psum_tech and grooming for each pixel based on the defined slope conditions
212 * @param[in] ta temperature map
213 * @param[in] rh relative humidity map
214 * @param[in] hs snow height map
215 * @param[in] date actual date
216*/
217void TechSnowA3D::setMeteo(const mio::Grid2DObject& ta, const mio::Grid2DObject& rh, const mio::Grid2DObject& hs, const mio::Date& date)
218{
219// Cleaning data from last timestep
220psum_tech = IOUtils::nodata;
221
222// Get actual time (hour and minutes)
223int date_hour, date_min;
224date.getTime(date_hour, date_min);
225
226//Period of technical snow production?
227if (date < earliest_production || date > end_season) return;
228
229std::cout << "[I] Production of technical snow has started\n";
230for (size_t ii=0; ii<skiRunsMap.size(); ii++) {
231// Check if slope or not -> if not, do nothing
232if (skiRunsMap(ii)==IOUtils::nodata) continue;
233// Check slope and section
234const size_t slope_nr = getSlopeNumber(skiRunsMap(ii));
235// Wet bulb temperature, in Celsius
236const double T_wet = IOUtils::K_TO_C( mio::Atmosphere::wetBulbTemperature(ta(ii), rh(ii), mean_elevation) );
237
238// Threshold wet bulb temperature
239double psum_technical = IOUtils::nodata;// [mm] Technical snow height in water equivalent
240
241if (T_wet < slope_conditions[slope_nr].wet_bulb_thresh) {
242// Check priority. Since we have cleaned up slope_conditions in setSlopeConditions(), it should always be valid
243const unsigned int priority = slope_conditions[slope_nr].priority - 1; //since vectors start at 0 and not 1!
244psum_technical = setPriority(date, priorities[ priority ], hs(ii), slope_conditions[slope_nr], date_hour);
245}
246
247// Technical snow production
248psum_tech(ii) = psum_technical;
249}
250std::cout << "[I] Production of technical snow has finished\n";
251}
252
253
254mio::Grid2DObject TechSnowA3D::getGrid(const SnGrids::Parameters& param) const
255{
256switch (param) {
257case SnGrids::GROOMING:
258return skiRunsMap; //the caller will decide if it is the right time for grooming or not
259case SnGrids::PSUM_TECH:
260return psum_tech;
261default:
262throw mio::InvalidArgumentException("The requested grid can not be provided by TechSnow", AT);
263}
264}

Archive Download this file

Revision: HEAD