Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/MeteoObj.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/MeteoObj.h>
19
20using namespace mio;
21using namespace std;
22
23/************************************************************
24 * static section *
25 ************************************************************/
26const size_t SnGrids::nrOfParameters = SnGrids::lastparam - SnGrids::firstparam + 1;
27std::vector<std::string> SnGrids::paramname;
28const bool SnGrids::__init = SnGrids::initStaticData();
29
30bool SnGrids::initStaticData()
31{
32//the order must be the same as in the enum
33paramname.push_back("TA");
34paramname.push_back("RH");
35paramname.push_back("VW");
36paramname.push_back("DW");
37paramname.push_back("ISWR");
38paramname.push_back("ISWR_DIFF");
39paramname.push_back("ISWR_DIR");
40paramname.push_back("ILWR");
41paramname.push_back("HS");
42paramname.push_back("PSUM");
43paramname.push_back("PSUM_PH");
44paramname.push_back("PSUM_TECH");
45paramname.push_back("GROOMING");
46paramname.push_back("TSG");
47paramname.push_back("TSS");
48paramname.push_back("TS0");
49paramname.push_back("TSNOW");
50paramname.push_back("TSNOW_AVG");
51paramname.push_back("RHOSNOW_AVG");
52paramname.push_back("SWE");
53paramname.push_back("RSNO");
54paramname.push_back("TOP_ALB");
55paramname.push_back("SURF_ALB");
56paramname.push_back("SP");
57paramname.push_back("RB");
58paramname.push_back("RG");
59paramname.push_back("N3");
60paramname.push_back("MS_SNOWPACK_RUNOFF");
61paramname.push_back("MS_SOIL_RUNOFF");
62paramname.push_back("MS_WATER");
63paramname.push_back("SFC_SUBL");
64paramname.push_back("MNS");
65paramname.push_back("STORE");
66paramname.push_back("GLACIER");
67paramname.push_back("GLACIER_EXPOSED");
68paramname.push_back("ET");
69paramname.push_back("TSOIL1");
70paramname.push_back("TSOIL2");
71paramname.push_back("TSOIL3");
72paramname.push_back("TSOIL4");
73paramname.push_back("TSOIL5");
74
75
76if (paramname.size()!=(SnGrids::lastparam+1))
77throw IOException("Wrong number of string representations for the SnGrids parameters! You forgot to update the \"paramname\" vector.", AT);
78
79return true;
80}
81
82const std::string& SnGrids::getParameterName(const size_t& parindex)
83{
84if (parindex >= SnGrids::nrOfParameters)
85throw IndexOutOfBoundsException("Trying to get name for parameter that does not exist", AT);
86
87return paramname[parindex];
88}
89
90size_t SnGrids::getParameterIndex(const std::string& parname)
91{
92for (size_t ii=0; ii<SnGrids::nrOfParameters; ii++) {
93if (paramname[ii] == parname) return ii;
94}
95
96return IOUtils::npos; //parameter not a part of SnGrids
97}
98
99
100/************************************************************
101 * MeteoObj *
102 ************************************************************/
103MeteoObj::MeteoObj(const mio::Config& in_config, const mio::DEMObject& in_dem)
104 : timer(), config(in_config), io(in_config), dem(in_dem),
105 ta(in_dem, IOUtils::nodata), rh(in_dem, IOUtils::nodata), psum(in_dem, IOUtils::nodata),
106 psum_ph(in_dem, IOUtils::nodata), vw(in_dem, IOUtils::nodata), dw(in_dem, IOUtils::nodata), p(in_dem, IOUtils::nodata), ilwr(in_dem, IOUtils::nodata),
107 sum_ta(), sum_rh(), sum_rh_psum(), sum_psum(), sum_psum_ph(), sum_vw(), sum_ilwr(),
108 vecMeteo(), date(), glaciers(NULL), count_sums(0), count_precip(0), skipWind(false) {}
109
110MeteoObj::~MeteoObj()
111{
112if (glaciers!=NULL) delete glaciers;
113}
114
115void MeteoObj::setSkipWind(const bool& i_skipWind) {
116skipWind = i_skipWind;
117}
118
119void MeteoObj::prepare(const mio::Date& in_date)
120{
121if (!MPIControl::instance().master()) // Only master reads data
122return;
123
124date = in_date;
125getMeteo(date);
126}
127
128void MeteoObj::get(const mio::Date& in_date, mio::Grid2DObject& out_ta, mio::Grid2DObject& out_rh, mio::Grid2DObject& out_psum,
129 mio::Grid2DObject& out_psum_ph, mio::Grid2DObject& out_vw, mio::Grid2DObject& out_dw, mio::Grid2DObject& out_p, mio::Grid2DObject& out_ilwr)
130{
131timer.restart(); //this method is called first, so we initiate the timing here
132
133if (MPIControl::instance().master()) {
134if (date.isUndef()) {
135date = in_date;
136getMeteo(date); //it will throw an exception if something goes wrong
137}
138if (in_date != date) {
139cerr << "[w] Meteo data was prepared for " << date.toString(Date::ISO);
140cerr << ", requested for " << in_date.toString(Date::ISO) << ", this is not optimal...\n";
141date = in_date;
142getMeteo(date); //it will throw an exception if something goes wrong
143}
144}
145
146//this acts as a barrier and forces MPI synchronization
147MPIControl::instance().broadcast(ta);
148MPIControl::instance().broadcast(rh);
149MPIControl::instance().broadcast(psum);
150MPIControl::instance().broadcast(psum_ph);
151MPIControl::instance().broadcast(vw);
152MPIControl::instance().broadcast(dw);
153MPIControl::instance().broadcast(p);
154MPIControl::instance().broadcast(ilwr);
155
156out_ta = ta;
157out_rh = rh;
158out_psum = psum;
159out_psum_ph = psum_ph;
160out_vw = vw;
161out_dw = dw;
162out_p = p;
163out_ilwr = ilwr;
164timer.stop();
165}
166
167void MeteoObj::get(const mio::Date& in_date, std::vector<mio::MeteoData>& o_vecMeteo)
168{
169timer.start();
170
171if (MPIControl::instance().master()) {
172if (date.isUndef()) {
173date = in_date;
174getMeteo(date); //it will throw an exception if something goes wrong
175}
176if (in_date != date) {
177cerr << "[w] Meteo data was prepared for " << date.toString(Date::ISO);
178cerr << ", requested for " << in_date.toString(Date::ISO) << ", this is not optimal...\n";
179date = in_date;
180getMeteo(date); //it will throw an exception if something goes wrong
181}
182}
183
184//this acts as a barrier and forces MPI synchronization
185MPIControl::instance().broadcast(vecMeteo);
186
187o_vecMeteo = vecMeteo;
188timer.stop();
189}
190
191double MeteoObj::getTiming() const
192{
193return timer.getElapsed();
194}
195
196void MeteoObj::checkInputsRequirements(std::vector<MeteoData>& vecData)
197{
198//This function checks that the necessary input data are available for the current timestamp
199unsigned int nb_ta=0, nb_iswr=0, nb_rh=0, nb_ilwr=0;
200unsigned int nb_ilwr_ta=0;
201
202if (vecData.empty())
203throw IOException("Vector of input meteo data is empty!", AT);
204
205for (size_t ii=0; ii<vecData.size(); ii++) {
206if (vecData[ii](MeteoData::TA) != IOUtils::nodata) nb_ta++;
207
208if (vecData[ii](MeteoData::ISWR) != IOUtils::nodata) nb_iswr++;
209
210if (vecData[ii](MeteoData::RH) != IOUtils::nodata) nb_rh++;
211
212if (vecData[ii](MeteoData::ILWR) != IOUtils::nodata) {
213nb_ilwr++;
214//We need ILWR and TA at the same location (so that the emissivity can be computed)
215//But since we rely on having one meteo1D station in the code, we also need ISWR at this place
216//(since we use the location of the measurement for computations with ISWR)
217if (vecData[ii](MeteoData::TA) != IOUtils::nodata && vecData[ii](MeteoData::ISWR) != IOUtils::nodata) nb_ilwr_ta++;
218}
219}
220
221if ( nb_ta==0 || nb_iswr==0 || nb_rh==0 || nb_ilwr==0) {
222printf("nb(ta)=%d nb(iswr)=%d nb(rh)=%d nb(ilwr)=%d\n",nb_ta, nb_iswr, nb_rh, nb_ilwr);
223 throw IOException("Not enough input meteo data on "+vecData[0].date.toString(Date::ISO), AT);
224}
225if ( nb_ilwr_ta==0 ) {
226cout << "[e] For emissivity calculation, at least one set of both TA and ILWR are needed at the same station!\n";
227throw IOException("Not enough input meteo data", AT);
228}
229}
230
231void MeteoObj::fillMeteoGrids(const Date& calcDate)
232{
233//fill the meteo parameter grids (of the AlpineControl object) using the data from the stations
234try {
235io.getMeteoData(calcDate, dem, MeteoData::PSUM, psum);
236io.getMeteoData(calcDate, dem, MeteoData::PSUM_PH, psum_ph);
237io.getMeteoData(calcDate, dem, MeteoData::RH, rh);
238io.getMeteoData(calcDate, dem, MeteoData::TA, ta);
239if (!skipWind) {
240io.getMeteoData(calcDate, dem, MeteoData::VW, vw);
241io.getMeteoData(calcDate, dem, MeteoData::DW, dw);
242}
243io.getMeteoData(calcDate, dem, MeteoData::P, p);
244io.getMeteoData(calcDate, dem, MeteoData::ILWR, ilwr);
245cout << "[i] 2D Interpolations done for " << calcDate.toString(Date::ISO) << "\n";
246} catch (long) {
247cout << "[e] at " << calcDate.toString(Date::ISO) << " Could not fill 2D meteo grids" << endl;
248} catch(std::exception& e) {
249cerr << e.what() << endl;
250throw;
251}
252}
253
254void MeteoObj::getMeteo(const Date& calcDate)
255{
256//Note: in case of MPI simulation only master node is responsible for file I/O
257if (!MPIControl::instance().master()) return;
258
259// Collect the Meteo values at each stations
260io.getMeteoData(calcDate, vecMeteo);
261checkInputsRequirements(vecMeteo);
262
263// Now fill the 2D Meteo Fields. Keep in mind that snowdrift might edit these fields
264fillMeteoGrids(calcDate);
265
266cout << "[i] Success reading/preparing meteo data for date: " << calcDate.toString(Date::ISO) << endl;
267}
268
269void MeteoObj::checkLapseRate(const std::vector<mio::MeteoData>& i_vecMeteo, const mio::MeteoData::Parameters& param)
270{
271std::vector<double> vecData, vecAltitudes;
272for (size_t ii=0; ii<i_vecMeteo.size(); ii++){
273const double& val = i_vecMeteo[ii](param);
274if (val != IOUtils::nodata){
275vecData.push_back( val );
276vecAltitudes.push_back( i_vecMeteo[ii].meta.position.getAltitude() );
277}
278}
279
280if (vecData.size()<2) return;
281if (param==MeteoData::PSUM) { //skip when there is no precip
282if (mio::Interpol2D::allZeroes(vecData)) return;
283}
284
285double A, B, R;
286std::string mesg;
287Interpol1D::NoisyLinRegression(vecAltitudes, vecData, A, B, R, mesg);
288const std::string date_str( i_vecMeteo[0].date.toString(Date::ISO) );
289const std::string param_str( MeteoData::getParameterName(param) );
290std::cout << "[check:Data_Lapse_Rate] " << date_str << " " << param_str << " " << std::fixed << std::setw(7) << std::setprecision(5) << A << " ";
291
292if (param==MeteoData::PSUM && A<-1e-3) { //when precip gradient is "too wrong", print values
293const std::streamsize prec = std::cout.precision();
294std::cout << "- ";
295for (size_t ii=0; ii<vecData.size(); ii++){
296std::cout << std::fixed << std::setw(4) << std::setprecision(2) << vecData[ii] << "@" << std::fixed << std::setw(6) << std::setprecision(1) << vecAltitudes[ii] << " ";
297}
298std::cout << std::setprecision(prec);
299std::cout .unsetf(ios_base::floatfield);
300}
301std::cout << "\n";
302}
303
304void MeteoObj::checkGridRange(const mio::Date& calcDate, const mio::Grid2DObject& grid, const mio::MeteoData::Parameters& param)
305{
306if (param==MeteoData::RH) {
307const double min = grid.grid2D.getMin();
308if (min<=0.05)
309std::cout << "[check:Grids_Range_Check] " << calcDate.toString(Date::ISO) << " Rh_min=" << min << "\n";
310} else if (param==MeteoData::VW) {
311const double max = grid.grid2D.getMax();
312if (max>40.)
313std::cout << "[check:Grids_Range_Check] " << calcDate.toString(Date::ISO) << " VW_max=" << max << "\n";
314} else
315throw IOException("Parameter '"+MeteoData::getParameterName(param)+"' not supported here", AT);
316}
317
318void MeteoObj::setGlacierMask(const Grid2DObject& glacierMask)
319{
320if (glacierMask.grid2D.getCount()>0) { //at least one pixel is glaciated...
321glaciers = new Glaciers(config, dem);
322glaciers->setGlacierMap( glacierMask );
323}
324}
325
326//this should only be called when "--nocompute" was set. So we consider that
327//most of the other modules have NOT been called.
328void MeteoObj::checkMeteoForcing(const mio::Date& calcDate)
329{
330if (calcDate != date) {
331cerr << "[w] Meteo data was prepared for " << date.toString(Date::ISO);
332cerr << ", requested for " << calcDate.toString(Date::ISO) << ", this is not optimal...\n";
333date = calcDate;
334getMeteo(date); //it will throw an exception if something goes wrong
335}
336
337if (glaciers!=NULL) {
338glaciers->correctTemperatures( ta );
339}
340
341//produce monthly gridded sums
342int year, month, day, hour, minute, second;
343calcDate.getDate(year, month, day, hour, minute, second);
344const bool startOfMonth = (day==1 && hour==0 && minute==0 && second==0);
345if (startOfMonth || sum_ta.empty()) { //we need to (re-)initialize the sums
346if (startOfMonth && !sum_ta.empty()) {
347sum_ta /= static_cast<double>(count_sums);
348//sum_rh /= static_cast<double>(count_sums);
349if (count_precip>0) sum_rh_psum /= static_cast<double>(count_precip);
350//if (count_precip>0) sum_psum /= static_cast<double>(count_precip);
351if (count_precip>0) sum_psum_ph /= static_cast<double>(count_precip);
352sum_vw /= static_cast<double>(count_sums);
353sum_ilwr /= static_cast<double>(count_sums);
354
355io.write2DGrid(sum_ta, MeteoGrids::TA, calcDate);
356//io.write2DGrid(sum_rh, MeteoGrids::RH, calcDate+1./(3600.*24));
357io.write2DGrid(sum_vw, MeteoGrids::VW, calcDate);
358io.write2DGrid(sum_ilwr, MeteoGrids::ILWR, calcDate);
359if (count_precip>0) io.write2DGrid(sum_rh_psum, MeteoGrids::RH, calcDate);
360if (count_precip>0) io.write2DGrid(sum_psum, MeteoGrids::PSUM, calcDate);
361if (count_precip>0) io.write2DGrid(sum_psum_ph, MeteoGrids::PSUM_PH, calcDate);
362}
363
364count_sums=0;
365count_precip=0;
366sum_ta.set(dem, 0.);
367//sum_rh.set(dem, 0.);
368sum_rh_psum.set(dem, 0.);
369sum_psum.set(dem, 0.);
370sum_psum_ph.set(dem, 0.);
371sum_vw.set(dem, 0.);
372sum_ilwr.set(dem, 0.);
373}
374
375count_sums++;
376sum_ta += ta;
377//sum_rh += rh;
378sum_vw += vw;
379sum_ilwr += ilwr;
380if (psum.grid2D.getMin()>0.) { //there are some precip
381sum_psum += psum;
382sum_psum_ph += psum_ph;
383sum_rh_psum += rh;
384count_precip++;
385}
386
387//check range in the grids
388checkGridRange(calcDate, rh, MeteoData::RH);
389checkGridRange(calcDate, vw, MeteoData::VW);
390
391//now get and print the lapse rates
392vector<MeteoData> vecMd;
393io.getMeteoData(calcDate, vecMd);
394checkLapseRate(vecMd, MeteoData::TA);
395checkLapseRate(vecMd, MeteoData::RH);
396checkLapseRate(vecMd, MeteoData::VW);
397checkLapseRate(vecMd, MeteoData::PSUM);
398}

Archive Download this file

Revision: HEAD