The original PrecipSNOW spatial interpolation algorithm has been reimplemented...

The original PrecipSNOW spatial interpolation algorithm has been reimplemented since the version implemented by Rob is unreliable

meteoio/InterpolationAlgorithms.cc

@@ -989,12 +989,11 @@ void SnowHNWInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
 	info.clear(); info.str("");
 
 	//retrieve optional arguments
-	std::string base_algo;
-	if (vecArgs.empty()){
-		base_algo=std::string("IDW_LAPSE");
-	} else if (vecArgs.size() == 1){
+	std::string base_algo("IDW_LAPSE");
+	const size_t nrArgs = vecArgs.size();
+	if (nrArgs == 1){
 		IOUtils::convertString(base_algo, vecArgs[0]);
-	} else { //incorrect arguments, throw an exception
+	} else if (nrArgs>1){ //incorrect arguments, throw an exception
 		throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
 	}
 
@@ -1012,7 +1011,8 @@ void SnowHNWInterpolation::calculate(const DEMObject& dem, Grid2DObject& grid)
 	mi.interpolate(date, dem, MeteoData::TA, ta);
 
 	//slope/curvature correction for solid precipitation
-	Interpol2D::SteepSlopeRedistribution(dem, ta, grid);
+	Interpol2D::PrecipSnow(dem, ta, grid);
+	//Interpol2D::SteepSlopeRedistribution(dem, ta, grid);
 	//Interpol2D::CurvatureCorrection(dem, ta, grid);
 }
 

meteoio/meteoStats/libinterpol2D.cc

@@ -555,6 +555,56 @@ void Interpol2D::SteepSlopeRedistribution(const DEMObject& dem, const Grid2DObje
 	}
 }
 
+/**
+* @brief Distribute precipitation in a way that reflects snow redistribution on the ground, according to (Huss, 2008)
+* This method modifies the solid precipitation distribution according to the local slope and curvature. See
+* <i>"Quantitative evaluation of different hydrological modelling approaches in a partly glacierized Swiss watershed"</i>, Magnusson et All., Hydrological Processes, 2010, under review.
+* and
+* <i>"Modelling runoff from highly glacierized alpine catchments in a changing climate"</i>, Huss et All., Hydrological Processes, <b>22</b>, 3888-3902, 2008.
+* @param dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis.
+* @param ta array of air temperatures used to determine if precipitation is rain or snow
+* @param grid 2D array of precipitation to fill
+* @author Florian Kobierska, Jan Magnusson and Mathias Bavay
+*/
+void Interpol2D::PrecipSnow(const DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid)
+{
+	if(!grid.isSameGeolocalization(dem)) {
+		throw IOException("Requested grid does not match the geolocalization of the DEM", AT);
+	}
+	const double dem_max_curvature=dem.max_curvature, dem_range_curvature=(dem.max_curvature-dem.min_curvature);
+
+	const size_t nrows = grid.getNy();
+	const size_t ncols = grid.getNx();
+	
+	for (size_t j=0;j<nrows;j++) {
+		for (size_t i=0;i<ncols;i++) {
+			// Get input data
+			const double slope = dem.slope(i, j);
+			const double curvature = dem.curvature(i, j);
+			double val = grid.grid2D(i, j);
+
+			if(ta.grid2D(i, j)<=Cst::t_water_freezing_pt) {
+				//we only modify the grid of precipitations if air temperature
+				//at this point is below or at freezing
+				if(slope==IOUtils::nodata || curvature==IOUtils::nodata) {
+					val = IOUtils::nodata;
+				} else if (slope>60.) {
+					//No snow precipitation happens for these slopes
+					val = 0.;
+				} else if (slope>40.) {
+					//Linear transition from no snow to 100% snow
+					val *= (60.-slope)/20.;
+				} //else: unchanged
+
+				if(val!=IOUtils::nodata && dem_range_curvature!=0.) {
+					//cf Huss
+					grid.grid2D(i, j) = val*(0.5-(curvature-dem_max_curvature)/dem_range_curvature);
+				}
+			}
+		}
+	}
+}
+
 //Compute the wind direction changes by the terrain, see Ryan, "a mathematical model for diagnosis
 //and prediction of surface winds in mountainous terrain", 1977, journal of applied meteorology, 16, 6
 /**

meteoio/meteoStats/libinterpol2D.h

@@ -60,6 +60,7 @@ class Interpol2D {
 		static void ListonWind(const DEMObject& i_dem, Grid2DObject& VW, Grid2DObject& DW);
 		static void CurvatureCorrection(DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid);
 		static void SteepSlopeRedistribution(const DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid);
+		static void PrecipSnow(const DEMObject& dem, const Grid2DObject& ta, Grid2DObject& grid);
 		static void ODKriging(const std::vector<double>& vecData,
 		                      const std::vector<StationData>& vecStations,
 		                      const DEMObject& dem, const Fit1D& variogram, Grid2DObject& grid);
@@ -94,7 +95,7 @@ class Interpol2D {
 
 		static void steepestDescentDisplacement(const DEMObject& dem, const Grid2DObject& grid, const size_t& ii, const size_t& jj, short &d_i_dest, short &d_j_dest);
 		static double depositAroundCell(const DEMObject& dem, const size_t& ii, const size_t& jj, const double& precip, Grid2DObject &grid);
-
+		
 		static void WinstralSX(const DEMObject& dem, const double& dmax, const double& in_bearing, Grid2DObject& grid);
 
 		//weighting methods