WSL/SLF GitLab Repository

ProcShade.cc 8.86 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/***********************************************************************************/
/*  Copyright 2012 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <fstream>
#include <sstream>
#include <errno.h>
#include <cstring>
#include <algorithm>

#include <meteoio/meteoFilters/ProcShade.h>
25
#include <meteoio/meteoLaws/Sun.h>
26
#include <meteoio/IOHandler.h>
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <meteoio/FileUtils.h>

using namespace std;

namespace mio {
//custom function for sorting cache_meteo_files
struct sort_pred {
	bool operator()(const std::pair<double,double> &left, const std::pair<double,double> &right) {
		if (left.first < right.first) return true;
		return false;
	}
};

40
const double ProcShade::diffuse_thresh = 15.; //below this threshold, not correction is performed since it will only be diffuse
41

42
ProcShade::ProcShade(const std::vector<std::string>& vec_args, const std::string& name, const Config &i_cfg)
43
        : ProcessingBlock(name), cfg(i_cfg), dem(), masks()
44
45
46
47
48
49
50
51
52
{
	parse_args(vec_args);
	properties.stage = ProcessingProperties::first; //for the rest: default values
}

void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
                        std::vector<MeteoData>& ovec)
{
	ovec = ivec;
53
	if (ovec.empty()) return;
54
	
55
56
	const std::string stationHash( ovec[0].meta.getHash() );
	SunObject Sun;
57
	
58
59
60
61
62
63
64
	//check if the station already has an associated mask, first as wildcard then by station hash
	std::map< std::string , std::vector< std::pair<double,double> > >::iterator mask = masks.find( "*" );
	if (mask==masks.end()) {
		//now look for our specific station hash
		mask = masks.find( stationHash );
		if (mask==masks.end()) {
			std::vector< std::pair<double,double> > tmp_mask;
65
			computeMask(dem, ovec[0].meta, tmp_mask);
66
67
68
			masks[ stationHash ] = tmp_mask;
			mask = masks.find( stationHash);
		}
69
70
	}
	
71
72
73
	double Md_prev = IOUtils::nodata;
	double julian_prev = 0.;
	for (size_t ii=0; ii<ovec.size(); ii++) { //now correct all timesteps
74
75
		double& tmp = ovec[ii](param);
		if (tmp == IOUtils::nodata) continue; //preserve nodata values
76
		if (tmp<diffuse_thresh) continue; //only diffuse radiation, there is nothing to correct
77
78
79
80

		const Coords position( ovec[ii].meta.position );
		Sun.setLatLon(position.getLat(), position.getLon(), position.getAltitude()); //if they are constant, nothing will be recomputed
		Sun.setDate(ovec[ii].date.getJulian(true), 0.); //quicker: we stick to gmt
81
		double sun_azi, sun_elev;
82
		Sun.position.getHorizontalCoordinates(sun_azi, sun_elev);
83
		
84
		const double mask_elev = getMaskElevation(mask->second, sun_azi);
85
		if (mask_elev>0 && mask_elev>sun_elev) { //the point is in the shade
86
87
88
89
90
91
92
93
			const double TA=ovec[ii](MeteoData::TA), RH=ovec[ii](MeteoData::RH), HS=ovec[ii](MeteoData::HS), RSWR=ovec[ii](MeteoData::RSWR);
			double ISWR=ovec[ii](MeteoData::ISWR);

			double albedo = .5;
			if (RSWR==IOUtils::nodata || ISWR==IOUtils::nodata || RSWR<=0 || ISWR<=0) {
				if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
					albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;

94
				if (ISWR==IOUtils::nodata && (RSWR!=IOUtils::nodata && HS!=IOUtils::nodata))
95
96
97
98
99
100
101
					ISWR = RSWR / albedo;
			} else {
				albedo = RSWR / ISWR;
				if (albedo>=1.) albedo=0.99;
				if (albedo<=0.) albedo=0.01;
			}

102
103
			const bool has_potRad = (ISWR!=IOUtils::nodata && TA!=IOUtils::nodata && RH!=IOUtils::nodata);
			if (has_potRad) 
104
				Sun.calculateRadiation(TA, RH, albedo);
105
106
			else 
				if (ovec[ii].date.getJulian(true) - julian_prev > 1.) continue; //no way to get ISWR and/or potential radiation, previous Md is too old
107
			
108
			const double Md = (has_potRad)? Sun.getSplitting(ISWR) : Md_prev; //fallback: use previous valid value
109
			tmp *= Md; //only keep the diffuse radiation, either on RSWR or ISWR
110
111
112
113
			if (has_potRad) {
				Md_prev = Md;
				julian_prev = ovec[ii].date.getJulian(true);
			}
114
115
116
117
118
		}
	}
}

//linear interpolation between the available points
119
double ProcShade::getMaskElevation(const std::vector< std::pair<double,double> > &mask, const double& azimuth) const
120
121
122
123
124
125
126
127
128
129
130
{
	const std::vector< std::pair<double, double> >::const_iterator next = std::upper_bound(mask.begin(), mask.end(), make_pair(azimuth, 0.), sort_pred()); //first element that is > azimuth
	
	double x1, y1, x2, y2;
	if (next!=mask.begin() && next!=mask.end()) { //normal case
		const size_t ii = next - mask.begin();
		x1 = mask[ii-1].first;
		y1 = mask[ii-1].second;
		x2 = mask[ii].first;
		y2 = mask[ii].second;
	} else {
131
132
133
134
		x1 = mask.back().first - 360.;
		y1 = mask.back().second;
		x2 = mask.front().first;
		y2 = mask.front().second;
135
136
137
138
139
140
141
142
143
144
145
146
147
148
	}
	
	const double a = (y2 - y1) / (x2 - x1);
	const double b = y2 - a * x2;
	
	return a*azimuth + b;
}

void ProcShade::readMask(const std::string& filter, const std::string& filename, std::vector< std::pair<double,double> > &o_mask)
{
	o_mask.clear();
	std::ifstream fin(filename.c_str());
	if (fin.fail()) {
		std::ostringstream ss;
149
		ss << "Filter " << filter << ": " << "error opening file \"" << filename << "\", possible reason: " << strerror(errno);
150
151
152
		throw AccessException(ss.str(), AT);
	}

153
	char eoln = FileUtils::getEoln(fin); //get the end of line character for the file
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192

	try {
		size_t lcount=0;
		double azimuth, value;
		do {
			lcount++;
			std::string line;
			getline(fin, line, eoln); //read complete line
			IOUtils::stripComments(line);
			IOUtils::trim(line);
			if (line.empty()) continue;

			std::istringstream iss(line);
			iss.setf(std::ios::fixed);
			iss.precision(std::numeric_limits<double>::digits10);
			iss >> std::skipws >> azimuth;
			if ( !iss || azimuth<0. || azimuth>360.) {
				std::ostringstream ss;
				ss << "Invalid azimuth in file " << filename << " at line " << lcount;
				throw InvalidArgumentException(ss.str(), AT);
			}
			iss >> std::skipws >> value;
			if ( !iss ){
				std::ostringstream ss;
				ss << "Invalid value in file " << filename << " at line " << lcount;
				throw InvalidArgumentException(ss.str(), AT);
			}

			const std::pair<double,double> tmp( azimuth, value );
			o_mask.push_back( tmp );
		} while (!fin.eof());
		fin.close();
	} catch (const std::exception&){
		if (fin.is_open()) {//close fin if open
			fin.close();
		}
		throw;
	}
	
193
	if (o_mask.empty()) throw InvalidArgumentException("In filter 'SHADE', no valid mask found in file '"+filename+"'", AT);
194
195
196
	std::sort(o_mask.begin(), o_mask.end(), sort_pred());
}

197
198
199
200
201
202
203
204
205
206
void ProcShade::computeMask(const DEMObject& i_dem, const StationData& sd, std::vector< std::pair<double,double> > &o_mask, const bool& dump_mask)
{
	o_mask.clear();
	Coords position( sd.position );
	if (!i_dem.gridify(position)) {
		const string msg = "In filter 'SHADE', station '"+sd.stationID+"' "+position.toString(Coords::LATLON)+" is not included in the DEM "+i_dem.llcorner.toString(Coords::LATLON);
		throw NoDataException(msg, AT);
	}
	
	i_dem.getHorizon(position, 10., o_mask); //by 10deg increments
207
208
	if (o_mask.empty()) throw InvalidArgumentException( "In filter 'SHADE', could not compute mask from DEM '"+i_dem.llcorner.toString(Coords::LATLON)+"'", AT);

209
210
211
212
213
214
215
216
	if (dump_mask) {
		std::cout << "Horizon mask for station '" << sd.stationID << "'\n";
		for (size_t ii=0; ii<o_mask.size(); ii++)
			std::cout << o_mask[ii].first << " " << o_mask[ii].second << "\n";
		std::cout << "\n";
	}
}

217
218
219
void ProcShade::parse_args(const std::vector<std::string>& vec_args)
{
	const size_t nrArgs = vec_args.size();
220
221
	if (nrArgs==0) { //compute from DEM
		IOHandler io(cfg);
222
		dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevations
223
		io.readDEM(dem);
224
	} else if (nrArgs==1) {
225
		const std::string root_path( cfg.getConfigRootDir() );
226
		//if this is a relative path, prefix the path with the current path
227
		const std::string in_filename( vec_args[0] );
228
		const std::string prefix = ( FileUtils::isAbsolutePath(in_filename) )? "" : root_path+"/";
229
230
		const std::string path( FileUtils::getPath(prefix+in_filename, true) );  //clean & resolve path
		const std::string filename( path + "/" + FileUtils::getFilename(in_filename) );
231
		std::vector< std::pair<double,double> > mask;
232
		readMask(getName(), filename, mask);
233
		masks["*"] = mask; //this mask is valid for ALL stations
234
235
236
237
238
239
	} else
		throw InvalidArgumentException("Wrong number of arguments for filter " + getName(), AT);

}

} //end namespace