FFSM++  1.1.0
French Forest Sector Model ++
Pixel.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 "Pixel.h"
23 #include "ThreadManager.h"
24 #include "Scheduler.h"
25 #include "Init.h"
26 
27 Pixel::Pixel(double ID_h, ThreadManager* MTHREAD_h): ID(ID_h)
28 {
29  MTHREAD=MTHREAD_h;
30  int nft = MTHREAD->MD->getForTypeIds().size();
31  vector<double> temp(nft,1);
32  //vector<double> temp2(nft,0);
33  spMods = temp;
34  avalCoef = temp;
35  //vMort = temp2;
36  //std::fill(v.begin(), v.end(), 0);
37 }
38 
40 {
41 }
42 
43 /**
44 The function return a vector of pointers to Pixels at the gived distance from the caller pixel.\\
45 The list start with those on the Top, then add those on the right, those on the bottom and those on the left. Finally it had the corner pixels (that are more far).\\
46 It takes into consideration borders correctly.
47 
48 Fully tested on internal points as well semi-border cases, border cases and corner cases. ALL OK.
49 
50 @param distLevel_h Distance in number of adiacent pixels. It has to be at least 1 (the function return an error if it is 0).
51 */
52 vector <Pixel *>
53 Pixel::getPixelsAtDistLevel(int distLevel_h) const{
54 
55  if (distLevel_h<1) {
56  msgOut(MSG_CRITICAL_ERROR, "getPixelsAtDistLevel() is defined for distances of at least 1 !");
57  }
58 
59  vector <Pixel *> toReturn;
60  int xNPixels = MTHREAD->GIS->getXNPixels();
61  int yNPixels = MTHREAD->GIS->getYNPixels();
62  int thisX = this->getX();
63  int thisY = this->getY();
64  int minX = max(0 , (thisX - distLevel_h)+1);
65  int maxX = min(xNPixels , thisX + distLevel_h);
66  int minY = max(0 , (thisY - distLevel_h)+1);
67  int maxY = min(yNPixels , thisY + distLevel_h);
68 
69  // getting the top pixels (corner exluded)...
70  if (thisY-distLevel_h >=0){
71  for(int i=minX;i<maxX;i++){
72  toReturn.push_back(MTHREAD->GIS->getPixel(i,thisY-distLevel_h));
73  }
74  }
75  // getting the right pixels (corner exluded)...
76  if (thisX+distLevel_h < xNPixels){
77  for(int i=minY;i<maxY;i++){
78  toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,i));
79  }
80  }
81  // getting the bottom pixels (corner exluded)...
82  if (thisY+distLevel_h < yNPixels){
83  for(int i=minX;i<maxX;i++){
84  toReturn.push_back(MTHREAD->GIS->getPixel(i,thisY+distLevel_h));
85  }
86  }
87  // getting the left pixels (corner exluded)...
88  if (thisX-distLevel_h >= 0){
89  for(int i=minY;i<maxY;i++){
90  toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,i));
91  }
92  }
93 
94  // getting the corners (correctly at the end, after already retrieved the other pixels..)...
95  // top-left..
96  if (thisX-distLevel_h >= 0 && thisY-distLevel_h >=0){
97  toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,thisY-distLevel_h));
98  }
99  // top-right..
100  if (thisX+distLevel_h < xNPixels && thisY-distLevel_h >=0){
101  toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,thisY-distLevel_h));
102  }
103  // bottom-right..
104  if (thisX+distLevel_h < xNPixels && thisY+distLevel_h <yNPixels){ // bug discovered 20070719
105  toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,thisY+distLevel_h));
106  }
107  // bottom-left..
108  if (thisX-distLevel_h >= 0 && thisY+distLevel_h <yNPixels){
109  toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,thisY+distLevel_h));
110  }
111  return toReturn;
112 }
113 
114 /*void //moved as inline function
115 Pixel::setValue(const string & layerName_h, const double & value_h){
116 
117  //tempValuePair.first = layerName_h; // type of first is string
118  //tempValuePair.second = value_h; // type of second is double
119  //tempValuePair = make_pair (layerName_h,value_h);
120  values.insert(pair<string, double>(layerName_h, value_h));
121  //values.insert(tempValuePair);
122 
123 
124 }*/
125 
126 /*
127 inline void
128 Pixel::setValue (const string& parName, const string& forName, const string& dClass, const int& year, const double& value_h){
129  values.insert(pair<string, double>(MTHREAD->GIS->pack(parName, forName, dClass, year), value_h));
130 }
131 */
132 
133 
134 void
135 Pixel::changeValue(const string &layerName_h, const double &value_h, const bool &setNoValueForZero){
136  map<string, double>::iterator p;
137  p=values.find(layerName_h);
138  if(p != values.end()){
139  if(setNoValueForZero && value_h == 0){
140  p->second = MTHREAD->GIS->getNoValue();
141  } else {
142  p->second = value_h;
143  }
144  } else {
145  msgOut(MSG_ERROR, "Coud not change pixel value for layer "+layerName_h+". Layer don't found.");
146  }
147  return;
148 }
149 
150 /*
151 void
152 Pixel::changeValue (const double &value_h, const string& parName, const string& forName, const string &dClass, const int &year, const bool &setNoValueForZero){
153  changeValue(MTHREAD->GIS->pack(parName, forName, dClass, year), value_h, setNoValueForZero);
154 }
155 */
156 
157 double
158 Pixel::getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue) const{
159  vIter=values.find(layerName_h);
160  if(vIter != values.end()) {
161  if(returnZeroForNoValue && vIter->second==MTHREAD->GIS->getNoValue()){
162  return 0.0;
163  } else {
164  return vIter->second;
165  }
166  } else {
167  msgOut(MSG_WARNING, "No layer \""+layerName_h+"\" found on pixel ("+i2s(getX())+","+i2s(getY())+"). Sure you didn't mispelled it?");
168  if(returnZeroForNoValue){
169  return 0.0;
170  } else {
171  return MTHREAD->GIS->getNoValue();
172  }
173  }
174 }
175 
176 /**
177 getMultiplier() returns the value of the multiplier as memorized in the spatialDataSubfolder layers or in the forData table.
178 It will looks for exact match or for lower years if available.
179 If no layers are available or the usePixelData option is not used, it will return 1.
180 If the tp_multiplier is asked for, it will adjusts for spatial variance coefficient.
181 If the mortCoef_multiplier is used and we are in the table settings it will adjust it by mortCoef_link.
182 */
183 double
184 Pixel::getMultiplier(const string& multiplierName, const string& forName, int year){
185 
186 
187  if(year==DATA_NOW){year = MTHREAD->SCD->getYear();}
188 
189 
190  double multiplierSpVar = (multiplierName == "tp_multiplier")?getSpModifier(forName):1.0;
191 
192  vector <string> modifiersFromTable = MTHREAD->MD->getStringVectorSetting("modifiersFromTable");
193 
194  if(std::find(modifiersFromTable.begin(), modifiersFromTable.end(), multiplierName) != modifiersFromTable.end()) {
195  // load multiplier from forData table..
196  int regId = getMyRegion()->getRegId();
197  double multiplier = MTHREAD->MD->getForData(multiplierName, regId, forName, "",year);
198  if (multiplierName == "mortCoef_multiplier"){
199  return pow(multiplier,MTHREAD->MD->getDoubleSetting("mortMultiplier_link",DATA_NOW))*multiplierSpVar; //Added to account that our multipliers are based on probability of presence and not on planted/managed forests, where mortality is somhow reduced
200  }
201  return multiplier*multiplierSpVar;
202 
203  } else {
204  // load multiplier from layer file..
205  return getSTData(multiplierName, forName, year, "", 1.0);
206 
207 // // return 1 if not using pixel mode
208 // if(!MTHREAD->MD->getBoolSetting("usePixelData")) return 1.0;
209 // string search_for = multiplierName+"#"+forName+"##"+i2s(year);
210 // map <string,double>::const_iterator i = values.upper_bound(search_for); //return the position always upper to the found one, even if it's an equal match.
211 // if(i!= values.begin()) i--; // this rewind the position to the one just before or equal
212 // const string& key = i->first;
213 // string search_base = search_for.substr(0,search_for.size()-4);
214 // if (key.compare(0, search_base.size(), search_base) == 0){
215 // //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
216 // //if(i->second != 1){
217 // // cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
218 // // exit(0);
219 // //}
220 // return i->second*multiplierSpVar;
221 // } else {
222 // //cout << "NOTM: " << search_for <<", "<< i->first << endl;
223 // return 1.0*multiplierSpVar;
224 // }
225 
226  }
227 }
228 
229 /**
230 getSTData() (standing for "get spatial temporal data") returns the value memorized in a spatial layer whose name follows the above convention
231 
232 parName#fType#dimension2#year
233 
234 It queries the layer framework with the following logic:
235 
236 @param parName Parameter name (exact match)
237 @param forName Forest type (exact match or look for "")
238 @param year Year (exact match or look for the highest lower value). Accept enum DATA_NOW
239 @param d2 Optional dimension (exact match, string. Default = "")
240 @param naToReturn Optional value to return if no match (default to secnario-specific novalue)
241 
242 
243 */
244 
245 double
246 Pixel::getSTData(const string& parName, const string& forName, int year, const string& d2, double naToReturn){
247 
248  if (naToReturn == RETNA) naToReturn = MTHREAD->MD->getDoubleSetting("noValue");
249 
250  // return NA if not using pixel mode
251  if(!MTHREAD->MD->getBoolSetting("usePixelData")) return naToReturn;
252 
253  // If year is DATA_NOW get current year
254  if (year == DATA_NOW) year = MTHREAD->SCD->getYear();
255 
256  // Looking for an exact match in forName
257  string search_for = parName+"#"+forName+"#"+d2+"#"+i2s(year);
258  map <string,double>::const_iterator i = values.upper_bound(search_for); //return the position always upper to the found one, even if it's an equal match.
259  if(i!= values.begin()) i--; // this rewind the position to the one just before or equal
260  const string& key = i->first;
261  string search_base = search_for.substr(0,search_for.size()-4);
262  if (key.compare(0, search_base.size(), search_base) == 0){
263  //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
264  //if(i->second != 1){
265  // cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
266  // exit(0);
267  //}
268  return i->second;
269  } else {
270  // An exact match in forName has not being found, look for no forName:
271  string search_for2 = parName+"#"+"#"+d2+"#"+i2s(year);
272  map <string,double>::const_iterator i2 = values.upper_bound(search_for2); //return the position always upper to the found one, even if it's an equal match.
273  if(i2!= values.begin()) i2--; // this rewind the position to the one just before or equal
274  const string& key2 = i2->first;
275  string search_base2 = search_for2.substr(0,search_for2.size()-4);
276  if (key2.compare(0, search_base2.size(), search_base2) == 0){
277  //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
278  //if(i->second != 1){
279  // cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
280  // exit(0);
281  //}
282  return i2->second;
283  } else {
284  //cout << "NOTM: " << search_for <<", "<< i->first << endl;
285  return naToReturn;
286  }
287  }
288 }
289 
290 
291 /**
292 The mortality returned is the increased yearly mortality due to any affecting pathogenes.
293 The function load the relevant pathogen mortality rule(s), for each of them check for how many years the phatogen is present with concentrations
294 above the threshold and returns the relavant increase in mortality (summing them in case of multiple pathogens).
295 
296 */
297 double
298 Pixel::getPathMortality(const string& forType, const string& dC, int year){
299  if(!MTHREAD->MD->getBoolSetting("usePathogenModule")) return 0.0;
300 
301  string debug=forType;
302  int initialOptYear = MTHREAD->MD->getIntSetting("initialOptYear");
303  int simulationYears = MTHREAD->MD->getIntSetting("simulationYears");
304 
305  int maxYear = initialOptYear + simulationYears;
306 
307  vector<pathRule*> pathRules = MTHREAD->MD->getPathMortalityRule(forType,dC);
308 
309  double pathMort = 0.0;
310  if(year==DATA_NOW){year = MTHREAD->SCD->getYear();}
311 
312  for(uint r=0;r<pathRules.size();r++){
313  string pathId=pathRules[r]->pathId;
314  double pres_min=pathRules[r]->pres_min;
315  vector<double> mortCoefficients=pathRules[r]->mortCoefficents;
316  double pathMort_thispath = 0.0;
317  for(uint y=year;y>(year-mortCoefficients.size());y--){
318  int i =year-y;
319  int y2 = y;
320  if(y>=maxYear){
321  y2=maxYear-1;
322  }
323 
324  string layerName="pathogen_pp#"+pathId+"#"+i2s(y2);
325  if(MTHREAD->GIS->layerExist(layerName)){
326  if (this->getDoubleValue(layerName,true)>= pres_min){
327  pathMort_thispath = mortCoefficients[i];
328  }
329  }
330  }
331  pathMort += pathMort_thispath;
332  }
333  return pathMort;
334 
335 }
336 
337 void
338 Pixel::correctInputMultiplier (const string& multiplierName, const string& forName, double coefficient){
339  string search_for = multiplierName+"#"+forName+"#";
340  for (std::map<string,double>::iterator it=values.lower_bound(search_for); it!=values.end(); ++it){
341  if (it->first.compare(0, search_for.size(), search_for) == 0){
342  //cout << ID << ";" << forName << ";" << coefficient << endl;
343  it->second = it->second * coefficient;
344  }
345  }
346 }
347 
348 double
349 Pixel::getDoubleValue (const string& parName, const string& forName, const string& dClass, const int& year, const bool& returnZeroForNoValue){
350  return getDoubleValue(MTHREAD->GIS->pack(parName, forName, dClass, year), returnZeroForNoValue);
351 }
352 
353 void
355 
356 }
357 
358 double
359 Pixel::getPastRegArea(const int& ft_idx, const int& year){
360  map <int,vector<double> >::const_iterator i=regArea.find(year);
361  if(i != regArea.end()) {
362  return i->second.at(ft_idx);
363  } else {
364  msgOut(MSG_ERROR, "Asking for a pastRegArea of a not-registered year. I don't have year "+i2s(year)+"!");
365  }
366 }
367 
368 void
369 Pixel::setPastRegArea(const double& value, const int& ft_idx, const int& year){
370  msgOut(MSG_CRITICAL_ERROR,"TODO");
371  /*map <int,vector<double> >::const_iterator i=regArea.find(year);
372  if(i != regArea.end()) {
373  // we already have this year, let's see if the vector is bif enough
374  int currsize = i->second.size();
375  for(j=0;j<ft_idx-currside;j++){
376 
377  }
378  return i->second.at(ft_idx);
379  } else {
380  // new year
381  }
382 
383 
384  pair<int,vector<double> newRegArea;
385  */
386 
387 
388 }
389 
390 void
391 Pixel::swap(const int& swap_what){
392  switch (swap_what){
393  case VAR_VOL:
394  vol_l = vol;
395  break;
396  case VAR_AREA:
397  area_l = area;
398  break;
399  default:
400  msgOut(MSG_CRITICAL_ERROR,"Don't know how to swap "+swap_what);
401  }
402 }
403 
404 
405 double
406 Pixel::getSpModifier(const string& ft){
407  vector<string>ftypes = MTHREAD->MD->getForTypeIds();
408  for (int i=0;i<ftypes.size();i++){
409  if (ftypes[i] == ft){
410  return spMods.at(i);
411  }
412  }
413  msgOut(MSG_CRITICAL_ERROR,"Asked spatial modifier for a forest type that doesn't exist");
414 
415 }
416 
418 Pixel::getMyRegion(const int& rLevel){
419  if(rLevel==2){
420  return l2region;
421  } else if (rLevel==1) {
422  return l2region->getParent();
423  } else {
424  msgOut(MSG_ERROR, "Requested a unknown level region code in getMyRegion().");
425  }
426 }
string pack(const string &parName, const string &forName, const string &dClass, const int &year) const
Definition: Gis.h:148
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 getSpModifier(const string &ft)
Definition: Pixel.cpp:406
int getIntSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1105
vector< vector< double > > area
Definition: Pixel.h:107
bool getBoolSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1117
void setPastRegArea(const double &value, const int &ft_idx, const int &year)
Definition: Pixel.cpp:369
int getRegId() const
Definition: ModelRegion.h:66
Forest types (struct)
Definition: ModelData.h:293
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
bool layerExist(const string &layerName_h, bool exactMatch=true) const
Return a pointer to a layer given its name.
Definition: Gis.cpp:552
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
ModelData * MD
the model data object
Definition: ThreadManager.h:72
int getX() const
Definition: Pixel.h:68
map< string, double >::const_iterator vIter
Definition: Pixel.h:148
void correctInputMultiplier(const string &multiplierName, const string &forName, double coefficient=1)
It apply a given coefficient to all the multipliers layers of a given ft.
Definition: Pixel.cpp:338
vector< double > avalCoef
Availability (of wood resources) coefficient. A [0,1] coefficient (new: by forest type) that reduces ...
Definition: Pixel.h:144
double getNoValue() const
Definition: Gis.h:133
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
Definition: ThreadManager.h:75
ModelRegion * l2region
Pointer to level 2 region where this pixel is.
Definition: Pixel.h:155
vector< vector< double > > vol_l
store the volumes of the previous year
Definition: Pixel.h:126
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
vector< Pixel * > getPixelsAtDistLevel(int distLevel_h) const
Return a vector of pixels at the specified distance (in levels, not in physical units) ...
Definition: Pixel.cpp:53
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
vector< string > getStringVectorSetting(const string &name_h, int reg=WORLD) const
Definition: ModelData.cpp:1129
ModelRegion * getMyRegion(const int &rLevel=2)
Definition: Pixel.cpp:418
void newYear()
Definition: Pixel.cpp:354
~Pixel()
Definition: Pixel.cpp:39
vector< pathRule * > getPathMortalityRule(const string &forType, const string &dC)
Return the pathogen mortality rule(s) associated with a given ft and dc (plural as more than a single...
Definition: ModelData.cpp:2034
int getYNPixels() const
Return the number of pixels on X.
Definition: Gis.h:130
Print an error message and stop the model.
Definition: BaseClass.h:62
vector< double > spMods
The sampled spatial modifiers (by forest type)
Definition: Pixel.h:154
double getSTData(const string &parName, const string &forName, int year, const string &d2="", double naToReturn=RETNA)
Definition: Pixel.cpp:246
double getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue=false) const
Return the value for a specific layer.
Definition: Pixel.cpp:158
const double getForData(const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW)
Definition: ModelData.cpp:1281
Pixel(double ID_h, ThreadManager *MTHREAD_h)
Definition: Pixel.cpp:27
int getYear()
Definition: Scheduler.h:49
double getPastRegArea(const int &ft_idx, const int &year)
Definition: Pixel.cpp:359
vector< vector< double > > vol
Definition: Pixel.h:90
double getMultiplier(const string &multiplierName, const string &forName, int year=DATA_NOW)
Definition: Pixel.cpp:184
Request the (scenario specific) NO VALUE to be returned.
Definition: BaseClass.h:79
vector< vector< double > > area_l
store the areas of the previous year
Definition: Pixel.h:127
Print a WARNING message.
Definition: BaseClass.h:60
vector< string > getForTypeIds(bool all=false)
By default it doesn&#39;t return forTypes used only as input.
Definition: ModelData.cpp:430
map< string, double > values
Map of values for each layer.
Definition: Pixel.h:147
double getPathMortality(const string &forType, const string &dC, int year=DATA_NOW)
Return the INCREASED mortality due to pathogen presence for a given ft and dc in a certain year (defa...
Definition: Pixel.cpp:298
void swap(const int &swap_what)
Assign to the delayed value the current values, e.g. vol_l = vol.
Definition: Pixel.cpp:391
Pixel * getPixel(int x_h, int y_h)
Definition: Gis.h:134
double getDoubleSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1109
int getXNPixels() const
Definition: Gis.h:129
ModelRegion * getParent()
Definition: ModelRegion.h:72
int getY() const
Definition: Pixel.h:69
map< int, vector< double > > regArea
Definition: Pixel.h:114