Code tagging and some usefull methods for debugging and the begining of a true...

Code tagging and some usefull methods for debugging and the begining of a true accumulator for time series (for any sampling rates).

meteoio/Meteo1DInterpolator.cc

@@ -108,4 +108,19 @@ string Meteo1DInterpolator::getInterpolationForParameter(const std::string& parn
 	return "linear"; //the default resampling is linear
 }
 
+std::ostream& operator<<(std::ostream& os, const Meteo1DInterpolator& Interpolator) {
+
+	os << "<Meteo1DInterpolator>\n";
+	for (unsigned int jj=0; jj<Interpolator.tasklist.size(); jj++){
+		os << MeteoData::getParameterName(jj) << "::" << Interpolator.tasklist[jj] << "\t";
+		for (unsigned int ii=0; ii<Interpolator.taskargs[jj].size(); ii++){
+			os << "ARGS: " << Interpolator.taskargs[jj][ii] << " ";
+		}
+		os << "\n";
+	}
+	os << "</Meteo1DInterpolator>\n";
+
+	return os;
+}
+
 } //namespace

meteoio/Meteo1DInterpolator.h

@@ -62,6 +62,8 @@ class Meteo1DInterpolator {
 		 */
 		unsigned int resampleData(const Date& date, std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
 
+		friend std::ostream& operator<<(std::ostream& os, const Meteo1DInterpolator& Interpolator);
+
  	private:
 		std::string getInterpolationForParameter(const std::string& parname, std::vector<std::string>& vecArguments);
 

meteoio/Meteo2DInterpolator.cc

@@ -145,14 +145,15 @@ void Meteo2DInterpolator::checkMinMax(const double& minval, const double& maxval
 
 	gridobj.grid2D.size(nx, ny);
 
-	for (unsigned int ii=0; ii<nx; ii++){
-		for (unsigned int jj=0; jj<ny; jj++){
-			if (gridobj.grid2D(ii,jj) == IOUtils::nodata){
+	for (unsigned int jj=0; jj<ny; jj++){
+		for (unsigned int ii=0; ii<nx; ii++){
+			double& value = gridobj.grid2D(ii,jj);
+			if (value == IOUtils::nodata){
 				//Do nothing
-			} else if (gridobj.grid2D(ii,jj) < minval){
-				gridobj.grid2D(ii,jj) = minval;
-			} else if (gridobj.grid2D(ii,jj) > maxval){
-				gridobj.grid2D(ii,jj) = maxval;
+			} else if (value < minval){
+				value = minval;
+			} else if (value > maxval){
+				value = maxval;
 			}
 		}
 	}

meteoio/MeteoProcessor.cc

@@ -30,7 +30,7 @@ void MeteoProcessor::processData(const Date& date, const std::vector<MeteoData>&
 	unsigned int startindex = IOUtils::npos, endindex = IOUtils::npos;
 	
 	//No need to operate on the raw data, a copy of relevant data will be stored in these vectors:
-	std::vector<MeteoData> vecWindowM;   
+	std::vector<MeteoData> vecWindowM;
 	std::vector<StationData> vecWindowS;
 
 	/*

meteoio/ResamplingAlgorithms.cc

@@ -55,6 +55,7 @@ bool ResamplingAlgorithms::initStaticData()
 {
 	algorithmMap["linear"]            = &ResamplingAlgorithms::LinearResampling;
 	algorithmMap["nearest_neighbour"] = &ResamplingAlgorithms::NearestNeighbour;
+	algorithmMap["accumulate"] = &ResamplingAlgorithms::Accumulate;
 
 	return true;
 }
@@ -87,10 +88,9 @@ const resamplingptr& ResamplingAlgorithms::getAlgorithm(const std::string& algon
  */
 
 void ResamplingAlgorithms::NearestNeighbour(const unsigned int& pos, const MeteoData::Parameters& paramindex,
-                                            const std::vector<std::string>& taskargs,
-                                            std::vector<MeteoData>& vecM, std::vector<StationData>& vecS)
+                                            const std::vector<std::string>& /*taskargs*/,
+                                            std::vector<MeteoData>& vecM, std::vector<StationData>& /*vecS*/)
 {
-	(void)taskargs; (void)vecS;
 	if (pos >= vecM.size())
 		throw IOException("The position of the resampled element is out of bounds", AT);
 
@@ -148,10 +148,9 @@ void ResamplingAlgorithms::NearestNeighbour(const unsigned int& pos, const Meteo
  * @endcode
  */
 void ResamplingAlgorithms::LinearResampling(const unsigned int& pos, const MeteoData::Parameters& paramindex,
-								    const std::vector<std::string>& taskargs,
-								    std::vector<MeteoData>& vecM, std::vector<StationData>& vecS)
+                                            const std::vector<std::string>& taskargs,
+                                            std::vector<MeteoData>& vecM, std::vector<StationData>& /*vecS*/)
 {
-	(void)vecS;
 	if (pos >= vecM.size())
 		throw IOException("The position of the resampled element is out of bounds", AT);
 
@@ -218,9 +217,92 @@ void ResamplingAlgorithms::LinearResampling(const unsigned int& pos, const Meteo
 	const double& val1 = vecM[indexP1].param(paramindex);
 	const double& val2 = vecM[indexP2].param(paramindex);
 	vecM[pos].param(paramindex) = Interpol1D::linearInterpolation(vecM[indexP1].date.getJulianDate(), val1,
-										                 vecM[indexP2].date.getJulianDate(), val2, 
-										                 vecM[pos].date.getJulianDate());
+	                              vecM[indexP2].date.getJulianDate(), val2,
+	                              vecM[pos].date.getJulianDate());
 }
 
+/**
+ * @brief Accumulation over a user given period.
+ * The input data is accumulated over a given time interval (given as filter argument, in minutes).
+ * This is for example needed for converting rain gauges measurements read every 10 minutes to
+ * hourly precipitation measurements. Remarks:
+ * - the accumulation period has to be provided as an argument (in seconds)
+ * @code
+ * HNW::filter1 = accumulate
+ * HNW::arg1	 = 3600
+ * @endcode
+ */
+void ResamplingAlgorithms::Accumulate(const unsigned int& pos, const MeteoData::Parameters& paramindex,
+                                      const std::vector<std::string>& taskargs,
+                                      std::vector<MeteoData>& vecM, std::vector<StationData>& vecS)
+{
+	//Not ready for usage yet...
+	throw IOException("Not finished...", AT);
+
+	unsigned int npos = pos;
+	if (pos >= vecM.size())
+		throw IOException("The position of the resampled element is out of bounds", AT);
+
+	//Get accumulation period
+	double accumulate_period;
+	if (taskargs.size()==1) {
+		IOUtils::convertString(accumulate_period, taskargs[0]);
+		if(accumulate_period<=0.) {
+			std::stringstream tmp;
+			tmp << "Invalid accumulation period (" << accumulate_period << ") ";
+			tmp << "for parameter " << MeteoData::getParameterName(paramindex);
+			throw InvalidArgumentException(tmp.str(), AT);
+		}
+	} else {
+		std::stringstream tmp;
+		tmp << "Please provide accumulation period (in seconds) for parameter " << MeteoData::getParameterName(paramindex);
+		throw InvalidArgumentException(tmp.str(), AT);
+	}
+	
+	//find start of accumulation period
+	bool found_start=false;
+	int ii=pos;
+	const Date Start(vecM[pos].date.getJulianDate() - accumulate_period/(24.*3600.));
+	for (; ii>=0; ii--) {
+		if(vecM[(unsigned int)ii].date<=Start) {
+			found_start=true;
+			break;
+		}
+	}
+
+	if(found_start==false) {
+		std::cerr << "[W] Could not accumulate " << MeteoData::getParameterName(paramindex);
+		std::cerr << ", not enough data for accumulation period\n";
+		return;
+	}
+
+	//resample the starting point
+	std::vector<std::string> dummy_args;
+	const unsigned int org_size = vecM.size();
+	LinearResampling(ii, paramindex, dummy_args, vecM, vecS);
+	if(vecM.size()>org_size)
+		npos++;
+
+	//start accumulating
+	double sum = 0.0;
+	bool exist=false;
+	while(vecM[ii].date<=vecM[npos].date) {
+		//we will at least enter once into this loop since accumulate_period>0
+		const double& val = vecM[npos].param(paramindex);
+		if (val != IOUtils::nodata) {
+			sum += val;
+			exist=true;
+		}
+		ii++;
+	}
+
+	//add last remaining point after interpolating it
+	//HACK: TODO
+
+
+	//check if at least one point has been summed. If not -> nodata
+}
+
+
 } //namespace
 

meteoio/ResamplingAlgorithms.h

@@ -48,11 +48,15 @@ class ResamplingAlgorithms {
 
 		//Available algorithms
 		static void LinearResampling(const unsigned int& position, const MeteoData::Parameters& paramindex,
-                                       const std::vector<std::string>& taskargs,
-                                       std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
+		                             const std::vector<std::string>& taskargs,
+		                             std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
 		static void NearestNeighbour(const unsigned int& position, const MeteoData::Parameters& paramindex, 
-                                       const std::vector<std::string>& taskargs,
-                                       std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
+		                             const std::vector<std::string>& taskargs,
+		                             std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
+
+		static void Accumulate(const unsigned int& pos, const MeteoData::Parameters& paramindex,
+		                       const std::vector<std::string>& taskargs,
+		                       std::vector<MeteoData>& vecM, std::vector<StationData>& vecS);
 		
  	private:
 		static std::map<std::string, resamplingptr> algorithmMap;

meteoio/marshal_meteoio.cc

@@ -313,7 +313,7 @@ void marshal_DOUBLE2D(POPBuffer &buf, DOUBLE2D &data,int maxsize, int flag, POPM
 		buf.Pack(&nx,1);
 		buf.Pack(&ny,1);
 		if (nx>0 && ny>0) {
-			for (unsigned int i=0;i<nx;i++) buf.Pack(&data(i,0),ny);
+			for (unsigned int i=0;i<ny;i++) buf.Pack(&data(0,i), nx);
 		}
 	} else {
 		unsigned int nx,ny;
@@ -321,14 +321,14 @@ void marshal_DOUBLE2D(POPBuffer &buf, DOUBLE2D &data,int maxsize, int flag, POPM
 		buf.UnPack(&ny,1);
 		if (nx>0 && ny>0) {
 			data.resize(nx,ny);
-			for (unsigned int i=0;i<nx;i++) buf.UnPack(&data(i,0),ny);
+			for (unsigned int i=0;i<ny;i++) buf.UnPack(&data(0,i),nx);
 		} else
 			data.clear();
 	}
 }
 
 void marshal_DOUBLE3D(POPBuffer &buf, DOUBLE3D &data,int maxsize, int flag, POPMemspool *temp)
-{
+{//HACK: this marshalling is sub-optimal!!
 	(void)maxsize;
 	(void)*temp;
 	if (flag & FLAG_MARSHAL) {
@@ -362,22 +362,20 @@ void marshal_INT2D(POPBuffer &buf, INT2D &data,int maxsize, int flag, POPMemspoo
 	(void)maxsize;
 	(void)*temp;
 	if (flag & FLAG_MARSHAL) {
-		unsigned int dim[2];
-		data.size(dim[0],dim[1]);
-		buf.Pack(dim,2);
-		unsigned int nx=dim[0];
-		unsigned int ny=dim[1];
+		unsigned int nx, ny;
+		data.size(nx,ny);
+		buf.Pack(&nx,1);
+		buf.Pack(&ny,1);
 		if (nx>0 && ny>0) {
-			for (unsigned int i=0;i<nx;i++) buf.Pack(&data(i,0),ny);
+			for (unsigned int i=0;i<ny;i++) buf.Pack(&data(0,i),nx);
 		}
 	} else {
-		unsigned int dim[2];
-		buf.UnPack(dim,2);
-		unsigned int nx=dim[0];
-		unsigned int ny=dim[1];
+		unsigned int nx,ny;
+		buf.UnPack(&nx,1);
+		buf.UnPack(&ny,1);
 		if (nx>0 && ny>0) {
 			data.resize(nx,ny);
-			for (unsigned int i=0;i<nx;i++) buf.UnPack(&data(i,0),ny);
+			for (unsigned int i=0;i<ny;i++) buf.UnPack(&data(0,i),nx);
 		} else
 			data.clear();
 	}

meteoio/plugins/ARPSIO.cc

@@ -348,6 +348,7 @@ void ARPSIO::readGridLayer(const std::string& parameter, const unsigned int& lay
 			}
 	}
 
+//HACK TODO swap ix/iy (for speed) and number cells from llcorner (cf ARCIO)
 	//read the data we are interested in
 	for (unsigned int ix = 0; ix < dimx; ix++) {
 		for (unsigned int iy = 0; iy < dimy; iy++) {