Rev 968 | Rev 1042 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
671 | werner | 2 | /******************************************************************************************** |
3 | ** iLand - an individual based forest landscape and disturbance model |
||
4 | ** http://iland.boku.ac.at |
||
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl |
||
6 | ** |
||
7 | ** This program is free software: you can redistribute it and/or modify |
||
8 | ** it under the terms of the GNU General Public License as published by |
||
9 | ** the Free Software Foundation, either version 3 of the License, or |
||
10 | ** (at your option) any later version. |
||
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, see <http://www.gnu.org/licenses/>. |
||
19 | ********************************************************************************************/ |
||
20 | |||
106 | Werner | 21 | #include "global.h" |
22 | #include "standloader.h" |
||
23 | |||
290 | werner | 24 | |
106 | Werner | 25 | #include "grid.h" |
26 | #include "model.h" |
||
189 | iland | 27 | #include "resourceunit.h" |
106 | Werner | 28 | #include "speciesset.h" |
29 | |||
30 | #include "helper.h" |
||
290 | werner | 31 | #include "random.h" |
137 | Werner | 32 | #include "expression.h" |
140 | Werner | 33 | #include "expressionwrapper.h" |
284 | werner | 34 | #include "environment.h" |
287 | werner | 35 | #include "csvfile.h" |
549 | werner | 36 | #include "mapgrid.h" |
699 | werner | 37 | #include "snapshot.h" |
106 | Werner | 38 | |
697 | werner | 39 | /** @class StandLoader |
40 | @ingroup tools |
||
41 | loads (initializes) trees for a "stand" from various sources. |
||
42 | StandLoader initializes trees on the landscape. It reads (usually) from text files, creates the |
||
43 | trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid). |
||
44 | |||
45 | See http://iland.boku.ac.at/initialize+trees |
||
46 | */ |
||
381 | werner | 47 | // provide a mapping between "Picus"-style and "iLand"-style species Ids |
384 | werner | 48 | QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17; |
155 | werner | 49 | QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy"; |
732 | werner | 50 | |
290 | werner | 51 | StandLoader::~StandLoader() |
52 | { |
||
53 | if (mRandom) |
||
54 | delete mRandom; |
||
732 | werner | 55 | if (mHeightGridResponse) |
56 | delete mHeightGridResponse; |
||
290 | werner | 57 | } |
106 | Werner | 58 | |
59 | |||
281 | werner | 60 | void StandLoader::copyTrees() |
61 | { |
||
62 | // we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards |
||
63 | const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid(); |
||
64 | ResourceUnit **p = ruGrid.begin(); |
||
65 | if (!p) |
||
66 | throw IException("Standloader: invalid resource unit pointer!"); |
||
67 | ++p; // skip the first... |
||
68 | const QVector<Tree> &tocopy = mModel->ru()->trees(); |
||
69 | for (; p!=ruGrid.end(); ++p) { |
||
70 | QRectF rect = (*p)->boundingBox(); |
||
71 | foreach(const Tree& tree, tocopy) { |
||
72 | Tree &newtree = (*p)->newTree(); |
||
73 | newtree = tree; // copy tree data... |
||
739 | werner | 74 | newtree.setPosition(tree.position()+rect.topLeft()); |
281 | werner | 75 | newtree.setRU(*p); |
76 | newtree.setNewId(); |
||
77 | } |
||
78 | } |
||
431 | werner | 79 | if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied."; |
281 | werner | 80 | } |
81 | |||
284 | werner | 82 | /** main routine of the stand setup. |
83 | */ |
||
106 | Werner | 84 | void StandLoader::processInit() |
85 | { |
||
86 | GlobalSettings *g = GlobalSettings::instance(); |
||
192 | werner | 87 | XmlHelper xml(g->settings().node("model.initialization")); |
191 | werner | 88 | |
281 | werner | 89 | QString copy_mode = xml.value("mode", "copy"); |
90 | QString type = xml.value("type", ""); |
||
192 | werner | 91 | QString fileName = xml.value("file", ""); |
106 | Werner | 92 | |
732 | werner | 93 | bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false); |
736 | werner | 94 | mHeightGridTries = xml.valueDouble("heightGrid.maxTries", 10.); |
732 | werner | 95 | QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function |
96 | if (height_grid_enabled) { |
||
97 | QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init"); |
||
98 | qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file; |
||
99 | |||
100 | QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false)); |
||
743 | werner | 101 | if (!p->isValid()) { |
102 | throw IException(QString("Error when loading grid with tree heights for stand initalization: file %1 not found or not valid.").arg(init_height_grid_file)); |
||
103 | } |
||
732 | werner | 104 | height_grid.swap(p); |
105 | mInitHeightGrid = height_grid.data(); |
||
106 | |||
107 | QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)"); |
||
108 | mHeightGridResponse = new Expression(expr); |
||
109 | mHeightGridResponse->setLinearizationEnabled(true); |
||
110 | } |
||
111 | |||
106 | Werner | 112 | Tree::resetStatistics(); |
281 | werner | 113 | |
114 | // one global init-file for the whole area: |
||
284 | werner | 115 | if (copy_mode=="single") { |
281 | werner | 116 | loadInitFile(fileName, type); |
299 | werner | 117 | evaluateDebugTrees(); |
281 | werner | 118 | return; |
284 | werner | 119 | } |
281 | werner | 120 | |
137 | Werner | 121 | |
281 | werner | 122 | // call a single tree init for each resource unit |
123 | if (copy_mode=="unit") { |
||
284 | werner | 124 | foreach( const ResourceUnit *const_ru, g->model()->ruList()) { |
125 | ResourceUnit *ru = const_cast<ResourceUnit*>(const_ru); |
||
126 | // set environment |
||
127 | g->model()->environment()->setPosition(ru->boundingBox().center()); |
||
128 | type = xml.value("type", ""); |
||
129 | fileName = xml.value("file", ""); |
||
319 | werner | 130 | if (fileName.isEmpty()) |
131 | continue; |
||
549 | werner | 132 | loadInitFile(fileName, type, 0, ru); |
431 | werner | 133 | if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees."; |
284 | werner | 134 | } |
299 | werner | 135 | evaluateDebugTrees(); |
284 | werner | 136 | return; |
281 | werner | 137 | } |
299 | werner | 138 | |
732 | werner | 139 | // map-modus: load a init file for each "polygon" in the standgrid |
549 | werner | 140 | if (copy_mode=="map") { |
141 | if (!g->model()->standGrid() || !g->model()->standGrid()->isValid()) |
||
142 | throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)")); |
||
143 | QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init"); |
||
144 | |||
145 | CSVFile map_file(map_file_name); |
||
146 | if (map_file.rowCount()==0) |
||
147 | throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name)); |
||
148 | int ikey = map_file.columnIndex("id"); |
||
149 | int ivalue = map_file.columnIndex("filename"); |
||
150 | if (ikey<0 || ivalue<0) |
||
151 | throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name)); |
||
152 | QString file_name; |
||
153 | for (int i=0;i<map_file.rowCount();i++) { |
||
154 | int key = map_file.value(i, ikey).toInt(); |
||
155 | if (key>0) { |
||
156 | file_name = map_file.value(i, ivalue).toString(); |
||
550 | werner | 157 | if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key; |
734 | werner | 158 | if (!file_name.isEmpty()) |
159 | loadInitFile(file_name, type, key, NULL); |
||
549 | werner | 160 | } |
161 | } |
||
732 | werner | 162 | mInitHeightGrid = 0; |
549 | werner | 163 | evaluateDebugTrees(); |
164 | return; |
||
165 | } |
||
904 | werner | 166 | |
167 | // standgrid mode: load one large init file |
||
168 | if (copy_mode=="standgrid") { |
||
169 | fileName = GlobalSettings::instance()->path(fileName, "init"); |
||
170 | if (!QFile::exists(fileName)) |
||
171 | throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName)); |
||
172 | QString content = Helper::loadTextFile(fileName); |
||
173 | // this processes the init file (also does the checking) and |
||
174 | // stores in a QHash datastrucutre |
||
175 | parseInitFile(content, fileName); |
||
176 | |||
177 | // setup the random distribution |
||
178 | QString density_func = xml.value("model.initialization.randomFunction", "1-x^2"); |
||
179 | if (logLevelInfo()) qDebug() << "density function:" << density_func; |
||
180 | if (!mRandom || (mRandom->densityFunction()!= density_func)) { |
||
181 | if (mRandom) |
||
182 | delete mRandom; |
||
183 | mRandom=new RandomCustomPDF(density_func); |
||
184 | if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func; |
||
185 | } |
||
186 | |||
187 | if (mStandInitItems.isEmpty()) |
||
188 | throw IException("StandLoader::processInit: 'mode' is 'standgrid' but the init file is either empty or contains no 'stand_id'-column."); |
||
189 | QHash<int, QVector<InitFileItem> >::const_iterator it = mStandInitItems.constBegin(); |
||
190 | while (it!=mStandInitItems.constEnd()) { |
||
191 | mInitItems = it.value(); // copy the items... |
||
192 | executeiLandInitStand(it.key()); |
||
193 | ++it; |
||
194 | } |
||
195 | qDebug() << "finished setup of trees."; |
||
196 | return; |
||
197 | |||
198 | } |
||
699 | werner | 199 | if (copy_mode=="snapshot") { |
200 | // load a snapshot from a file |
||
201 | Snapshot shot; |
||
549 | werner | 202 | |
699 | werner | 203 | QString input_db = GlobalSettings::instance()->path(fileName); |
204 | shot.loadSnapshot(input_db); |
||
205 | return; |
||
206 | } |
||
281 | werner | 207 | throw IException("StandLoader::processInit: invalid initalization.mode!"); |
208 | } |
||
209 | |||
967 | werner | 210 | void StandLoader::processAfterInit() |
211 | { |
||
212 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.initialization")); |
||
213 | |||
214 | QString mode = xml.value("mode", "copy"); |
||
215 | if (mode=="standgrid") { |
||
216 | // load a file with saplings per polygon |
||
217 | QString filename = xml.value("saplingFile", ""); |
||
218 | if (filename.isEmpty()) |
||
219 | return; |
||
220 | filename = GlobalSettings::instance()->path(filename, "init"); |
||
221 | if (!QFile::exists(filename)) |
||
222 | throw IException(QString("load-sapling-ini-file: file '%1' does not exist.").arg(filename)); |
||
223 | CSVFile init_file(filename); |
||
224 | int istandid = init_file.columnIndex("stand_id"); |
||
225 | if (istandid==-1) |
||
226 | throw IException("Sapling-Init: the init file contains no 'stand_id' column (required in 'standgrid' mode)."); |
||
227 | |||
228 | int stand_id = -99999; |
||
229 | int ilow = -1, ihigh = 0; |
||
230 | int total = 0; |
||
231 | for (int i=0;i<init_file.rowCount();++i) { |
||
232 | int row_stand = init_file.value(i, istandid).toInt(); |
||
233 | if (row_stand != stand_id) { |
||
234 | if (stand_id>=0) { |
||
235 | // process stand |
||
236 | ihigh = i-1; // up to the last |
||
237 | total += loadSaplingsLIF(stand_id, init_file, ilow, ihigh); |
||
238 | } |
||
239 | ilow = i; // mark beginning of new stand |
||
240 | stand_id = row_stand; |
||
241 | } |
||
242 | } |
||
243 | if (stand_id>=0) |
||
244 | total += loadSaplingsLIF(stand_id, init_file, ilow, init_file.rowCount()-1); // the last stand |
||
245 | qDebug() << "initialization of sapling: total created:" << total; |
||
246 | |||
247 | } |
||
248 | |||
249 | } |
||
250 | |||
281 | werner | 251 | void StandLoader::evaluateDebugTrees() |
252 | { |
||
137 | Werner | 253 | // evaluate debugging |
254 | QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree"); |
||
299 | werner | 255 | int counter=0; |
137 | Werner | 256 | if (!dbg_str.isEmpty()) { |
733 | werner | 257 | if (dbg_str == "debugstamp") { |
734 | werner | 258 | qDebug() << "debug_tree = debugstamp: try touching all trees..."; |
733 | werner | 259 | // try to force an error if a stamp is invalid |
260 | AllTreeIterator at(GlobalSettings::instance()->model()); |
||
753 | werner | 261 | double total_offset=0.; |
733 | werner | 262 | while (Tree *t=at.next()) { |
263 | total_offset += t->stamp()->offset(); |
||
734 | werner | 264 | if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex())) |
265 | qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!"; |
||
733 | werner | 266 | } |
734 | werner | 267 | qDebug() << "debug_tree = debugstamp: try touching all trees finished..."; |
733 | werner | 268 | return; |
269 | } |
||
270 | TreeWrapper tw; |
||
271 | Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables |
||
137 | Werner | 272 | AllTreeIterator at(GlobalSettings::instance()->model()); |
273 | double result; |
||
274 | while (Tree *t = at.next()) { |
||
140 | Werner | 275 | tw.setTree(t); |
137 | Werner | 276 | result = dexp.execute(); |
299 | werner | 277 | if (result) { |
137 | Werner | 278 | t->enableDebugging(); |
299 | werner | 279 | counter++; |
280 | } |
||
137 | Werner | 281 | } |
732 | werner | 282 | qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees."; |
137 | Werner | 283 | } |
106 | Werner | 284 | } |
285 | |||
904 | werner | 286 | |
732 | werner | 287 | /// load a single init file. Calls loadPicusFile() or loadiLandFile() |
288 | /// @param fileName file to load |
||
289 | /// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution" |
||
549 | werner | 290 | int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru) |
281 | werner | 291 | { |
284 | werner | 292 | QString pathFileName = GlobalSettings::instance()->path(fileName, "init"); |
293 | if (!QFile::exists(pathFileName)) |
||
294 | throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName)); |
||
281 | werner | 295 | |
549 | werner | 296 | if (type=="picus" || type=="single") { |
297 | if (stand_id>0) |
||
298 | throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type); |
||
294 | werner | 299 | return loadPicusFile(pathFileName, ru); |
549 | werner | 300 | } |
384 | werner | 301 | if (type=="iland" || type=="distribution") |
549 | werner | 302 | return loadiLandFile(pathFileName, ru, stand_id); |
281 | werner | 303 | |
384 | werner | 304 | throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type); |
281 | werner | 305 | } |
306 | |||
393 | werner | 307 | int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru) |
106 | Werner | 308 | { |
294 | werner | 309 | QString content = Helper::loadTextFile(fileName); |
310 | if (content.isEmpty()) { |
||
311 | qDebug() << "file not found: " + fileName; |
||
393 | werner | 312 | return 0; |
294 | werner | 313 | } |
393 | werner | 314 | return loadSingleTreeList(content, ru, fileName); |
294 | werner | 315 | } |
316 | |||
389 | werner | 317 | /** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting. |
318 | returns the number of loaded trees. |
||
319 | */ |
||
320 | int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName) |
||
294 | werner | 321 | { |
106 | Werner | 322 | if (!ru) |
323 | ru = mModel->ru(); |
||
284 | werner | 324 | Q_ASSERT(ru!=0); |
106 | Werner | 325 | |
284 | werner | 326 | QPointF offset = ru->boundingBox().topLeft(); |
327 | SpeciesSet *speciesSet = ru->speciesSet(); // of default RU |
||
328 | |||
384 | werner | 329 | QString my_content(content); |
330 | // cut out the <trees> </trees> part if present |
||
331 | if (content.contains("<trees>")) { |
||
332 | QRegExp rx(".*<trees>(.*)</trees>.*"); |
||
333 | rx.indexIn(content, 0); |
||
334 | if (rx.capturedTexts().count()<1) |
||
389 | werner | 335 | return 0; |
384 | werner | 336 | my_content = rx.cap(1).trimmed(); |
337 | } |
||
338 | |||
948 | werner | 339 | CSVFile infile; |
340 | infile.loadFromString(my_content); |
||
384 | werner | 341 | |
138 | Werner | 342 | |
948 | werner | 343 | int iID = infile.columnIndex("id"); |
344 | int iX = infile.columnIndex("x"); |
||
345 | int iY = infile.columnIndex("y"); |
||
346 | int iBhd = infile.columnIndex("bhdfrom"); |
||
384 | werner | 347 | if (iBhd<0) |
948 | werner | 348 | iBhd = infile.columnIndex("dbh"); |
384 | werner | 349 | double height_conversion = 100.; |
948 | werner | 350 | int iHeight =infile.columnIndex("treeheight"); |
384 | werner | 351 | if (iHeight<0) { |
948 | werner | 352 | iHeight = infile.columnIndex("height"); |
384 | werner | 353 | height_conversion = 1.; // in meter |
354 | } |
||
948 | werner | 355 | int iSpecies = infile.columnIndex("species"); |
356 | int iAge = infile.columnIndex("age"); |
||
204 | werner | 357 | if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1) |
948 | werner | 358 | throw IException(QString("Initfile %1 is not valid!\nRequired columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName)); |
106 | Werner | 359 | |
138 | Werner | 360 | double dbh; |
384 | werner | 361 | bool ok; |
389 | werner | 362 | int cnt=0; |
384 | werner | 363 | QString speciesid; |
948 | werner | 364 | for (int i=1;i<infile.rowCount();i++) { |
365 | dbh = infile.value(i, iBhd).toDouble(); |
||
285 | werner | 366 | |
948 | werner | 367 | //if (dbh<5.) |
368 | // continue; |
||
369 | |||
106 | Werner | 370 | QPointF f; |
371 | if (iX>=0 && iY>=0) { |
||
948 | werner | 372 | f.setX( infile.value(i, iX).toDouble() ); |
373 | f.setY( infile.value(i, iY).toDouble() ); |
||
106 | Werner | 374 | f+=offset; |
285 | werner | 375 | |
106 | Werner | 376 | } |
285 | werner | 377 | // position valid? |
378 | if (!mModel->heightGrid()->valueAt(f).isValid()) |
||
379 | continue; |
||
380 | Tree &tree = ru->newTree(); |
||
381 | tree.setPosition(f); |
||
247 | werner | 382 | if (iID>=0) |
948 | werner | 383 | tree.setId(infile.value(i, iID).toInt()); |
106 | Werner | 384 | |
138 | Werner | 385 | tree.setDbh(dbh); |
948 | werner | 386 | tree.setHeight(infile.value(i,iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary |
384 | werner | 387 | |
948 | werner | 388 | speciesid = infile.value(i,iSpecies).toString(); |
270 | werner | 389 | int picusid = speciesid.toInt(&ok); |
390 | if (ok) { |
||
384 | werner | 391 | int idx = picusSpeciesIds.indexOf(picusid); |
270 | werner | 392 | if (idx==-1) |
393 | throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid)); |
||
138 | Werner | 394 | speciesid = iLandSpeciesIds[idx]; |
270 | werner | 395 | } |
138 | Werner | 396 | Species *s = speciesSet->species(speciesid); |
397 | if (!ru || !s) |
||
389 | werner | 398 | throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid)); |
388 | werner | 399 | tree.setSpecies(s); |
107 | Werner | 400 | |
388 | werner | 401 | ok = true; |
402 | if (iAge>=0) |
||
948 | werner | 403 | tree.setAge(infile.value(i, iAge).toInt(&ok), tree.height()); // this is a *real* age |
388 | werner | 404 | if (iAge<0 || !ok || tree.age()==0) |
405 | tree.setAge(0, tree.height()); // no real tree age available |
||
406 | |||
138 | Werner | 407 | tree.setRU(ru); |
106 | Werner | 408 | tree.setup(); |
389 | werner | 409 | cnt++; |
106 | Werner | 410 | } |
389 | werner | 411 | return cnt; |
106 | Werner | 412 | //qDebug() << "loaded init-file contained" << lines.count() <<"lines."; |
413 | //qDebug() << "lines: " << lines; |
||
414 | } |
||
287 | werner | 415 | |
393 | werner | 416 | /** initialize trees on a resource unit based on dbh distributions. |
417 | use a fairly clever algorithm to determine tree positions. |
||
418 | see http://iland.boku.ac.at/initialize+trees |
||
419 | @param content tree init file (including headers) in a string |
||
420 | @param ru resource unit |
||
421 | @param fileName source file name (for error reporting) |
||
422 | @return number of trees added |
||
423 | */ |
||
549 | werner | 424 | int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName) |
287 | werner | 425 | { |
904 | werner | 426 | int total_count = parseInitFile(content, fileName, ru); |
427 | if (total_count==0) |
||
428 | return 0; |
||
429 | |||
430 | |||
431 | // setup the random distribution |
||
432 | QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2"); |
||
433 | if (logLevelInfo()) qDebug() << "density function:" << density_func; |
||
434 | if (!mRandom || (mRandom->densityFunction()!= density_func)) { |
||
435 | if (mRandom) |
||
436 | delete mRandom; |
||
437 | mRandom=new RandomCustomPDF(density_func); |
||
438 | if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func; |
||
439 | } |
||
440 | if (stand_id>0) { |
||
441 | // execute stand based initialization |
||
442 | executeiLandInitStand(stand_id); |
||
443 | } else { |
||
444 | // exeucte the initialization based on single resource units |
||
445 | executeiLandInit(ru); |
||
446 | ru->cleanTreeList(); |
||
447 | } |
||
448 | return total_count; |
||
449 | |||
450 | } |
||
451 | |||
452 | int StandLoader::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru) |
||
453 | { |
||
287 | werner | 454 | if (!ru) |
455 | ru = mModel->ru(); |
||
456 | Q_ASSERT(ru!=0); |
||
457 | SpeciesSet *speciesSet = ru->speciesSet(); // of default RU |
||
458 | Q_ASSERT(speciesSet!=0); |
||
459 | |||
431 | werner | 460 | //DebugTimer t("StandLoader::loadiLandFile"); |
294 | werner | 461 | CSVFile infile; |
462 | infile.loadFromString(content); |
||
463 | |||
287 | werner | 464 | int icount = infile.columnIndex("count"); |
465 | int ispecies = infile.columnIndex("species"); |
||
466 | int idbh_from = infile.columnIndex("dbh_from"); |
||
467 | int idbh_to = infile.columnIndex("dbh_to"); |
||
468 | int ihd = infile.columnIndex("hd"); |
||
469 | int iage = infile.columnIndex("age"); |
||
311 | werner | 470 | int idensity = infile.columnIndex("density"); |
287 | werner | 471 | if (icount<0 || ispecies<0 || idbh_from<0 || idbh_to<0 || ihd<0 || iage<0) |
472 | throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName)); |
||
473 | |||
904 | werner | 474 | int istandid = infile.columnIndex("stand_id"); |
294 | werner | 475 | mInitItems.clear(); |
904 | werner | 476 | mStandInitItems.clear(); |
477 | |||
287 | werner | 478 | InitFileItem item; |
384 | werner | 479 | bool ok; |
393 | werner | 480 | int total_count = 0; |
287 | werner | 481 | for (int row=0;row<infile.rowCount();row++) { |
968 | werner | 482 | item.count = infile.value(row, icount).toDouble(); |
549 | werner | 483 | total_count += item.count; |
484 | item.dbh_from = infile.value(row, idbh_from).toDouble(); |
||
485 | item.dbh_to = infile.value(row, idbh_to).toDouble(); |
||
486 | item.hd = infile.value(row, ihd).toDouble(); |
||
775 | werner | 487 | if (item.hd==0. || item.dbh_from / 100. * item.hd < 4.) |
782 | werner | 488 | qWarning() << QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ; |
489 | //throw IException(QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ); |
||
549 | werner | 490 | ok = true; |
491 | if (iage>=0) |
||
492 | item.age = infile.value(row, iage).toInt(&ok); |
||
493 | if (iage<0 || !ok) |
||
494 | item.age = 0; |
||
384 | werner | 495 | |
549 | werner | 496 | item.species = speciesSet->species(infile.value(row, ispecies).toString()); |
497 | if (idensity>=0) |
||
498 | item.density = infile.value(row, idensity).toDouble(); |
||
499 | else |
||
500 | item.density = 0.; |
||
971 | werner | 501 | if (item.density<-1) |
549 | werner | 502 | throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.") |
503 | .arg(item.density) |
||
504 | .arg(fileName) |
||
505 | .arg(row)); |
||
506 | if (!item.species) { |
||
507 | throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.") |
||
508 | .arg(infile.value(row, ispecies).toString()) |
||
509 | .arg(fileName) |
||
510 | .arg(row)); |
||
511 | } |
||
904 | werner | 512 | if (istandid>=0) { |
513 | int standid = infile.value(row,istandid).toInt(); |
||
514 | mStandInitItems[standid].push_back(item); |
||
515 | } else { |
||
516 | mInitItems.push_back(item); |
||
517 | } |
||
287 | werner | 518 | } |
393 | werner | 519 | return total_count; |
287 | werner | 520 | |
521 | } |
||
522 | |||
904 | werner | 523 | |
549 | werner | 524 | int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id) |
294 | werner | 525 | { |
526 | if (!QFile::exists(fileName)) |
||
527 | throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName)); |
||
528 | QString content = Helper::loadTextFile(fileName); |
||
549 | werner | 529 | return loadDistributionList(content, ru, stand_id, fileName); |
294 | werner | 530 | } |
531 | |||
287 | werner | 532 | // evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions. |
533 | // e.g. 12 = centerpixel, 0: upper left corner, ... |
||
534 | int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4, |
||
535 | 1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9}; |
||
290 | werner | 536 | int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15, |
537 | 6,18,2,10,4,24,12,0,8,14,20,22}; |
||
538 | |||
539 | |||
549 | werner | 540 | |
287 | werner | 541 | // sort function |
542 | bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2) |
||
311 | werner | 543 | { |
544 | return s1.second < s2.second; |
||
545 | } |
||
546 | |||
549 | werner | 547 | struct SInitPixel { |
732 | werner | 548 | double basal_area; // accumulated basal area |
549 | QPoint pixelOffset; // location of the pixel |
||
550 | ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to |
||
551 | double h_max; // predefined maximum height at given pixel (if available from LIDAR or so) |
||
971 | werner | 552 | bool locked; // pixel is dedicated to a single species |
553 | SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.), locked(false) {} |
||
549 | werner | 554 | }; |
555 | |||
556 | bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2) |
||
557 | { |
||
558 | return s1.basal_area < s2.basal_area; |
||
559 | } |
||
560 | |||
971 | werner | 561 | bool sortInitPixelUnlocked(const SInitPixel &s1, const SInitPixel &s2) |
562 | { |
||
563 | return !s1.locked && s2.locked; |
||
564 | } |
||
549 | werner | 565 | |
311 | werner | 566 | /** |
567 | */ |
||
549 | werner | 568 | |
290 | werner | 569 | void StandLoader::executeiLandInit(ResourceUnit *ru) |
287 | werner | 570 | { |
290 | werner | 571 | |
287 | werner | 572 | QPointF offset = ru->boundingBox().topLeft(); |
573 | QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset); |
||
574 | |||
575 | // a multimap holds a list for all trees. |
||
576 | // key is the index of a 10x10m pixel within the resource unit |
||
577 | QMultiMap<int, int> tree_map; |
||
549 | werner | 578 | //QHash<int,SInitPixel> tcount; |
579 | |||
287 | werner | 580 | QVector<QPair<int, double> > tcount; // counts |
581 | for (int i=0;i<100;i++) |
||
582 | tcount.push_back(QPair<int,double>(i,0.)); |
||
583 | |||
584 | int key; |
||
311 | werner | 585 | double rand_val, rand_fraction; |
312 | werner | 586 | int total_count = 0; |
290 | werner | 587 | foreach(const InitFileItem &item, mInitItems) { |
311 | werner | 588 | rand_fraction = fabs(double(item.density)); |
287 | werner | 589 | for (int i=0;i<item.count;i++) { |
590 | // create trees |
||
591 | int tree_idx = ru->newTreeIndex(); |
||
592 | Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree |
||
707 | werner | 593 | tree.setDbh(nrandom(item.dbh_from, item.dbh_to)); |
287 | werner | 594 | tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height |
595 | tree.setSpecies(item.species); |
||
384 | werner | 596 | if (item.age<=0) |
388 | werner | 597 | tree.setAge(0,tree.height()); |
381 | werner | 598 | else |
388 | werner | 599 | tree.setAge(item.age, tree.height()); |
287 | werner | 600 | tree.setRU(ru); |
601 | tree.setup(); |
||
312 | werner | 602 | total_count++; |
311 | werner | 603 | |
604 | // calculate random value. "density" is from 1..-1. |
||
605 | rand_val = mRandom->get(); |
||
606 | if (item.density<0) |
||
607 | rand_val = 1. - rand_val; |
||
705 | werner | 608 | rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction); |
311 | werner | 609 | |
287 | werner | 610 | // key: rank of target pixel |
611 | // first: index of target pixel |
||
612 | // second: sum of target pixel |
||
311 | werner | 613 | key = limit(int(100*rand_val), 0, 99); // get from random number generator |
287 | werner | 614 | tree_map.insert(tcount[key].first, tree_idx); // store tree in map |
615 | tcount[key].second+=tree.basalArea(); // aggregate the basal area for each 10m pixel |
||
312 | werner | 616 | if ( (total_count < 20 && i%2==0) |
617 | || (total_count<100 && i%10==0 ) |
||
618 | || (i%30==0) ) { |
||
287 | werner | 619 | qSort(tcount.begin(), tcount.end(), sortPairLessThan); |
312 | werner | 620 | } |
287 | werner | 621 | } |
622 | qSort(tcount.begin(), tcount.end(), sortPairLessThan); |
||
623 | } |
||
624 | |||
625 | int bits, index, pos; |
||
626 | int c; |
||
627 | QList<int> trees; |
||
628 | QPoint tree_pos; |
||
289 | werner | 629 | |
287 | werner | 630 | for (int i=0;i<100;i++) { |
631 | trees = tree_map.values(i); |
||
632 | c = trees.count(); |
||
290 | werner | 633 | QPointF pixel_center = ru->boundingBox().topLeft() + QPointF((i/10)*10. + 5., (i%10)*10. + 5.); |
634 | if (!mModel->heightGrid()->valueAt(pixel_center).isValid()) { |
||
635 | // no trees on that pixel: let trees die |
||
636 | foreach(int tree_idx, trees) { |
||
637 | ru->trees()[tree_idx].die(); |
||
638 | } |
||
639 | continue; |
||
640 | } |
||
641 | |||
287 | werner | 642 | bits = 0; |
643 | index = -1; |
||
290 | werner | 644 | double r; |
287 | werner | 645 | foreach(int tree_idx, trees) { |
646 | if (c>18) { |
||
647 | index = (index + 1)%25; |
||
648 | } else { |
||
649 | int stop=1000; |
||
290 | werner | 650 | index = 0; |
287 | werner | 651 | do { |
290 | werner | 652 | //r = drandom(); |
653 | //if (r<0.5) // skip position with a prob. of 50% -> adds a little "noise" |
||
654 | // index++; |
||
655 | //index = (index + 1)%25; // increase and roll over |
||
656 | |||
657 | // search a random position |
||
658 | r = drandom(); |
||
659 | index = limit(int(25 * r*r), 0, 24); // use rnd()^2 to search for locations -> higher number of low indices (i.e. 50% of lookups in first 25% of locations) |
||
287 | werner | 660 | } while (isBitSet(bits, index)==true && stop--); |
661 | if (!stop) |
||
662 | qDebug() << "executeiLandInit: found no free bit."; |
||
663 | setBit(bits, index, true); // mark position as used |
||
664 | } |
||
290 | werner | 665 | // get position from fixed lists (one for even, one for uneven resource units) |
666 | pos = ru->index()%2?evenlist[index]:unevenlist[index]; |
||
287 | werner | 667 | tree_pos = offsetIdx // position of resource unit |
668 | + QPoint(5*(i/10), 5*(i%10)) // relative position of 10x10m pixel |
||
669 | + QPoint(pos/5, pos%5); // relative position within 10x10m pixel |
||
670 | //qDebug() << tree_no++ << "to" << index; |
||
671 | ru->trees()[tree_idx].setPosition(tree_pos); |
||
672 | } |
||
673 | } |
||
674 | } |
||
549 | werner | 675 | |
676 | |||
550 | werner | 677 | |
678 | // Initialization routine based on a stand map. |
||
679 | // Basically a list of 10m pixels for a given stand is retrieved |
||
680 | // and the filled with the same procedure as the resource unit based init |
||
681 | // see http://iland.boku.ac.at/initialize+trees |
||
549 | werner | 682 | void StandLoader::executeiLandInitStand(int stand_id) |
683 | { |
||
684 | |||
685 | const MapGrid *grid = GlobalSettings::instance()->model()->standGrid(); |
||
901 | werner | 686 | if (mCurrentMap) |
687 | grid = mCurrentMap; |
||
549 | werner | 688 | |
732 | werner | 689 | // get a list of positions of all pixels that belong to our stand |
549 | werner | 690 | QList<int> indices = grid->gridIndices(stand_id); |
550 | werner | 691 | if (indices.isEmpty()) { |
692 | qDebug() << "stand" << stand_id << "not in project area. No init performed."; |
||
693 | return; |
||
549 | werner | 694 | } |
695 | // a multiHash holds a list for all trees. |
||
696 | // key is the location of the 10x10m pixel |
||
697 | QMultiHash<QPoint, int> tree_map; |
||
698 | QList<SInitPixel> pixel_list; // working list of all 10m pixels |
||
732 | werner | 699 | pixel_list.reserve(indices.size()); |
549 | werner | 700 | |
701 | foreach (int i, indices) { |
||
702 | SInitPixel p; |
||
703 | p.pixelOffset = grid->grid().indexOf(i); // index in the 10m grid |
||
732 | werner | 704 | p.resource_unit = GlobalSettings::instance()->model()->ru( grid->grid().cellCenterPoint(p.pixelOffset)); |
705 | if (mInitHeightGrid) |
||
706 | p.h_max = mInitHeightGrid->grid().constValueAtIndex(p.pixelOffset); |
||
549 | werner | 707 | pixel_list.append(p); |
708 | } |
||
709 | double area_factor = grid->area(stand_id) / 10000.; |
||
710 | |||
732 | werner | 711 | int key=0; |
549 | werner | 712 | double rand_val, rand_fraction; |
713 | int total_count = 0; |
||
732 | werner | 714 | int total_tries = 0; |
715 | int total_misses = 0; |
||
716 | if (mInitHeightGrid && !mHeightGridResponse) |
||
717 | throw IException("executeiLandInitStand: trying to initialize with height grid but without response function."); |
||
971 | werner | 718 | |
719 | Species *last_locked_species=0; |
||
549 | werner | 720 | foreach(const InitFileItem &item, mInitItems) { |
971 | werner | 721 | if (item.density>1.) { |
722 | // special case with single-species-area |
||
723 | if (total_count==0) { |
||
724 | // randomize the pixels |
||
725 | for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it) |
||
726 | it->basal_area = drandom(); |
||
727 | qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan); |
||
728 | for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it) |
||
729 | it->basal_area = 0.; |
||
730 | } |
||
731 | |||
732 | if (item.species != last_locked_species) { |
||
733 | last_locked_species=item.species; |
||
734 | qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelUnlocked); |
||
735 | } |
||
736 | } else { |
||
737 | qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan); |
||
738 | last_locked_species=0; |
||
739 | } |
||
740 | rand_fraction = item.density; |
||
549 | werner | 741 | int count = item.count * area_factor + 0.5; // round |
732 | werner | 742 | double init_max_height = item.dbh_to/100. * item.hd; |
549 | werner | 743 | for (int i=0;i<count;i++) { |
744 | |||
732 | werner | 745 | bool found = false; |
746 | int tries = mHeightGridTries; |
||
747 | while (!found &&--tries) { |
||
748 | // calculate random value. "density" is from 1..-1. |
||
971 | werner | 749 | if (item.density <= 1.) { |
750 | rand_val = mRandom->get(); |
||
751 | if (item.density<0) |
||
752 | rand_val = 1. - rand_val; |
||
753 | rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction); |
||
754 | } else { |
||
755 | // limited area: limit potential area using the "density" input parameter |
||
756 | rand_val = drandom() * qMin(item.density/100., 1.); |
||
757 | } |
||
732 | werner | 758 | ++total_tries; |
549 | werner | 759 | |
732 | werner | 760 | // key: rank of target pixel |
761 | key = limit(int(pixel_list.count()*rand_val), 0, pixel_list.count()-1); // get from random number generator |
||
549 | werner | 762 | |
732 | werner | 763 | if (mInitHeightGrid) { |
764 | // calculate how good the selected pixel fits w.r.t. the predefined height |
||
765 | double p_value = pixel_list[key].h_max>0.?mHeightGridResponse->calculate(init_max_height/pixel_list[key].h_max):0.; |
||
766 | if (drandom() < p_value) |
||
767 | found = true; |
||
768 | } else { |
||
769 | found = true; |
||
770 | } |
||
971 | werner | 771 | if (!last_locked_species && pixel_list[key].locked) |
772 | found = false; |
||
732 | werner | 773 | } |
774 | if (tries<0) ++total_misses; |
||
549 | werner | 775 | |
776 | // create a tree |
||
777 | ResourceUnit *ru = pixel_list[key].resource_unit; |
||
778 | int tree_idx = ru->newTreeIndex(); |
||
779 | Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree |
||
780 | tree.setDbh(nrandom(item.dbh_from, item.dbh_to)); |
||
781 | tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height |
||
782 | tree.setSpecies(item.species); |
||
783 | if (item.age<=0) |
||
784 | tree.setAge(0,tree.height()); |
||
785 | else |
||
786 | tree.setAge(item.age, tree.height()); |
||
787 | tree.setRU(ru); |
||
788 | tree.setup(); |
||
789 | total_count++; |
||
790 | |||
791 | // store in the multiHash the position of the pixel and the tree_idx in the resepctive resource unit |
||
792 | tree_map.insert(pixel_list[key].pixelOffset, tree_idx); |
||
793 | pixel_list[key].basal_area+=tree.basalArea(); // aggregate the basal area for each 10m pixel |
||
971 | werner | 794 | if (last_locked_species) |
795 | pixel_list[key].locked = true; |
||
549 | werner | 796 | |
797 | // resort list |
||
971 | werner | 798 | if (last_locked_species==0 && ((total_count < 20 && i%2==0) |
549 | werner | 799 | || (total_count<100 && i%10==0 ) |
971 | werner | 800 | || (i%30==0)) ) { |
549 | werner | 801 | qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan); |
802 | } |
||
803 | } |
||
804 | } |
||
732 | werner | 805 | if (total_misses>0 || total_tries > total_count) { |
744 | werner | 806 | if (logLevelInfo()) qDebug() << "init for stand" << stand_id << "treecount:" << total_count << ", tries:" << total_tries << ", misses:" << total_misses << ", %miss:" << qRound(total_misses*100 / (double)total_count); |
732 | werner | 807 | } |
549 | werner | 808 | |
809 | int bits, index, pos; |
||
810 | int c; |
||
811 | QList<int> trees; |
||
812 | QPoint tree_pos; |
||
813 | |||
814 | foreach(const SInitPixel &p, pixel_list) { |
||
815 | trees = tree_map.values(p.pixelOffset); |
||
816 | c = trees.count(); |
||
817 | bits = 0; |
||
818 | index = -1; |
||
819 | double r; |
||
820 | foreach(int tree_idx, trees) { |
||
821 | if (c>18) { |
||
822 | index = (index + 1)%25; |
||
823 | } else { |
||
824 | int stop=1000; |
||
825 | index = 0; |
||
826 | do { |
||
827 | // search a random position |
||
828 | r = drandom(); |
||
829 | index = limit(int(25 * r*r), 0, 24); // use rnd()^2 to search for locations -> higher number of low indices (i.e. 50% of lookups in first 25% of locations) |
||
830 | } while (isBitSet(bits, index)==true && stop--); |
||
831 | if (!stop) |
||
832 | qDebug() << "executeiLandInit: found no free bit."; |
||
833 | setBit(bits, index, true); // mark position as used |
||
834 | } |
||
835 | // get position from fixed lists (one for even, one for uneven resource units) |
||
836 | pos = p.resource_unit->index()%2?evenlist[index]:unevenlist[index]; |
||
837 | tree_pos = p.pixelOffset * cPxPerHeight; // convert to LIF index |
||
600 | werner | 838 | tree_pos += QPoint(pos/cPxPerHeight, pos%cPxPerHeight); |
550 | werner | 839 | |
549 | werner | 840 | p.resource_unit->trees()[tree_idx].setPosition(tree_pos); |
734 | werner | 841 | // test if tree position is valid.. |
842 | if (!GlobalSettings::instance()->model()->grid()->isIndexValid(tree_pos)) |
||
843 | qDebug() << "Standloader: invalid position!"; |
||
549 | werner | 844 | } |
845 | } |
||
971 | werner | 846 | if (logLevelInfo()) |
847 | qDebug() << "init for stand" << stand_id << "with area" << "area (m2)" << grid->area(stand_id) << "count of 10m pixels:" << indices.count() << "initialized trees:" << total_count; |
||
550 | werner | 848 | |
549 | werner | 849 | } |
600 | werner | 850 | |
851 | /// a (hacky) way of adding saplings of a certain age to a stand defined by 'stand_id'. |
||
852 | int StandLoader::loadSaplings(const QString &content, int stand_id, const QString &fileName) |
||
853 | { |
||
777 | werner | 854 | Q_UNUSED(fileName); |
603 | werner | 855 | const MapGrid *stand_grid; |
856 | if (mCurrentMap) |
||
857 | stand_grid = mCurrentMap; // if set |
||
858 | else |
||
859 | stand_grid = GlobalSettings::instance()->model()->standGrid(); // default |
||
860 | |||
600 | werner | 861 | QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels |
862 | if (indices.isEmpty()) { |
||
863 | qDebug() << "stand" << stand_id << "not in project area. No init performed."; |
||
864 | return -1; |
||
865 | } |
||
866 | double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare) |
||
867 | |||
868 | // parse the content of the init-file |
||
869 | // species |
||
870 | CSVFile init; |
||
871 | init.loadFromString(content); |
||
872 | int ispecies = init.columnIndex("species"); |
||
873 | int icount = init.columnIndex("count"); |
||
874 | int iheight = init.columnIndex("height"); |
||
875 | int iage = init.columnIndex("age"); |
||
876 | if (ispecies==-1 || icount==-1) |
||
877 | throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!"); |
||
878 | |||
879 | const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet(); |
||
880 | double height, age; |
||
881 | int total = 0; |
||
882 | for (int row=0;row<init.rowCount();++row) { |
||
883 | int pxcount = init.value(row, icount).toDouble() * area_factor + 0.5; // no. of pixels that should be filled (sapling grid is the same resolution as the lif-grid) |
||
884 | const Species *species = set->species(init.value(row, ispecies).toString()); |
||
885 | if (!species) |
||
886 | throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString())); |
||
887 | height = iheight==-1?0.05: init.value(row, iheight).toDouble(); |
||
888 | age = iage==-1?1:init.value(row,iage).toDouble(); |
||
889 | |||
890 | int misses = 0; |
||
891 | int hits = 0; |
||
892 | while (hits < pxcount) { |
||
893 | int rnd_index = irandom(0, indices.count()-1); |
||
894 | QPoint offset=stand_grid->grid().indexOf(indices[rnd_index]); |
||
895 | ResourceUnit *ru = GlobalSettings::instance()->model()->ru(stand_grid->grid().cellCenterPoint(offset)); |
||
896 | // |
||
897 | offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates |
||
898 | int in_p = irandom(0, cPxPerHeight*cPxPerHeight-1); // index of lif-pixel |
||
899 | offset += QPoint(in_p / cPxPerHeight, in_p % cPxPerHeight); |
||
900 | if (ru->saplingHeightForInit(offset) > height) { |
||
901 | misses++; |
||
902 | } else { |
||
903 | // ok |
||
904 | hits++; |
||
951 | werner | 905 | ru->resourceUnitSpecies(species).changeSapling().addSapling(offset, height, age); |
600 | werner | 906 | } |
907 | if (misses > 3*pxcount) { |
||
908 | qDebug() << "tried to add" << pxcount << "saplings at stand" << stand_id << "but failed in finding enough free positions. Added" << hits << "and stopped."; |
||
909 | break; |
||
910 | } |
||
911 | } |
||
912 | total += hits; |
||
913 | |||
914 | } |
||
915 | return total; |
||
916 | } |
||
966 | werner | 917 | |
967 | werner | 918 | bool LIFValueHigher(const float *a, const float *b) |
966 | werner | 919 | { |
967 | werner | 920 | return *a > *b; |
966 | werner | 921 | } |
922 | |||
967 | werner | 923 | int StandLoader::loadSaplingsLIF(int stand_id, const CSVFile &init, int low_index, int high_index) |
966 | werner | 924 | { |
925 | const MapGrid *stand_grid; |
||
926 | if (mCurrentMap) |
||
927 | stand_grid = mCurrentMap; // if set |
||
928 | else |
||
929 | stand_grid = GlobalSettings::instance()->model()->standGrid(); // default |
||
930 | |||
967 | werner | 931 | if (!stand_grid->isValid(stand_id)) |
932 | return 0; |
||
933 | |||
966 | werner | 934 | QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels |
935 | if (indices.isEmpty()) { |
||
936 | qDebug() << "stand" << stand_id << "not in project area. No init performed."; |
||
967 | werner | 937 | return 0; |
966 | werner | 938 | } |
939 | // prepare space for LIF-pointers (2m Pixel) |
||
940 | QVector<float*> lif_ptrs; |
||
967 | werner | 941 | lif_ptrs.reserve(indices.size() * cPxPerHeight * cPxPerHeight); |
966 | werner | 942 | for (int l=0;l<indices.size();++l){ |
943 | QPoint offset=stand_grid->grid().indexOf(indices[l]); |
||
944 | offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates |
||
967 | werner | 945 | for (int y=0;y<cPxPerHeight;++y) |
946 | for(int x=0;x<cPxPerHeight;++x) |
||
947 | lif_ptrs.push_back( GlobalSettings::instance()->model()->grid()->ptr(offset.x()+x, offset.y()+y) ); |
||
966 | werner | 948 | } |
949 | // sort based on LIF-Value |
||
967 | werner | 950 | std::sort(lif_ptrs.begin(), lif_ptrs.end(), LIFValueHigher); // higher: highest values first |
966 | werner | 951 | |
952 | |||
953 | double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare) |
||
954 | |||
955 | // parse the content of the init-file |
||
956 | // species |
||
957 | int ispecies = init.columnIndex("species"); |
||
958 | int icount = init.columnIndex("count"); |
||
959 | int iheight = init.columnIndex("height"); |
||
960 | int iheightfrom = init.columnIndex("height_from"); |
||
961 | int iheightto = init.columnIndex("height_to"); |
||
962 | int iage = init.columnIndex("age"); |
||
963 | int itopage = init.columnIndex("age4m"); |
||
964 | int iminlif = init.columnIndex("min_lif"); |
||
965 | if ((iheightfrom==-1) ^ (iheightto==-1)) |
||
966 | throw IException("Error while loading saplings: height not correctly provided. Use either 'height' or 'height_from' and 'height_to'."); |
||
967 | if (ispecies==-1 || icount==-1) |
||
968 | throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!"); |
||
969 | |||
970 | const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet(); |
||
971 | double height, age; |
||
972 | int total = 0; |
||
967 | werner | 973 | for (int row=low_index;row<=high_index;++row) { |
974 | double pxcount = init.value(row, icount).toDouble() * area_factor; // no. of pixels that should be filled (sapling grid is the same resolution as the lif-grid) |
||
966 | werner | 975 | const Species *species = set->species(init.value(row, ispecies).toString()); |
976 | if (!species) |
||
977 | throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString())); |
||
978 | height = iheight==-1?0.05: init.value(row, iheight).toDouble(); |
||
979 | age = iage==-1?1:init.value(row,iage).toDouble(); |
||
980 | double age4m = itopage==-1?10:init.value(row, itopage).toDouble(); |
||
981 | double height_from = iheightfrom==-1?-1.: init.value(row, iheightfrom).toDouble(); |
||
982 | double height_to = iheightto==-1?-1.: init.value(row, iheightto).toDouble(); |
||
983 | double min_lif = iminlif==-1?1.: init.value(row, iminlif).toDouble(); |
||
984 | // find LIF-level in the pixels |
||
985 | int min_lif_index = 0; |
||
986 | for (QVector<float*>::ConstIterator it=lif_ptrs.constBegin(); it!=lif_ptrs.constEnd(); ++it, ++min_lif_index) |
||
967 | werner | 987 | if (**it <= min_lif) |
966 | werner | 988 | break; |
989 | |||
990 | if (min_lif_index < pxcount) { |
||
991 | // not enough LIF pixels available |
||
992 | min_lif_index = pxcount; // try the brightest pixels (ie with the largest value for the LIF) |
||
993 | } |
||
994 | |||
967 | werner | 995 | double hits = 0.; |
966 | werner | 996 | while (hits < pxcount) { |
997 | int rnd_index = irandom(0, min_lif_index-1); |
||
998 | if (iheightfrom!=-1) { |
||
999 | height = limit(nrandom(height_from, height_to), 0.05,4.); |
||
1000 | age = qMax(qRound(height/4. * age4m),1); // assume a linear relationship between height and age |
||
1001 | } |
||
1002 | QPoint offset = GlobalSettings::instance()->model()->grid()->indexOf(lif_ptrs[rnd_index]); |
||
1003 | |||
1004 | ResourceUnit *ru = GlobalSettings::instance()->model()->ru(GlobalSettings::instance()->model()->grid()->cellCenterPoint(offset)); |
||
1005 | ru->resourceUnitSpecies(species).changeSapling().addSapling(offset, height, age); |
||
967 | werner | 1006 | hits += ru->resourceUnitSpecies(species).sapling().representedStemNumber(height); |
1007 | |||
1008 | |||
966 | werner | 1009 | } |
1010 | total += pxcount; |
||
1011 | |||
1012 | } |
||
1013 | return total; |
||
1014 | } |