Alpine3D

Alpine3D Svn Source Tree

Root/trunk/tools/glaciers.cc

1/* This is a test program to run single time steps simulations of the katabatic flow temperature correction
2 * compile with something like:
3 * gcc glaciers.cc -g -I ~/usr/include/ -o glaciers -lmeteoio -lalpine3d -L ~/usr/lib/ -lstdc++
4 */
5
6#include <meteoio/MeteoIO.h>
7#include <alpine3d/Alpine3D.h>
8
9using namespace mio;
10
11const Grid2DObject computeGlacierMap(const Grid2DObject& lus)
12{
13Grid2DObject map(lus); //only glaciated pixels should be IOUtils::nodata
14for (size_t ii=0; ii<lus.size(); ii++) {
15if (map(ii)==IOUtils::nodata) {
16map(ii) = 0;
17continue;
18}
19const int lus_code = ((int)(map(ii) + 0.0001) - 10000) / 100;
20if (lus_code==14)
21map(ii) = IOUtils::nodata;
22}
23
24return map;
25}
26
27void setPOI(const DEMObject& dem, std::vector<Coords>& pts, std::vector<StationData> &vecStation)
28{
29vecStation.clear();
30
31for(unsigned int i=0; i<pts.size(); i++) {
32dem.gridify(pts[i]);
33StationData station;
34std::stringstream sname, sid;
35sname << "Point_" << pts[i].getGridI() << "_" << pts[i].getGridJ();
36station.stationName = sname.str();
37sid << "pt_" << i+1;
38station.stationID = sid.str();
39station.position = pts[i];
40const double altitude = dem.grid2D( pts[i].getGridI(), pts[i].getGridJ() );
41station.position.setAltitude(altitude, false);
42vecStation.push_back(station);
43}
44}
45
46int main(int /*argc*/, char** argv) {
47Config cfg( "io.ini" );
48cfg.addKey("GRID2DPATH", "Output", "./");
49IOManager io( cfg );
50
51DEMObject dem;
52dem.setUpdatePpt(DEMObject::SLOPE);
53io.readDEM(dem);
54
55//read and initialize the POIs
56std::vector<StationData> vecStation;
57std::vector<Coords> pts;
58io.readPOI(pts);
59setPOI(dem, pts, vecStation);
60
61Grid2DObject lus;
62io.read2DGrid(lus, "arolla_rad.lus");
63
64Glaciers glacier(cfg, dem);
65glacier.setGlacierMap( computeGlacierMap(lus) );
66/*Grid2DObject alt, dist;
67glacier.getGrids(alt, dist);
68io.write2DGrid(computeGlacierMap(lus), "lus.asc");
69io.write2DGrid(alt, "alt.asc");
70io.write2DGrid(dist, "dist.asc");
71exit(0);*/
72
73const double TZ = cfg.get("TIME_ZONE", "Input");
74Date d1, d2;
75IOUtils::convertString(d1, argv[1], TZ);
76IOUtils::convertString(d2, argv[2], TZ);
77const double Tstep = 1./24.;
78
79//loop over the timesteps and fill the POIs
80std::vector< std::vector<MeteoData> > vecvecMeteo;
81vecvecMeteo.insert(vecvecMeteo.begin(), vecStation.size(), std::vector<MeteoData>()); //allocation for the vectors
82
83for(; d1<=d2; d1+=Tstep) { //time loop
84Grid2DObject ta;
85io.getMeteoData(d1, dem, MeteoData::TA, ta);
86Grid2DObject ta_corr(ta);
87glacier.correctTemperatures(ta_corr);
88
89for(unsigned int ii=0; ii<pts.size(); ii++) { //data extraction for 1 station at a time
90MeteoData meteo(d1, vecStation[ii]);
91meteo(MeteoData::TA) = ta.grid2D( pts[ii].getGridI(), pts[ii].getGridJ());
92meteo.addParameter("TA_corr");
93meteo("TA_corr") = ta_corr.grid2D( pts[ii].getGridI(), pts[ii].getGridJ());
94vecvecMeteo[ii].push_back( meteo );
95}
96}
97
98//write out the POIs' timeseries
99std::cout << "Writing output data" << std::endl;
100io.writeMeteoData(vecvecMeteo);
101}
102
103
104
105
106
107
108

Archive Download this file

Revision: HEAD