FFSM++  1.1.0
French Forest Sector Model ++
Layers.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2015 by Laboratoire d'Economie Forestière *
3  * http://ffsm-project.org *
4  * *
5  * This program is free software; you can redistribute it and/or modify *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 3 of the License, or *
8  * (at your option) any later version, given the compliance with the *
9  * exceptions listed in the file COPYING that is distribued together *
10  * with this file. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program; if not, write to the *
19  * Free Software Foundation, Inc., *
20  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
21  ***************************************************************************/
22 #include <QtCore>
23 
24 #include <math.h>
25 #include <algorithm>
26 
27 #include "Layers.h"
28 #include "Gis.h"
29 #include "ThreadManager.h"
30 #include "Scheduler.h"
31 
32 Layers::Layers(ThreadManager* MTHREAD_h, string name_h, string label_h, bool isInteger_h, bool dynamicContent_h, string fullFilename_h, bool display_h)
33 {
34  MTHREAD=MTHREAD_h;
35  name = name_h;
36  label = label_h;
37  isInteger = isInteger_h;
38  dynamicContent = dynamicContent_h;
39  fullFileName = fullFilename_h;
40  display = display_h;
41 }
42 
44 {
45 }
46 
47 void
48 Layers::addLegendItem(int ID_h, string label_h, int rColor_h, int gColor_h, int bColor_h, double minValue_h, double maxValue_h){
49 
50  for (uint i=0;i<legendItems.size();i++){
51  if (legendItems.at(i).ID == ID_h){
52  msgOut(MSG_ERROR, "Trying to add a legend item that already exist on this layer (layer: "+label+" - legend label: "+label_h+")");
53  //cout << "ID: "<<ID_h<<" Label: "<<label_h<<" minValue: "<<minValue_h << " maxValue: "<<maxValue_h<<endl;
54  return;
55  }
56  }
57 
58  LegendItems ITEM;
59  ITEM.ID = ID_h;
60  ITEM.label = label_h;
61  ITEM.rColor = rColor_h;
62  ITEM.gColor = gColor_h;
63  ITEM.bColor = bColor_h;
64  ITEM.minValue = minValue_h;
65  ITEM.maxValue = maxValue_h;
66  ITEM.cashedCount=0;
67  legendItems.push_back(ITEM);
68 
69 }
70 
71 void
72 Layers::addLegendItems(vector<LegendItems> legendItems_h){
73  vector <LegendItems> toAdd;
74  for(uint i=0; i<legendItems_h.size();i++){
75  bool existing = false;
76  for (uint j=0;j<legendItems.size();j++){
77  if(legendItems_h[i].ID == legendItems[j].ID){
78  existing = true;
79  break;
80  }
81  }
82  if(existing){
83  msgOut(MSG_WARNING, "Legend item "+i2s(legendItems_h[i].ID)+" non added on layer "+this->name+" as already existing.");
84  } else {
85  toAdd.push_back(legendItems_h[i]);
86  }
87  }
88  legendItems.insert( legendItems.end(), toAdd.begin(), toAdd.end() );
89 }
90 
91 
92 /**
93 Used in the init stage, this function take as input the real map code as just read from the map file, and filter it according to the reclassification rules.
94 @see ReclassRules
95 */
96 double
98  bool check =false;
99  std::vector <double> cumPVector;
100  std::vector <double> outCodesVector;
101  double cumP = 0;
102  double returnCode=0;
103 
104  for(uint i=0; i<reclassRulesVector.size(); i++){
105  if (reclassRulesVector.at(i).inCode == code_h){
106  check = true;
107  cumP += reclassRulesVector.at(i).p;
108  cumPVector.push_back(cumP);
109  outCodesVector.push_back(reclassRulesVector.at(i).outCode);
110  }
111  }
112  if (!check) {return code_h;}
113  if (cumP <= 0.99999999 || cumP >= 1.00000001){msgOut(MSG_CRITICAL_ERROR,"the sum of land use reclassification rules is not 1 for at least one input code (input code: "+d2s(code_h)+"; cumP: "+d2s(cumP)+")");}
114  double random;
115  //srand(time(NULL)); // this would re-initialise the random seed
116  random = ((double)rand() / ((double)(RAND_MAX)+(double)(1)) );
117  for(uint i=0; i<cumPVector.size(); i++){
118  if (random <= cumPVector.at(i)){
119  returnCode = outCodesVector.at(i);
120  break;
121  }
122  }
123  return returnCode;
124 }
125 
126 /**
127 This function take as input the value stored in the pixel for the specific layer, loops over the legend item and find the one that match it, returning its color.
128 <br>If the layer is of type integer, the match is agains legendItem IDs, otherwise we compare the legendItem ranges.
129 @see LegendItems
130 */
131 QColor
132 Layers::getColor(double ID_h){
133  QColor nocolor(255,255,255);
134  if (ID_h == MTHREAD->GIS->getNoValue()){
135  return nocolor;
136  }
137  if (isInteger){
138  for(uint i=0; i<legendItems.size(); i++){
139  if (legendItems.at(i).ID == ((int)ID_h)){
140  QColor color(legendItems.at(i).rColor, legendItems.at(i).gColor, legendItems.at(i).bColor);
141  return color;
142  }
143  }
144  return nocolor;
145  }
146  else {
147  for(uint i=0; i<legendItems.size(); i++){
148  if (ID_h < legendItems.at(i).maxValue && ID_h >= legendItems.at(i).minValue){
149  QColor color(legendItems.at(i).rColor, legendItems.at(i).gColor, legendItems.at(i).bColor);
150  return color;
151  }
152  }
153  return nocolor;
154  }
155 }
156 /**
157 This function take as input the value stored in the pixel for the specific layer, loops over the legend item and find the one that match it, returning its label.
158 <br>If the layer is of type integer, the match is agains legendItem IDs, otherwise we compare the legendItem ranges.
159 @see LegendItems
160 */
161 string
162 Layers::getCategory(double ID_h){
163  if (ID_h == MTHREAD->GIS->getNoValue()){
164  return "";
165  }
166  if (isInteger){
167  for(uint i=0; i<legendItems.size(); i++){
168  if (legendItems.at(i).ID == ((int)ID_h)){
169  return legendItems.at(i).label;
170  }
171  }
172  return "";
173  }
174  else {
175  for(uint i=0; i<legendItems.size(); i++){
176  if (ID_h < legendItems.at(i).maxValue && ID_h >= legendItems.at(i).minValue){
177  return legendItems.at(i).label;
178  }
179  }
180  return "";
181  }
182 }
183 
184 
185 
186 
187 void
189 
190  for (uint i=0; i<legendItems.size(); i++){
191  legendItems.at(i).cashedCount=0; //initialized with 0 values...
192  }
193  double totPixels = MTHREAD->GIS->getXyNPixels();
194  double pixelValue;
195  for (uint j=0;j<totPixels;j++){
196  pixelValue = MTHREAD->GIS->getPixel(j)->getDoubleValue(name);
197  if (isInteger){
198  for(uint i=0; i<legendItems.size(); i++){
199  if (legendItems.at(i).ID == ((int)pixelValue)){
200  legendItems.at(i).cashedCount++;
201  break;
202  }
203  }
204  }
205  else {
206  for(uint i=0; i<legendItems.size(); i++){
207  if (pixelValue < legendItems.at(i).maxValue && pixelValue >= legendItems.at(i).minValue){
208  legendItems.at(i).cashedCount++;
209  break;
210  }
211  }
212  }
213  }
214  if (debug){
215  msgOut(MSG_INFO, "Layer statistics - Count by Legend items");
216  msgOut(MSG_INFO, "Layer name: "+label);
217  msgOut(MSG_INFO, "Total plots: "+ d2s(totPixels));
218  for(uint i=0;i<legendItems.size();i++){
219  msgOut(MSG_INFO, legendItems.at(i).label+": "+i2s(legendItems.at(i).cashedCount));
220  }
221  }
222 }
223 void
225 
226 
227  vector <double> origValues;
228  int maskValue = -MTHREAD->GIS->getNoValue();
229  double totPixels = MTHREAD->GIS->getXyNPixels();
230  for (uint i=0;i<totPixels;i++){
231  double pxValue= MTHREAD->GIS->getPixel(i)->getDoubleValue(name);
232  if(pxValue != MTHREAD->GIS->getNoValue()){
233  origValues.push_back(pxValue);
234  MTHREAD->GIS->getPixel(i)->changeValue(name,maskValue);
235  }
236  }
237  random_shuffle(origValues.begin(), origValues.end()); // randomize the elements of the array.
238 
239  for (uint i=0;i<totPixels;i++){
240  double pxValue= MTHREAD->GIS->getPixel(i)->getDoubleValue(name);
241  if(pxValue != MTHREAD->GIS->getNoValue()){
242  double toChangeValue = origValues.at(origValues.size()-1);
243  //cout << toChangeValue << endl;
244  origValues.pop_back();
245  MTHREAD->GIS->getPixel(i)->changeValue(name,toChangeValue);
246  }
247  }
248 
249 }
250 void
252 
253  if(MTHREAD->MD->getIntSetting("outputLevel",DATA_NOW)<OUTVL_MAPS) return;
254  //if(!display || !dynamicContent) return; // already checked in Gis::printLayers() and allowed only for the first year
255  string mapBaseDirectory = MTHREAD->MD->getBaseDirectory()+MTHREAD->MD->getOutputDirectory()+"maps/";
256  string mapGridOutputDirectory = mapBaseDirectory+"asciiGrids/";
257  string catsOutputDirectory = mapBaseDirectory+"cats/";
258  string coloursOutputDirectory = mapBaseDirectory+"colr/";
259 
260  string mapFilename = mapGridOutputDirectory +name+ "_" +i2s(MTHREAD->SCD->getYear()) +"_" +MTHREAD->getScenarioName();
261  string catsFilename = catsOutputDirectory +name+ "_" +i2s(MTHREAD->SCD->getYear()) +"_" +MTHREAD->getScenarioName();
262  string coloursFilename = coloursOutputDirectory +name+ "_" +i2s(MTHREAD->SCD->getYear()) +"_" +MTHREAD->getScenarioName();
263  string filenameListIntLayers = mapBaseDirectory+"integerListLayers/"+MTHREAD->getScenarioName();
264  string filenameListFloatLayers = mapBaseDirectory+"floatListLayers/"+MTHREAD->getScenarioName();
265 
266  // printing the map...
267  string header;
268  if(MTHREAD->MD->getIntSetting("mapOutputFormat",DATA_NOW) == 1){ // GRASS ASCII Grid
269  header = "north: " + d2s(MTHREAD->GIS->getGeoTopY()) + "\n"
270  + "south: " + d2s(MTHREAD->GIS->getGeoBottomY()) + "\n"
271  + "east: " + d2s(MTHREAD->GIS->getGeoRightX()) + "\n"
272  + "west: " + d2s(MTHREAD->GIS->getGeoLeftX()) + "\n"
273  + "rows: " + i2s(MTHREAD->GIS->getYNPixels()) + "\n"
274  + "cols: " + i2s(MTHREAD->GIS->getXNPixels()) + "\n"
275  + "null: " + d2s(MTHREAD->GIS->getNoValue()) + "\n";
276 
277  } else if(MTHREAD->MD->getIntSetting("mapOutputFormat",DATA_NOW) == 2){
278  header = "ncols: " + i2s(MTHREAD->GIS->getXNPixels()) + "\n"
279  + "lrows: " + i2s(MTHREAD->GIS->getYNPixels()) + "\n"
280  + "xllcornel: " + d2s(MTHREAD->GIS->getGeoLeftX()) + "\n"
281  + "yllcorner: " + d2s(MTHREAD->GIS->getGeoBottomY()) + "\n"
282  + "cellsize: " + d2s(MTHREAD->GIS->getXMetersByPixel()) + "\n"
283  + "nodata_value: " + d2s(MTHREAD->GIS->getNoValue()) + "\n";
285  msgOut(MSG_ERROR, "The X resolution is different to the Y resolution. I am exporting the map in ArcInfo ASCII Grid format using the X resolution, but be aware that it is incorrect, as this format doesn't support different X-Y resolutions.");
286  }
287 
288  } else {
289  msgOut(MSG_ERROR,"Map not print for unknow output type.");
290  }
291 
292  ofstream outm; //out map
293  outm.open(mapFilename.c_str(), ios::out); //ios::app to append..
294  if (!outm){ msgOut(MSG_ERROR,"Error in opening the file "+mapFilename+".");}
295  outm << header << "\n";
296 
297  for (int i=0;i<MTHREAD->GIS->getYNPixels();i++){
298  for (int j=0;j<MTHREAD->GIS->getXNPixels();j++){
299  outm << MTHREAD->GIS->getPixel(j, i)->getDoubleValue(name) << " ";
300  }
301  outm << "\n";
302  }
303  outm.close();
304 
305  //printing the cat file
306  ofstream outc; //out category file
307  outc.open(catsFilename.c_str(), ios::out); //ios::app to append..
308  if (!outc){ msgOut(MSG_ERROR,"Error in opening the file "+catsFilename+".");}
309  outc << "# " << name << "_-_" << i2s(MTHREAD->SCD->getYear()) << "\n\n\n";
310  outc << "0.00 0.00 0.00 0.00"<<"\n";
311 
312  if (isInteger){
313  for(uint i=0;i<legendItems.size();i++){
314  outc << legendItems[i].ID << ":"<< legendItems[i].label << "\n";
315  }
316  }
317  else {
318  for(uint i=0;i<legendItems.size();i++){
319  outc << legendItems[i].minValue << ":"<< legendItems[i].maxValue << ":"<< legendItems[i].label << "\n";
320  }
321  }
322 
323  //printing the colour legend file
324  ofstream outcl; //out colour file
325  outcl.open(coloursFilename.c_str(), ios::out); //ios::app to append..
326  if (!outcl){ msgOut(MSG_ERROR,"Error in opening the file "+coloursFilename+".");}
327  outcl << "% " << name << "_-_" << i2s(MTHREAD->SCD->getYear()) << "\n\n\n";
328 
329  if (isInteger){
330  for(uint i=0;i<legendItems.size();i++){
331  outcl << legendItems[i].ID << ":"<< legendItems[i].rColor << ":" << legendItems[i].gColor << ":" << legendItems[i].bColor << "\n";
332  }
333  }
334  else {
335  for(uint i=0;i<legendItems.size();i++){
336  outcl << legendItems[i].minValue << ":"<< legendItems[i].rColor << ":" << legendItems[i].gColor << ":" << legendItems[i].bColor << " "<< legendItems[i].maxValue << ":"<< legendItems[i].rColor << ":" << legendItems[i].gColor << ":" << legendItems[i].bColor << "\n";
337  }
338  }
339 
340  // adding the layer to the list of saved layers..
341  ofstream outList;
342  if (isInteger){
343  outList.open(filenameListIntLayers.c_str(), ios::app); // append !!!
344  outList << name << "_" << MTHREAD->SCD->getYear() << "_" << MTHREAD->getScenarioName() << "\n";
345  }
346  else {
347  outList.open(filenameListFloatLayers.c_str(), ios::app); // append !!!
348  outList << name << "_" << MTHREAD->SCD->getYear() << "_" << MTHREAD->getScenarioName() << "\n";
349  }
350  outList.close();
351 }
352 
353 void
355 
356  if(!display || !dynamicContent) return;
357 
358  int xNPixels = MTHREAD->GIS->getXNPixels();
359  int subXR = MTHREAD->GIS->getSubXR();
360  int subXL = MTHREAD->GIS->getSubXL();
361  int subYT = MTHREAD->GIS->getSubYT();
362  int subYB = MTHREAD->GIS->getSubYB();
363 
364  string mapBaseDirectory = MTHREAD->MD->getBaseDirectory()+MTHREAD->MD->getOutputDirectory()+"maps/bitmaps/";
365  string mapFilename = mapBaseDirectory +name+ "_" +i2s(MTHREAD->SCD->getYear()) +"_" +MTHREAD->getScenarioName()+".png";
366 
367  QImage image = QImage(subXR-subXL+1, subYB-subYT+1, QImage::Format_RGB32);
368  image.fill(qRgb(255, 255, 255));
369  for (int countRow=subYT;countRow<subYB;countRow++){
370  for (int countColumn=subXL;countColumn<subXR;countColumn++){
371  double value = MTHREAD->GIS->getPixel(countRow*xNPixels+countColumn)->getDoubleValue(name);
372  QColor color = this->getColor(value);
373  image.setPixel(countColumn-subXL,countRow-subYT,color.rgb());
374  }
375  }
376  image.save(mapFilename.c_str());
377 }
void addLegendItems(vector< LegendItems > legendItems_h)
Definition: Layers.cpp:72
double getYMetersByPixel() const
Definition: Gis.h:141
Print an ERROR message, but don&#39;t stop the model.
Definition: BaseClass.h:61
The required data is for the current year.
Definition: BaseClass.h:73
double getGeoLeftX() const
Definition: Gis.h:138
int getIntSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1105
double getGeoBottomY() const
Definition: Gis.h:137
void print()
Print the layer content as an ASCII grid map with its companion files (classification and colors)...
Definition: Layers.cpp:251
string i2s(const int &int_h) const
integer to string conversion
Definition: BaseClass.cpp:219
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
void changeValue(const string &layerName_h, const double &value_h, const bool &setNoValueForZero=false)
Change the value of an existing layerMTHREAD->GIS->pack(parName, forName, dClass, year)...
Definition: Pixel.cpp:135
string label
Label of the layer (spaces allowed)
Definition: Layers.h:99
ModelData * MD
the model data object
Definition: ThreadManager.h:72
string getBaseDirectory() const
Definition: ModelData.h:118
string name
ID of the layer (no spaces allowed)
Definition: Layers.h:98
vector< LegendItems > legendItems
Vector of legend items.
Definition: Layers.h:104
double getNoValue() const
Definition: Gis.h:133
double maxValue
Definition: Layers.h:122
Output verbosity level print (also) the maps in ascii grid format.
Definition: BaseClass.h:87
bool isInteger
Type of the layer (true==integer, false==double. If true, on each legend item: minValue==maxValue==ID...
Definition: Layers.h:100
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
Definition: ThreadManager.h:75
void msgOut(const int &msgCode_h, const string &msg_h, const bool &refreshGUI_h=true) const
Overloaded function to print the output log.
Definition: BaseClass.cpp:50
QColor getColor(double ID_h)
Evaluates all the legend items to find the one that match the input code, and return its color as a Q...
Definition: Layers.cpp:132
double getXyNPixels() const
Return the number of pixels on Y.
Definition: Gis.h:131
int getSubYT() const
Definition: Gis.h:144
Gis * GIS
GIS information and methods.
Definition: ThreadManager.h:73
Thread manager. Responsable to manage the main thread and "speak" with the GUI.
Definition: ThreadManager.h:65
void randomShuffle()
For some sensitivity analisys, random the values for this layer for not-empty values (only integer la...
Definition: Layers.cpp:224
double getXMetersByPixel() const
Definition: Gis.h:140
string fullFileName
Filename of the associated dataset (map)
Definition: Layers.h:103
void addLegendItem(int ID_h, string label_h, int rColor_h, int gColor_h, int bColor_h, double minValue_h, double maxValue_h)
Add a legend item.
Definition: Layers.cpp:48
int getYNPixels() const
Return the number of pixels on X.
Definition: Gis.h:130
Layers(ThreadManager *MTHREAD_h, string name_h, string label_h, bool isInteger_h, bool dynamicContent_h, string fullFilename_h, bool display_h=true)
In the constructor we set the main layer properties.
Definition: Layers.cpp:32
Print an error message and stop the model.
Definition: BaseClass.h:62
double getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue=false) const
Return the value for a specific layer.
Definition: Pixel.cpp:158
string d2s(const double &double_h) const
double to string conversion
Definition: BaseClass.cpp:229
int getYear()
Definition: Scheduler.h:49
double filterExogenousDataset(double code_h)
Used to reclassify the land use map for "generic" categories.
Definition: Layers.cpp:97
Legend items.
Definition: Layers.h:115
int bColor
Definition: Layers.h:120
int ID
Definition: Layers.h:116
void printBinMap()
Print a binary reppresentation of the data (a standard image, e.g. a .png file). It prints only the s...
Definition: Layers.cpp:354
double minValue
Definition: Layers.h:121
int getSubYB() const
Definition: Gis.h:145
Print a WARNING message.
Definition: BaseClass.h:60
int gColor
Definition: Layers.h:119
int getSubXR() const
Definition: Gis.h:143
string getScenarioName()
void countMyPixels(bool debug=false)
Count the pixels going to each legend item and print them if debug==true.
Definition: Layers.cpp:188
string getCategory(double ID_h)
Evaluates all the legend items to find the one that match the input code, and return its label...
Definition: Layers.cpp:162
int getSubXL() const
Definition: Gis.h:142
Pixel * getPixel(int x_h, int y_h)
Definition: Gis.h:134
int cashedCount
count the pixels whitin a item range
Definition: Layers.h:123
string label
Definition: Layers.h:117
string getOutputDirectory() const
Return a vector of objects that together provide the specified resource in the specified quantity...
Definition: ModelData.h:113
int getXNPixels() const
Definition: Gis.h:129
vector< ReclassRules > reclassRulesVector
Vector of initial reclassification rules.
Definition: Layers.h:105
double getGeoTopY() const
Return a pixel pointer from its ID.
Definition: Gis.h:136
Print an INFO message.
Definition: BaseClass.h:59
double getGeoRightX() const
Definition: Gis.h:139
~Layers()
Definition: Layers.cpp:43
bool display
Normally true, but some layers used to just keep data shoudn&#39;t be normally processed.
Definition: Layers.h:102
int rColor
Definition: Layers.h:118
bool dynamicContent
True if the content may change during simulation year.
Definition: Layers.h:101