The DEMObject was not properly buffered (the update flag was not honored),...

The DEMObject was not properly buffered (the update flag was not honored), this has been fixed. The A3DIO plugin now implements readStationData. A few speed improvements have been brought to A3DIO and SMETIO by using the "reserve" method of the vectors (so that the vectors don't have to reallocate memory constantly).

configure.ac

@@ -50,6 +50,7 @@ then
 
 	AX_CFLAGS_GCC_OPTION(-g, OPTIM)
 	AX_CFLAGS_GCC_OPTION(-O3, OPTIM)
+	AX_CFLAGS_GCC_OPTION(-ftree-vectorize, OPTIM)
 	AC_SUBST(OPTIM)
 
 	AX_CXXFLAGS_GCC_OPTION(-combine, MAKE_OPTIM)

meteoio/A3DIO.cc

@@ -57,6 +57,7 @@ namespace mio {
  */
 
 const double A3DIO::plugin_nodata = -9999.0; //plugin specific nodata value
+const unsigned int A3DIO::buffer_reserve = 23*24*2; //kind of average size of a buffer for optimizing vectors
 
 //Main constructor
 A3DIO::A3DIO(const std::string& configfile) : IOInterface(NULL), cfg(configfile)
@@ -132,10 +133,40 @@ void A3DIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& data,
 	}
 }
 
-void A3DIO::readStationData(const Date&, std::vector<StationData>&)
+void A3DIO::readStationData(const Date& timestamp, std::vector<StationData>& vecStation)
 {
 	//Nothing so far
-	throw IOException("Nothing implemented here", AT);
+	//throw IOException("Nothing implemented here", AT);
+	vecStation.clear();
+
+	//read 1D station and put it into vecStation
+	std::string file_1d;
+	cfg.getValue("METEOPATH", "Input", file_1d);
+	file_1d += "/meteo1d.txt";
+
+	if (!IOUtils::fileExists(file_1d)) {
+		throw FileNotFoundException(file_1d, AT);
+	}
+
+	fin.clear();
+	fin.open (file_1d.c_str(), std::ifstream::in);
+	if (fin.fail()) {
+		throw FileAccessException(file_1d,AT);
+	}
+
+	StationData sd;
+	read1DStation(file_1d, sd);
+	vecStation.push_back(sd);
+	cleanup();
+
+	//read 2D stations and add them to vecStation
+	std::vector<StationData> tmpvecS;
+	read2DStations(timestamp, tmpvecS);
+	cleanup();
+
+	for(unsigned int i=0; i<tmpvecS.size(); i++) {
+		vecStation.push_back(tmpvecS[i]);
+	}
 }
 
 void A3DIO::readMeteoData(const Date& dateStart, const Date& dateEnd, std::vector< std::vector<MeteoData> >& vecMeteo)
@@ -166,7 +197,7 @@ void A3DIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
 		read2DMeteo(vecMeteo, vecStation);
 	} catch(std::exception& e){
 		std::cerr << "[E] No meteo2d data found or error while reading it, using only Meteo1D data: "
-			<< std::endl << "\t" << e.what() << std::endl;
+		          << std::endl << "\t" << e.what() << std::endl;
 	}
 }
 
@@ -186,36 +217,13 @@ void A3DIO::convertUnits(MeteoData& meteo)
 	}
 }
 
