Alpine3D

Alpine3D Svn Source Tree

Root/trunk/tools/pt_extract.cc

1#include <meteoio/MeteoIO.h>
2// compile with something like:
3// gcc pt_extract.cc -I ~/usr/include/ -o pt_extract -lmeteoio -L ~/usr/lib/ -lstdc++
4//Then run it without any arguments from the path where the hard-coded GRID2DPATH, etc work
5//Since it creates internally its own ini file, there is no need to have one on the disk.
6
7using namespace std;
8using namespace mio; //The MeteoIO namespace is called mio
9
10//edit this function to configure pt_extract
11void setINI(Config& cfg) {
12const string pts_file = "./input/surface-grids/arolla.poi";
13const string dem_file = "./input/surface-grids/arolla.dem";
14
15cfg.addKey("GRID2D", "Input", "ARC");
16cfg.addKey("GRID2DPATH", "Input", "./output/grids");
17cfg.addKey("DEM", "Input", "ARC");
18cfg.addKey("DEMFILE", "Input", dem_file);
19cfg.addKey("COORDSYS", "Input", "CH1903");
20cfg.addKey("TZ", "Input", "1");
21cfg.addKey("POI", "Input", "SMET");
22cfg.addKey("POIFILE", "Input", pts_file);
23cfg.addKey("METEO", "Output", "SMET");
24cfg.addKey("METEOPATH", "Output", "./output");
25cfg.addKey("TZ", "Output", "1");
26cfg.addKey("COORDSYS", "Output", "CH1903");
27}
28
29///////////////////////////// no need to edit anything here below ////////////////////////////
30
31Date getDate(const string& name, Config& cfg) {
32double TZ;
33cfg.getValue("TZ", "Input", TZ);
34Date date;
35if( IOUtils::convertString(date, name, TZ)==false ) {
36stringstream ss;
37ss << "Can not parse date in \"" << name << "\"";
38throw ConversionFailedException(ss.str(), AT);
39}
40return date;
41}
42
43void fillMeteo(MeteoData& meteo, const Grid2DObject& grid, const string& fname, const Coords& pt) {
44size_t beg = fname.find_last_of(".");
45const string ext = fname.substr(beg+1, fname.length()-beg);
46
47if(ext.compare("sdp")==0) meteo(MeteoData::HS) = grid.grid2D( pt.getGridI(), pt.getGridJ());
48if(ext.compare("tss")==0) meteo(MeteoData::TSS) = grid.grid2D( pt.getGridI(), pt.getGridJ());
49if(ext.compare("ta")==0) meteo(MeteoData::TA) = grid.grid2D( pt.getGridI(), pt.getGridJ());
50if(ext.compare("lwr")==0) meteo(MeteoData::ILWR) = grid.grid2D( pt.getGridI(), pt.getGridJ());
51if(ext.compare("swr")==0) meteo(MeteoData::ISWR) = grid.grid2D( pt.getGridI(), pt.getGridJ());
52if(ext.compare("soil")==0) meteo(MeteoData::TSG) = grid.grid2D( pt.getGridI(), pt.getGridJ());
53if(ext.compare("rh")==0) meteo(MeteoData::RH) = grid.grid2D( pt.getGridI(), pt.getGridJ());
54if(ext.compare("psum")==0) meteo(MeteoData::PSUM) = grid.grid2D( pt.getGridI(), pt.getGridJ());
55}
56
57int main(int /*argc*/, char** /*argv*/) {
58Config cfg;
59setINI(cfg);
60
61//create iomanager
62IOManager io(cfg);
63Grid2DObject grid;
64
65//get all files from directory
66const string path = cfg.get("GRID2DPATH", "Input");
67const string pattern = "20";
68list<string> dirlist;
69FileUtils::readDirectory(path, dirlist, pattern);
70dirlist.sort();
71
72//get special points and create meteo vectors
73vector<Coords> pts;
74io.readPOI(pts);
75DEMObject dem;
76io.readDEM(dem);
77io.read2DGrid(grid, *dirlist.begin());
78vector< vector<MeteoData> > vecvecMeteo;
79vector<StationData> vecStation;
80for(unsigned int i=0; i<pts.size(); i++) {
81grid.gridify(pts[i]);
82StationData station;
83stringstream sname, sid;
84sname << "Point_" << pts[i].getGridI() << "_" << pts[i].getGridJ();
85station.stationName = sname.str();
86sid << "pt_" << i+1;
87station.stationID = sid.str();
88station.position = pts[i];
89const double altitude = dem.grid2D( pts[i].getGridI(), pts[i].getGridJ() );
90station.position.setAltitude(altitude, false);
91vecStation.push_back(station);
92}
93vecvecMeteo.insert(vecvecMeteo.begin(), vecStation.size(), vector<MeteoData>()); //allocation for the vectors
94
95vector<Date> vecDate;
96list<string>::iterator itr;
97for( itr = dirlist.begin(); itr != dirlist.end(); itr++ ) {
98io.read2DGrid(grid, *itr); //itr does not contain the path
99const Date date = getDate(*itr, cfg);
100
101for(unsigned int i=0; i<pts.size(); i++) { //data extraction for 1 station at a time
102vector<MeteoData>& vecMeteo = vecvecMeteo[i];
103bool found=false;
104unsigned int j=0;
105for(; j<vecMeteo.size(); j++) {
106if(vecMeteo[j].date==date) {
107found=true;
108break;
109}
110}
111
112if(found==false) {
113MeteoData meteo;
114meteo.meta = vecStation[i];
115meteo.date = date;
116vecMeteo.push_back(meteo);
117}
118
119fillMeteo(vecMeteo[j], grid, *itr, pts[i]);
120}
121}
122
123const string outpath = cfg.get("METEOPATH", "Output");
124io.writeMeteoData(vecvecMeteo, outpath);
125
126return 0;
127}

Archive Download this file

Revision: HEAD