The first algorithm implementing (Winstral, 2002) is now being tested. This...

The first algorithm implementing (Winstral, 2002) is now being tested. This seems to work but the mass balance still has to be validated. The IOManager can now be constructed directly with a io.ini file name. Some code clean up has been performed in InterpolationAlgorithms.

meteoio/IOManager.cc

@@ -22,6 +22,17 @@ using namespace std;
 
 namespace mio {
 
+IOManager::IOManager(const std::string& filename_in) : cfg(filename_in), rawio(cfg), bufferedio(rawio, cfg),
+                                            meteoprocessor(cfg), interpolator(cfg), dataGenerator(cfg),
+                                            v_params(), v_coords(), v_stations(),
+                                            proc_properties(), point_cache(), filtered_cache(),
+                                            fcache_start(Date(0.0, 0.)), fcache_end(Date(0.0, 0.)), //this should not matter, since 0 is still way back before any real data...
+                                            processing_level(IOManager::filtered | IOManager::resampled | IOManager::generated),
+                                            virtual_stations(false), skip_virtual_stations(false)
+{
+	initIOManager();
+}
+
 IOManager::IOManager(const Config& i_cfg) : cfg(i_cfg), rawio(cfg), bufferedio(rawio, cfg),
                                             meteoprocessor(cfg), interpolator(cfg), dataGenerator(cfg),
                                             v_params(), v_coords(), v_stations(),
@@ -29,6 +40,11 @@ IOManager::IOManager(const Config& i_cfg) : cfg(i_cfg), rawio(cfg), bufferedio(r
                                             fcache_start(Date(0.0, 0.)), fcache_end(Date(0.0, 0.)), //this should not matter, since 0 is still way back before any real data...
                                             processing_level(IOManager::filtered | IOManager::resampled | IOManager::generated),
                                             virtual_stations(false), skip_virtual_stations(false)
+{
+	initIOManager();
+}
+
+void IOManager::initIOManager()
 {
 	meteoprocessor.getWindowSize(proc_properties);
 	interpolator.setIOManager(*this); //because "*this" does not necessarily exist in the initialization list...

meteoio/IOManager.h

@@ -38,6 +38,7 @@ class IOManager {
 			num_of_levels = 1 << 4
 		};
 
+		IOManager(const std::string& filename_in);
 		IOManager(const Config& i_cfg);
 
 		//Legacy support to support functionality of the IOInterface superclass:
@@ -213,6 +214,7 @@ class IOManager {
 		void clear_cache();
 
 	private:
+		void initIOManager();
 		void initVirtualStations();
 		void fill_filtered_cache();
 		bool read_filtered_cache(const Date& start_date, const Date& end_date,

meteoio/InterpolationAlgorithms.cc

@@ -51,6 +51,8 @@ InterpolationAlgorithm* AlgorithmFactory::getAlgorithm(const std::string& i_algo
 		return new ILWRAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
 	} else if (algoname == "WIND_CURV"){// wind velocity interpolation (using a heuristic terrain effect)
 		return new SimpleWindInterpolationAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
+	} else if (algoname == "WINSTRAL"){// Winstral wind exposure factor
+		return new WinstralAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
 	} else if (algoname == "ODKRIG"){// ordinary kriging
 		return new OrdinaryKrigingAlgorithm(i_mi, i_vecArgs, i_algoname, iom);
 	} else if (algoname == "ODKRIG_LAPSE"){// ordinary kriging with lapse rate
@@ -369,7 +371,7 @@ LocalIDWLapseAlgorithm::LocalIDWLapseAlgorithm(Meteo2DInterpolator& i_mi, const
 	if (vecArgs.size() == 1) { //compute lapse rate on a reduced data set
 		IOUtils::convertString(nrOfNeighbors, vecArgs[0]);
 	} else { //incorrect arguments, throw an exception
-		throw InvalidArgumentException("Please provide the number of nearest neighbors to use for the IDW_LAPSE algorithm", AT);
+		throw InvalidArgumentException("Please provide the number of nearest neighbors to use for the "+algo+" algorithm", AT);
 	}
 }
 
@@ -567,11 +569,93 @@ void SimpleWindInterpolationAlgorithm::calculate(const DEMObject& dem, Grid2DObj
 }
 
 
+double WinstralAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
+{
+	//This algorithm is only valid for HNW (we could add HS later)
+	if (in_param!=MeteoData::HNW)
+		return 0.0;
+
+	date = i_date;
+	param = in_param;
+	nrOfMeasurments = getData(date, param, vecData, vecMeta);
+
+	if (nrOfMeasurments==0)
+		return 0.0;
+
+	return 0.9;
+}
+
+void WinstralAlgorithm::initGrid(const DEMObject& dem, Grid2DObject& grid)
+{
+	//retrieve optional arguments
+	std::string base_algo;
+	if (vecArgs.empty()){
+		base_algo=std::string("IDW_LAPSE");
+	} else if (vecArgs.size() == 1){
+		IOUtils::convertString(base_algo, vecArgs[0]);
+	} else { //incorrect arguments, throw an exception
+		throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
+	}
+
+	//initialize precipitation grid with user supplied algorithm (IDW_LAPSE by default)
+	IOUtils::toUpper(base_algo);
+	vector<string> vecArgs2;
+	mi.getArgumentsForAlgorithm(MeteoData::getParameterName(param), base_algo, vecArgs2);
+	auto_ptr<InterpolationAlgorithm> algorithm(AlgorithmFactory::getAlgorithm(base_algo, mi, vecArgs2, iomanager));
+	algorithm->getQualityRating(date, param);
+	algorithm->calculate(dem, grid);
+	info << algorithm->getInfo();
+}
+
+void WinstralAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
+{
+	info.clear(); info.str("");
+	initGrid(dem, grid);
+
+	//get synoptic wind direction
+	const double synoptic_bearing = 0.; //HACK
+	//two options:
+	// 1) locate the stations in DEM and check if they are higher than their surroundings within a given radius
+	// 2) simply compute a mean or median direction
+	// (2) can be used on all the stations selected in (1)
+
+	//compute wind exposure factor
+	Grid2DObject Sx;
+	Interpol2D::WinstralSX(dem, 300., synoptic_bearing, Sx);
+
+	//use the Sx to tweak the precipitation field
+	//get the scaling parameters
+	const double min_sx = Sx.grid2D.getMin(); //negative
+	const double max_sx = Sx.grid2D.getMax(); //positive
+	double sum_erosion=0., sum_deposition=0.;
+	//erosion: fully eroded at min_sx
+	for(size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
+		const double sx = Sx(ii);
+		double &val = grid(ii);
+		if (sx<0.) {
+			sum_erosion += val*sx/min_sx;
+			val *= 1. - sx/min_sx;
+		}
+		else sum_deposition += val;
+	}
+
+	//deposition: garantee mass balance conservation
+	const double ratio = sum_erosion/sum_deposition;
+	for(size_t ii=0; ii<Sx.getNx()*Sx.getNy(); ii++) {
+		const double sx = Sx(ii);
+		double &val = grid(ii);
+		if (sx>0.) {
+			val *= (1. + sx/max_sx) * ratio;
+		}
+	}
+}
+
+
 std::string USERInterpolation::getGridFileName() const
 {
 //HACK: use read2DGrid(grid, MeteoGrid::Parameters, Date) instead?
 	if (vecArgs.size() != 1){
-		throw InvalidArgumentException("Please provide the path to the grids for the USER interpolation algorithm", AT);
+		throw InvalidArgumentException("Please provide the path to the grids for the "+algo+" interpolation algorithm", AT);
 	}
 	const std::string ext(".asc");
 	const std::string& grid_path = vecArgs[0];
@@ -594,7 +678,7 @@ double USERInterpolation::getQualityRating(const Date& i_date, const MeteoData::
 	filename = getGridFileName();
 
 	if (!IOUtils::validFileName(filename)) {
-		cerr << "[E] Invalid grid filename for USER interpolation algorithm: " << filename << "\n";
+		cerr << "[E] Invalid grid filename for "+algo+" interpolation algorithm: " << filename << "\n";
 		return 0.0;
 	}
 	if(IOUtils::fileExists(filename)) {
@@ -628,7 +712,7 @@ double SnowHNWInterpolation::getQualityRating(const Date& i_date, const MeteoDat
 void SnowHNWInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
 {
 	info.clear(); info.str("");
-	info.clear(); info.str("");
+
 	//retrieve optional arguments
 	std::string base_algo;
 	if (vecArgs.empty()){
@@ -636,7 +720,7 @@ void SnowHNWInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
 	} else if (vecArgs.size() == 1){
 		IOUtils::convertString(base_algo, vecArgs[0]);
 	} else { //incorrect arguments, throw an exception
-		throw InvalidArgumentException("Wrong number of arguments supplied for the HNW_SNOW algorithm", AT);
+		throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
 	}
 
 	//initialize precipitation grid with user supplied algorithm (IDW_LAPSE by default)

meteoio/InterpolationAlgorithms.h

@@ -354,7 +354,7 @@ class ILWRAlgorithm : public InterpolationAlgorithm {
 
 /**
  * @class SimpleWindInterpolationAlgorithm
- * @brief Curvature/slope influenced  wind interpolation algorithm.
+ * @brief Curvature/slope influenced wind interpolation algorithm.
  * This is an implementation of the method described in (Liston & Elder, 2006): the wind speed and direction are
  * spatially interpolated using IDWLapseAlgorithm for the wind speed and using the user defined wind direction
  * interpolation algorithm. Then, the wind speed and direction fields are altered by wind weighting factors
@@ -373,6 +373,28 @@ class SimpleWindInterpolationAlgorithm : public InterpolationAlgorithm {
 		std::vector<double> vecDataVW, vecDataDW; ///<vectors of extracted VW and DW
 };
 
+/**
+ * @class WinstralAlgorithm
+ * @brief DEM-based wind-exposure interpolation algorithm.
+ * This is an implementation of the method described in Winstral, Elder, & Davis,
+ * <i>"Spatial snow modeling of wind-redistributed snow using terrain-based parameters"</i>, 2002,
+ * Journal of Hydrometeorology, <b>3(5)</b>, 524-538.
+ * The DEM is used to compute wind exposure factors that are used to alter the precipitation fields.
+ * It is usually a good idea to provide a DEM that also contain the accumulated snow height in order
+ * to get a progrfessive softening of the terrain features.
+ */
+class WinstralAlgorithm : public InterpolationAlgorithm {
+	public:
+		WinstralAlgorithm(Meteo2DInterpolator& i_mi,
+					const std::vector<std::string>& i_vecArgs,
+					const std::string& i_algo, IOManager& iom)
+			: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
+		virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
+		virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
+	private:
+		void initGrid(const DEMObject& dem, Grid2DObject& grid);
+};
+
 /**
  * @class USERInterpolation
  * @brief Reads user provided gridded data on the disk.

meteoio/meteostats/libinterpol2D.cc

@@ -24,6 +24,8 @@
 #include <meteoio/MathOptim.h> //math optimizations
 #include <meteoio/ResamplingAlgorithms2D.h> //for Winstral
 
+#include <meteoio/Timer.h> //HACK temporary for benchamrks
+
 using namespace std;
 
 namespace mio {
@@ -562,6 +564,8 @@ void Interpol2D::SteepSlopeRedistribution(const DEMObject& dem, const Grid2DObje
 //compute the Winstral sx factor for one single direction and one single point (ii,jj) in dem up to dmax distance
 double Interpol2D::WinstralSX_core(const Grid2DObject& dem, const double& dmax, const double& bearing, const size_t& i, const size_t& j)
 {
+	const double dmin = 20.; //cells closer than dmin don't play any role
+	const double inv_dmin = 1./dmin;
 	const double inv_dmax = 1./dmax;
 	const double alpha_rad = bearing*Cst::to_rad;
 	const double ref_altitude = dem(i, j);
@@ -569,8 +573,9 @@ double Interpol2D::WinstralSX_core(const Grid2DObject& dem, const double& dmax,
 	const int ii = static_cast<int>(i), jj = static_cast<int>(j);
 	const int ncols = static_cast<int>(dem.ncols), nrows = static_cast<int>(dem.nrows);
 
-	double max_tan_sx = 0.;
 	int ll=ii, mm=jj;
+
+	double max_tan_sx = 0.;
 	size_t nb_cells = 0;
 	while( !(ll<0 || ll>ncols-1 || mm<0 || mm>nrows-1) ) {
 		const double altitude = dem(ll, mm);
@@ -579,6 +584,7 @@ double Interpol2D::WinstralSX_core(const Grid2DObject& dem, const double& dmax,
 			//compute local sx
 			const double delta_elev = altitude - ref_altitude;
 			const double inv_distance = Optim::invSqrt( cellsize_sq*(Optim::pow2(ll-ii) + Optim::pow2(mm-jj)) );
+			if(inv_distance>inv_dmin) continue; //don't consider cells closer than dmin
 			if(inv_distance<inv_dmax) break; //stop if distance>dmax
 
 			const double tan_sx = delta_elev*inv_distance;
@@ -599,6 +605,8 @@ double Interpol2D::WinstralSX_core(const Grid2DObject& dem, const double& dmax,
 //Get the distance-weighted average of Sx for one single direction and one single point (ii,jj) in dem up to dmax distance
 double Interpol2D::AvgSX_core(const Grid2DObject& dem, const Grid2DObject& sx, const double& dmax, const double& bearing, const size_t& i, const size_t& j)
 {
+	const double dmin = 20.; //cells closer than dmin don't play any role
+	const double inv_dmin = 1./dmin;
 	const double inv_dmax = 1./dmax;
 	const double alpha_rad = bearing*Cst::to_rad;
 	const double cellsize_sq = Optim::pow2(dem.cellsize);
@@ -614,6 +622,7 @@ double Interpol2D::AvgSX_core(const Grid2DObject& dem, const Grid2DObject& sx, c
 		if(altitude==mio::IOUtils::nodata) continue; //jump over nodata cells
 		if( !(ll==ii && mm==jj) ) {
 			const double inv_distance = Optim::invSqrt( cellsize_sq*(Optim::pow2(ll-ii) + Optim::pow2(mm-jj)) );
+			if(inv_distance>inv_dmin) continue; //don't consider cells closer than dmin
 			if(inv_distance<inv_dmax) break; //stop if distance>dmax
 
 			sum_sx += sx(ll,mm);
@@ -697,15 +706,17 @@ void Interpol2D::Winstral_deposition(const DEMObject& dem, const double& dmax, c
 
 	//now remove deposition that happens too far from erosion
 	if(scale<1.) sd = ResamplingAlgorithms2D::BilinearResampling( sd, 1./scale ); //we should be back at sx dimensions
-	for(size_t jj = 0; jj<sd.nrows; jj++) {
-		for(size_t ii = 0; ii<sd.ncols; ii++) {
+
+	const size_t jj_max = min(sd.nrows, grid.nrows); //HACK
+	const size_t ii_max = min(sd.ncols, grid.ncols);
+	for(size_t jj = 0; jj<jj_max; jj++) {
+		for(size_t ii = 0; ii<ii_max; ii++) {
 			if(sx(ii,jj)<=0)
 				grid(ii,jj) = sx(ii,jj);
 			else if(sd(ii,jj)<0.) grid(ii,jj) = sx(ii,jj) - sd(ii,jj);
 			else grid(ii,jj) = 0.;
 		}
 	}
-
 }
 
 void Interpol2D::WinstralSB(const DEMObject& dem, const double& dmax, const double& sepdist, const double& in_bearing, Grid2DObject& grid)

meteoio/meteostats/libinterpol2D.h

@@ -64,6 +64,8 @@ class Interpol2D {
 		                      const std::vector<StationData>& vecStations,
 		                      const DEMObject& dem, const Fit1D& variogram, Grid2DObject& grid);
 
+		static void WinstralSX(const DEMObject& dem, const double& dmax, const double& in_bearing, Grid2DObject& grid);
+
 	private:
 		//generic functions
 		static double InvHorizontalDistance(const double& X1, const double& Y1, const double& X2, const double& Y2);
@@ -91,7 +93,6 @@ class Interpol2D {
 		static double WinstralSX_core(const Grid2DObject& dem, const double& dmax, const double& bearing, const size_t& ii, const size_t& jj);
 		static double AvgSX_core(const Grid2DObject& dem, const Grid2DObject& sx, const double& dmax, const double& bearing, const size_t& ii, const size_t& jj);
 
-		static void WinstralSX(const DEMObject& dem, const double& dmax, const double& in_bearing, Grid2DObject& grid);
 		static void WinstralSB(const DEMObject& dem, const double& dmax, const double& sepdist, const double& in_bearing, Grid2DObject& grid);
 		static void Winstral_deposition(const DEMObject& dem, const double& dmax, const double& in_bearing, Grid2DObject& grid);