-void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd,
-				std::vector< std::vector<MeteoData> >& vecMeteo,
-				std::vector< std::vector<StationData> >& vecStation)
+void A3DIO::read1DStation(std::string& file_1d, StationData& sd)
 {
 	double latitude=IOUtils::nodata, longitude=IOUtils::nodata,
-		xcoord=IOUtils::nodata, ycoord=IOUtils::nodata, altitude=IOUtils::nodata;
-	std::string tmp="", line="";
-	Date tmp_date;
-	std::vector<std::string> tmpvec;
+	       xcoord=IOUtils::nodata, ycoord=IOUtils::nodata, altitude=IOUtils::nodata;
 	std::map<std::string, std::string> header; // A map to save key value pairs of the file header
-	MeteoData tmpdata;
-	StationData sd;
-	bool eofreached = false;
-
-	cfg.getValue("METEOPATH", "Input", tmp);
-	tmp += "/meteo1d.txt";
-
-	if (!IOUtils::fileExists(tmp)) {
-		throw FileNotFoundException(tmp, AT);
-	}
-
-	fin.clear();
-	fin.open (tmp.c_str(), std::ifstream::in);
-	if (fin.fail()) {
-		throw FileAccessException(tmp,AT);
-	}
-
-	char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
-
-	//Go through file, save key value pairs
+	
+	//read and parse the header
 	try {
 		//Read in station meta data
 		IOUtils::readKeyValueHeader(header, fin, 5, "="); //Read in 5 lines as header
@@ -239,19 +247,56 @@ void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd,
 		try {
 			location.check();
 		} catch(...) {
-			std::cerr << "[E] Error in geographic coordinates in file " << tmp << " trapped at " << AT << std::endl;
+			std::cerr << "[E] Error in geographic coordinates in file " << file_1d << " trapped at " << AT << std::endl;
 			throw;
 		}
 
-		sd.setStationData(location, "", "");
+		sd.setStationData(location, "meteo1d", "Meteo1D station");
+	} catch(...) {
+		std::cout << "[E] " << AT << ": "<< std::endl;
+		cleanup();
+		throw;
+	}
+}
+
+void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd,
+				std::vector< std::vector<MeteoData> >& vecMeteo,
+				std::vector< std::vector<StationData> >& vecStation)
+{
+	std::string file_1d="", line="";
+	Date tmp_date;
+	
+	MeteoData tmpdata;
+	StationData sd;
+	bool eofreached = false;
+
+	cfg.getValue("METEOPATH", "Input", file_1d);
+	file_1d += "/meteo1d.txt";
+
+	if (!IOUtils::fileExists(file_1d)) {
+		throw FileNotFoundException(file_1d, AT);
+	}
+
+	fin.clear();
+	fin.open (file_1d.c_str(), std::ifstream::in);
+	if (fin.fail()) {
+		throw FileAccessException(file_1d,AT);
+	}
+
+	char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
+
+	//get station metadata
+	read1DStation(file_1d, sd);
 
+	//Go through file, save key value pairs
+	try {
 		//Read one line, construct Date object and see whether date is greater or equal than the date_in object
 		IOUtils::skipLines(fin, 1, eoln); //skip rest of line
 
 		//Loop going through the data sequentially until dateStart is found
 		do {
 			getline(fin, line, eoln); //read complete line
-			eofreached = readMeteoDataLine(line, tmpdata, tmp);
+			eofreached = readMeteoDataLine(line, tmpdata, file_1d);
 			//tmpdata.cleanData();
 			convertUnits(tmpdata);
 
@@ -265,6 +310,8 @@ void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd,
 		} else if ((tmpdata.date <= dateEnd)  && (!eofreached)) {
 			vecMeteo.push_back( std::vector<MeteoData>() );
 			vecStation.push_back( std::vector<StationData>() );
+			vecMeteo[0].reserve(buffer_reserve);
+			vecStation[0].reserve(buffer_reserve);
 		}
 
 		while ((tmpdata.date <= dateEnd)  && (!eofreached)) {
@@ -273,7 +320,7 @@ void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd,
 			vecStation[0].push_back(sd);
 
 			getline(fin, line, eoln); //read complete line
-			eofreached = readMeteoDataLine(line, tmpdata, tmp);
+			eofreached = readMeteoDataLine(line, tmpdata, file_1d);
 			//tmpdata.cleanData();
 			convertUnits(tmpdata);
 		}
@@ -292,6 +339,7 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
 	Date tmp_date;
 	int tmp_ymdh[4];
 	std::vector<std::string> tmpvec;
+	tmpvec.reserve(6);
 	double tmp_values[6];
 
 	if (IOUtils::readLineToVec(line, tmpvec) != 10) {
@@ -326,6 +374,29 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
 }
 
 
+void A3DIO::read2DStations(const Date& timestamp, std::vector<StationData>& vecStation)
+{
+	unsigned int stations=0;
+	std::vector<std::string> filenames = std::vector<std::string>();
+	std::map<std::string, unsigned int> hashStations = std::map<std::string, unsigned int>();
+	
+	constructMeteo2DFilenames(timestamp, timestamp, filenames);
+	stations = getNrOfStations(filenames, hashStations);
+	vecStation.resize(stations); //we receive an empty vector, so we need to allocate the room it needs
+
+	try {
+		for (unsigned int ii=0; ii<filenames.size(); ii++){
+			read2DMeteoHeader(filenames[ii], hashStations, vecStation);
+		}
+	} catch(...) {
+		//clear all 2D meteo data if error occurs
+		if (vecStation.size() > 1)
+			vecStation.erase(vecStation.begin()+1, vecStation.end());
+		cleanup();
+		throw;
+	}
+}
+
 /*
   Preamble: Files are in METEOFILE directory. 4 types of files:
   prec????.txt == hnw
@@ -374,6 +445,8 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
 		for (unsigned int jj=0; jj<tmpvecS.size(); jj++){
 			vecMeteo.push_back( std::vector<MeteoData>() );
 			vecStation.push_back( std::vector<StationData>() );
+			vecStation[jj+1].reserve(buffer_reserve);
+			vecMeteo[jj+1].reserve(buffer_reserve);
 			for (unsigned int ii=0; ii<vecMeteo[0].size(); ii++){
 				//NOTE: there needs to be the same amount of 1D and 2D data
 				vecStation[jj+1].push_back(tmpvecS[jj]);

meteoio/A3DIO.h

@@ -68,8 +68,11 @@ class A3DIO : public IOInterface {
 
 	private:
 		static const double plugin_nodata; //plugin specific nodata value, e.g. -9999
+		static const unsigned int buffer_reserve; //for optimizing vectors (ie: avoid unecessary resizing)
+		void read1DStation(std::string& file_1d, StationData& sd);
 		void read1DMeteo(const Date& dateStart, const Date& dateEnd, 
 		                 std::vector< std::vector<MeteoData> >&, std::vector< std::vector<StationData> >&); ///< No buffering
+		void read2DStations(const Date& timestamp, std::vector<StationData>& vecStation);
 		void read2DMeteo(std::vector< std::vector<MeteoData> >&, std::vector< std::vector<StationData> >&); ///< No buffering
 
 		void constructMeteo2DFilenames(const Date& _startDate, const Date& _endDate, std::vector<std::string>& _filenames);

meteoio/BufferedIOHandler.cc

@@ -58,12 +58,22 @@ void BufferedIOHandler::read2DGrid(Grid2DObject& _grid2Dobj, const std::string&
 void BufferedIOHandler::readDEM(DEMObject& _grid2Dobj)
 {
 	std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:DEM");
-	if (it != mapBufferedGrids.end()) { //already in map
-		_grid2Dobj = (*it).second; 
+	if (it != mapBufferedGrids.end()) {
+		//already in map. If the update properties have changed,
+		//we copy the ones given in input and force the update of the object
+		const DEMObject::update_type in_ppt = (DEMObject::update_type)_grid2Dobj.getUpdatePpt();
+		_grid2Dobj = (*it).second;
+		const DEMObject::update_type buff_ppt = (DEMObject::update_type)_grid2Dobj.getUpdatePpt();
+		if(in_ppt!=buff_ppt) {
+			_grid2Dobj.setUpdatePpt(in_ppt);
+			_grid2Dobj.update();
+		}
 		return;
 	}
 	
 	DEMObject tmpgrid2D;
+	 //copy the updating policy of the destination
+	tmpgrid2D.setUpdatePpt((DEMObject::update_type)_grid2Dobj.getUpdatePpt());
 	iohandler.readDEM(tmpgrid2D);
 	mapBufferedGrids["/:DEM"] = tmpgrid2D;
 	_grid2Dobj = tmpgrid2D;

meteoio/Coords.cc

@@ -558,6 +558,7 @@ short int Coords::getEPSG() const {
 * @param epsg epsg code
 */
 void Coords::setEPSG(const short int epsg) {
+//TODO: get rid of the zone letter. This is not part of the standard and redundant (and messy)
 	bool found=false;
 	std::string coord_sys, coord_param;
 

meteoio/plugins/SMETIO.cc

@@ -47,6 +47,7 @@ namespace mio {
  */
 
 const std::string SMETIO::smet_version = "0.99";
+const unsigned int SMETIO::buffer_reserve = 23*24*2; //kind of average size of a buffer for optimizing vectors
 map<string, MeteoData::Parameters> SMETIO::mapParameterByName;
 const bool SMETIO::__init = SMETIO::initStaticData();
 
@@ -188,11 +189,11 @@ void SMETIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_ou
 void SMETIO::readStationData(const Date&, std::vector<StationData>& vecStation)
 {//big HACK: this is a barbaric code duplication!! Plus it should support coordinates in the data
 //ie: it should use the given date!
-	unsigned int startindex=0, endindex=vecFiles.size();
 	vecStation.clear();
+	vecStation.reserve(nr_stations);
 
 	//Now loop through all requested stations, open the respective files and parse them
-	for (unsigned int ii=startindex; ii<endindex; ii++){
+	for (unsigned int ii=0; ii<nr_stations; ii++){
 		bool isAscii = true;
 		string filename = vecFiles.at(ii); //filename of current station
 		
@@ -257,6 +258,8 @@ void SMETIO::parseInputOutputSection()
 		counter++;
 	} while (filename != "");
 
+	nr_stations = counter - 1;
+
 	//Parse output section: extract info on whether to write ASCII or BINARY format, gzipped or not
 	outpath = "";
 	outputIsAscii = true; 
@@ -314,6 +317,8 @@ void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
 		
 		vecMeteo = vector< vector<MeteoData> >(vecFiles.size());
 		vecStation = vector< vector<StationData> >(vecFiles.size());
+		vecMeteo.reserve(nr_stations);
+		vecStation.reserve(nr_stations);
 	}
 
 	//Now loop through all requested stations, open the respective files and parse them
@@ -372,12 +377,15 @@ void SMETIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
 }
 
 void SMETIO::readDataBinary(const char&, const std::string&, const double& timezone,
-					    const StationData& sd, const std::vector<std::string>& vecDataSequence,
-					    const Date& dateStart, const Date& dateEnd,
-					    std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
+                            const StationData& sd, const std::vector<std::string>& vecDataSequence,
+                            const Date& dateStart, const Date& dateEnd,
+                            std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
 {
 	const unsigned int nrOfColumns = vecDataSequence.size();
 
+	vecMeteo.reserve(buffer_reserve);
+	vecStation.reserve(buffer_reserve);
+
 	while (!fin.eof()){
 		MeteoData md;
 		StationData tmpsd = sd;
@@ -437,14 +445,18 @@ void SMETIO::readDataBinary(const char&, const std::string&, const double& timez
 }
 
 void SMETIO::readDataAscii(const char& eoln, const std::string& filename, const double& timezone,
-					   const StationData& sd, const std::vector<std::string>& vecDataSequence,
-					   const Date& dateStart, const Date& dateEnd,
-					   std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
+                           const StationData& sd, const std::vector<std::string>& vecDataSequence,
+                           const Date& dateStart, const Date& dateEnd,
+                           std::vector<MeteoData>& vecMeteo, std::vector<StationData>& vecStation)
 {
 	string line = "";
 	vector<string> tmpvec;
 	const unsigned int nrOfColumns = vecDataSequence.size();
 
+	tmpvec.reserve(nrOfColumns);
+	vecMeteo.reserve(buffer_reserve);
+	vecStation.reserve(buffer_reserve);
+
 	while (!fin.eof()){
 		getline(fin, line, eoln);
 		IOUtils::stripComments(line);
@@ -542,6 +554,7 @@ void SMETIO::readHeader(const char& eoln, const std::string& filename, bool& loc
 	if(epsg!=IOUtils::snodata) {
 		sd.position.setEPSG(epsg);
 	}
+	
 	IOUtils::getValueForKey(mapHeader, "easting", easting, IOUtils::nothrow);
 	if (easting != IOUtils::nodata){ //HACK
 		IOUtils::getValueForKey(mapHeader, "northing", northing);

meteoio/plugins/SMETIO.h

@@ -62,6 +62,7 @@ class SMETIO : public IOInterface {
 
 	private:
 		static const std::string smet_version;
+		static const unsigned int buffer_reserve; //for optimizing vectors (ie: avoid unecessary resizing)
 		static std::map<std::string, MeteoData::Parameters> mapParameterByName; ///<Associate name and meteo parameter
 		static const bool __init;    ///<helper variable to enable the init of static collection data
 		static bool initStaticData();///<initialize the static map meteoparamname 
@@ -94,6 +95,7 @@ class SMETIO : public IOInterface {
 		void checkSignature(const std::vector<std::string>& vecSignature, const std::string& filename, bool& isAscii);
 		void setFormatting(const MeteoData::Parameters& paramindex);
 
+		unsigned int nr_stations; //number of stations to read from
 		std::vector<std::string> vecFiles;  //read from the Config [Input] section
 		std::string outpath;                //read from the Config [Output] section
 		bool outputIsAscii, outputIsGzipped;//read from the Config [Output] section