Rev 1104 | Rev 1167 | 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 | |||
373 | werner | 21 | #include "seeddispersal.h" |
22 | |||
23 | #include "globalsettings.h" |
||
24 | #include "model.h" |
||
808 | werner | 25 | #include "debugtimer.h" |
373 | werner | 26 | #include "helper.h" |
391 | werner | 27 | #include "species.h" |
989 | werner | 28 | #ifdef ILAND_GUI |
373 | werner | 29 | #include <QtGui/QImage> |
989 | werner | 30 | #endif |
373 | werner | 31 | |
32 | /** @class SeedDispersal |
||
697 | werner | 33 | @ingroup core |
373 | werner | 34 | The class encapsulates the dispersal of seeds of one species over the whole landscape. |
697 | werner | 35 | The dispersal algortihm operate on grids with a 20m resolution. |
373 | werner | 36 | |
697 | werner | 37 | See http://iland.boku.ac.at/dispersal |
38 | |||
373 | werner | 39 | */ |
764 | werner | 40 | |
41 | Grid<float> *SeedDispersal::mExternalSeedBaseMap = 0; |
||
42 | QHash<QString, QVector<double> > SeedDispersal::mExtSeedData; |
||
43 | int SeedDispersal::mExtSeedSizeX = 0; |
||
44 | int SeedDispersal::mExtSeedSizeY = 0; |
||
45 | |||
373 | werner | 46 | SeedDispersal::~SeedDispersal() |
47 | { |
||
48 | if (isSetup()) { |
||
49 | |||
50 | } |
||
51 | } |
||
52 | |||
391 | werner | 53 | // ************ Setup ************** |
54 | |||
55 | /** setup of the seedmaps. |
||
56 | This sets the size of the seed map and creates the seed kernel (species specific) |
||
57 | */ |
||
373 | werner | 58 | void SeedDispersal::setup() |
59 | { |
||
391 | werner | 60 | if (!GlobalSettings::instance()->model() |
61 | || !GlobalSettings::instance()->model()->heightGrid() |
||
62 | || !mSpecies) |
||
373 | werner | 63 | return; |
391 | werner | 64 | |
65 | const float seedmap_size = 20.f; |
||
373 | werner | 66 | // setup of seed map |
67 | mSeedMap.clear(); |
||
391 | werner | 68 | mSeedMap.setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
373 | werner | 69 | mSeedMap.initialize(0.); |
764 | werner | 70 | mExternalSeedMap.clear(); |
391 | werner | 71 | mIndexFactor = int(seedmap_size) / cPxSize; // ratio seed grid / lip-grid: |
550 | werner | 72 | if (logLevelInfo()) qDebug() << "Seed map setup. Species:"<< mSpecies->id() << "kernel-size: " << mSeedMap.sizeX() << "x" << mSeedMap.sizeY() << "pixels."; |
373 | werner | 73 | |
445 | werner | 74 | if (mSpecies==0) |
75 | throw IException("Setup of SeedDispersal: Species not defined."); |
||
76 | |||
802 | werner | 77 | if (fmod(GlobalSettings::instance()->settings().valueDouble("model.world.buffer",0),seedmap_size) != 0.) |
78 | throw IException("SeedDispersal:setup(): The buffer (model.world.buffer) must be a integer multiple of the seed pixel size (currently 20m, e.g. 20,40,60,...))."); |
||
79 | |||
415 | werner | 80 | // settings |
445 | werner | 81 | mTM_occupancy = 1.; // is currently constant |
82 | mSpecies->treeMigKernel(mTM_as1, mTM_as2, mTM_ks); // set values |
||
83 | mTM_fecundity_cell = mSpecies->fecundity_m2() * seedmap_size*seedmap_size * mTM_occupancy; // scale to production for the whole cell |
||
84 | mNonSeedYearFraction = mSpecies->nonSeedYearFraction(); |
||
415 | werner | 85 | |
445 | werner | 86 | createKernel(mKernelSeedYear, mTM_fecundity_cell); |
415 | werner | 87 | |
88 | // the kernel for non seed years looks similar, but is simply linearly scaled down |
||
89 | // using the species parameter NonSeedYearFraction. |
||
90 | // the central pixel still gets the value of 1 (i.e. 100% probability) |
||
445 | werner | 91 | createKernel(mKernelNonSeedYear, mTM_fecundity_cell*mNonSeedYearFraction); |
415 | werner | 92 | |
472 | werner | 93 | // debug info |
94 | mDumpSeedMaps = GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false); |
||
95 | if (mDumpSeedMaps) { |
||
1102 | werner | 96 | QString path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 97 | Helper::saveToTextFile(QString("%1/seedkernelYes_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelSeedYear)); |
98 | Helper::saveToTextFile(QString("%1/seedkernelNo_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelNonSeedYear)); |
||
417 | werner | 99 | } |
472 | werner | 100 | // external seeds |
481 | werner | 101 | mHasExternalSeedInput = false; |
491 | werner | 102 | mExternalSeedBuffer = 0; |
103 | mExternalSeedDirection = 0; |
||
836 | werner | 104 | mExternalSeedBackgroundInput = 0.; |
472 | werner | 105 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.externalSeedEnabled",false)) { |
764 | werner | 106 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) { |
107 | // external seed input specified by sectors and around the project area (seedbelt) |
||
108 | setupExternalSeedsForSpecies(mSpecies); |
||
109 | } else { |
||
110 | // external seeds specified fixedly per cardinal direction |
||
111 | // current species in list?? |
||
112 | mHasExternalSeedInput = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSpecies").contains(mSpecies->id()); |
||
113 | QString dir = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSource").toLower(); |
||
114 | // encode cardinal positions as bits: e.g: "e,w" -> 6 |
||
115 | mExternalSeedDirection += dir.contains("n")?1:0; |
||
116 | mExternalSeedDirection += dir.contains("e")?2:0; |
||
117 | mExternalSeedDirection += dir.contains("s")?4:0; |
||
118 | mExternalSeedDirection += dir.contains("w")?8:0; |
||
837 | werner | 119 | QStringList buffer_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBuffer").split(QRegExp("([^\\.\\w]+)")); |
764 | werner | 120 | int index = buffer_list.indexOf(mSpecies->id()); |
121 | if (index>=0) { |
||
122 | mExternalSeedBuffer = buffer_list[index+1].toInt(); |
||
123 | qDebug() << "enabled special buffer for species" <<mSpecies->id() << ": distance of" << mExternalSeedBuffer << "pixels = " << mExternalSeedBuffer*20. << "m"; |
||
124 | } |
||
836 | werner | 125 | |
126 | // background seed rain (i.e. for the full landscape), use regexp |
||
837 | werner | 127 | QStringList background_input_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBackgroundInput").split(QRegExp("([^\\.\\w]+)")); |
836 | werner | 128 | index = background_input_list.indexOf(mSpecies->id()); |
129 | if (index>=0) { |
||
130 | mExternalSeedBackgroundInput = background_input_list[index+1].toDouble(); |
||
131 | qDebug() << "enabled background seed input (for full area) for species" <<mSpecies->id() << ": p=" << mExternalSeedBackgroundInput; |
||
132 | } |
||
133 | |||
764 | werner | 134 | if (mHasExternalSeedInput) |
135 | qDebug() << "External seed input enabled for" << mSpecies->id(); |
||
491 | werner | 136 | } |
472 | werner | 137 | } |
415 | werner | 138 | |
373 | werner | 139 | // setup of seed kernel |
391 | werner | 140 | // const int max_radius = 15; // pixels |
141 | // |
||
142 | // mSeedKernel.clear(); |
||
143 | // mSeedKernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
||
144 | // mKernelOffset = max_radius; |
||
145 | // // filling of the kernel.... for simplicity: a linear kernel |
||
146 | // QPoint center = QPoint(mKernelOffset, mKernelOffset); |
||
147 | // const double max_dist = max_radius * seedmap_size; |
||
148 | // for (float *p=mSeedKernel.begin(); p!=mSeedKernel.end();++p) { |
||
149 | // double d = mSeedKernel.distance(center, mSeedKernel.indexOf(p)); |
||
150 | // *p = qMax( 1. - d / max_dist, 0.); |
||
151 | // } |
||
373 | werner | 152 | |
153 | |||
154 | // randomize seed map.... set 1/3 to "filled" |
||
375 | werner | 155 | //for (int i=0;i<mSeedMap.count(); i++) |
156 | // mSeedMap.valueAtIndex(mSeedMap.randomPosition()) = 1.; |
||
373 | werner | 157 | |
158 | |||
375 | werner | 159 | // QImage img = gridToImage(mSeedMap, true, -1., 1.); |
160 | // img.save("seedmap.png"); |
||
161 | |||
162 | // img = gridToImage(mSeedMap, true, -1., 1.); |
||
764 | werner | 163 | // img.save("seedmap_e.png"); |
373 | werner | 164 | } |
165 | |||
764 | werner | 166 | void SeedDispersal::setupExternalSeeds() |
167 | { |
||
168 | mExternalSeedBaseMap = 0; |
||
169 | if (!GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) |
||
170 | return; |
||
171 | |||
978 | werner | 172 | DebugTimer t("setup of external seed maps."); |
764 | werner | 173 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.settings.seedDispersal.seedBelt")); |
1102 | werner | 174 | int seedbelt_width =xml.valueInt(".width",10); |
764 | werner | 175 | // setup of sectors |
176 | // setup of base map |
||
177 | const float seedmap_size = 20.f; |
||
178 | mExternalSeedBaseMap = new Grid<float>; |
||
179 | mExternalSeedBaseMap->setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
||
180 | mExternalSeedBaseMap->initialize(0.); |
||
181 | if (mExternalSeedBaseMap->count()*4 != GlobalSettings::instance()->model()->heightGrid()->count()) |
||
947 | werner | 182 | throw IException("error in setting up external seeds: the width and height of the project area need to be a multiple of 20m when external seeds are enabled."); |
764 | werner | 183 | // make a copy of the 10m height grid in lower resolution and mark pixels that are forested and outside of |
184 | // the project area. |
||
185 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) |
||
186 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
765 | werner | 187 | bool val = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isForestOutside(); |
764 | werner | 188 | mExternalSeedBaseMap->valueAtIndex(x,y) = val?1.f:0.f; |
189 | if(GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isValid()) |
||
190 | mExternalSeedBaseMap->valueAtIndex(x,y) = -1.f; |
||
191 | } |
||
836 | werner | 192 | QString path = GlobalSettings::instance()->path(GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath")); |
765 | werner | 193 | |
194 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
||
989 | werner | 195 | #ifdef ILAND_GUI |
765 | werner | 196 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
197 | img.save(path + "/seedbeltmap_before.png"); |
||
989 | werner | 198 | #endif |
765 | werner | 199 | } |
764 | werner | 200 | // img.save("seedmap.png"); |
201 | // now scan the pixels of the belt: paint all pixels that are close to the project area |
||
202 | // we do this 4 times (for all cardinal direcitons) |
||
203 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
204 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
205 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
206 | continue; |
||
765 | werner | 207 | int look_forward = std::min(x + seedbelt_width, mExternalSeedBaseMap->sizeX()-1); |
1106 | werner | 208 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 209 | // fill pixels |
210 | for(; x<look_forward;++x) { |
||
211 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
212 | if (v==1.f) v=2.f; |
||
213 | } |
||
214 | } |
||
215 | } |
||
216 | } |
||
217 | // right to left |
||
218 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
219 | for (int x=mExternalSeedBaseMap->sizeX();x>=0;--x) { |
||
220 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
221 | continue; |
||
222 | int look_forward = std::max(x - seedbelt_width, 0); |
||
1106 | werner | 223 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 224 | // fill pixels |
225 | for(; x>look_forward;--x) { |
||
226 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
227 | if (v==1.f) v=2.f; |
||
228 | } |
||
229 | } |
||
230 | } |
||
231 | } |
||
232 | // up and down *** |
||
233 | // from top to bottom |
||
234 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
235 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
236 | |||
237 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
238 | continue; |
||
765 | werner | 239 | int look_forward = std::min(y + seedbelt_width, mExternalSeedBaseMap->sizeY()-1); |
764 | werner | 240 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
241 | // fill pixels |
||
242 | for(; y<look_forward;++y) { |
||
243 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
244 | if (v==1.f) v=2.f; |
||
245 | } |
||
246 | } |
||
247 | } |
||
248 | } |
||
249 | // bottom to top *** |
||
250 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
251 | for (int x=mExternalSeedBaseMap->sizeX();x>=0;--x) { |
||
252 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
253 | continue; |
||
254 | int look_forward = std::max(y - seedbelt_width, 0); |
||
255 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
||
256 | // fill pixels |
||
257 | for(; y>look_forward;--y) { |
||
258 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
259 | if (v==1.f) v=2.f; |
||
260 | } |
||
261 | } |
||
262 | } |
||
263 | } |
||
765 | werner | 264 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
989 | werner | 265 | #ifdef ILAND_GUI |
765 | werner | 266 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
972 | werner | 267 | img.save(path + "/seedbeltmap_after.png"); |
989 | werner | 268 | #endif |
765 | werner | 269 | } |
764 | werner | 270 | mExtSeedData.clear(); |
1102 | werner | 271 | int sectors_x = xml.valueInt("sizeX",0); |
272 | int sectors_y = xml.valueInt("sizeY",0); |
||
764 | werner | 273 | if(sectors_x<1 || sectors_y<1) |
274 | throw IException(QString("setup of external seed dispersal: invalid number of sectors x=%1 y=%3").arg(sectors_x).arg(sectors_y)); |
||
275 | QDomElement elem = xml.node("."); |
||
276 | for(QDomNode n = elem.firstChild(); !n.isNull(); n = n.nextSibling()) { |
||
277 | if (n.nodeName().startsWith("species")) { |
||
278 | QStringList coords = n.nodeName().split("_"); |
||
279 | if (coords.count()!=3) |
||
280 | throw IException("external seed species definition is not valid: " + n.nodeName()); |
||
281 | int x = coords[1].toInt(); |
||
282 | int y = coords[2].toInt(); |
||
283 | if (x<0 || x>=sectors_x || y<0 || y>=sectors_y) |
||
284 | throw IException(QString("invalid sector for specifiing external seed input (x y): %1 %2 ").arg(x).arg(y) ); |
||
285 | int index = y*sectors_x + x; |
||
286 | |||
287 | QString text = xml.value("." + n.nodeName()); |
||
288 | qDebug() << "processing element " << n.nodeName() << "x,y:" << x << y << text; |
||
289 | // we assume pairs of name and fraction |
||
290 | QStringList species_list = text.split(" "); |
||
291 | for (int i=0;i<species_list.count();++i) { |
||
292 | QVector<double> &space = mExtSeedData[species_list[i]]; |
||
293 | if (space.isEmpty()) |
||
294 | space.resize(sectors_x*sectors_y); // are initialized to 0s |
||
295 | double fraction = species_list[++i].toDouble(); |
||
296 | space[index] = fraction; |
||
297 | } |
||
298 | } |
||
299 | } |
||
300 | mExtSeedSizeX = sectors_x; |
||
301 | mExtSeedSizeY = sectors_y; |
||
302 | qDebug() << "setting up of external seed maps finished"; |
||
303 | } |
||
304 | |||
305 | void SeedDispersal::finalizeExternalSeeds() |
||
306 | { |
||
307 | if (mExternalSeedBaseMap) |
||
308 | delete mExternalSeedBaseMap; |
||
309 | mExternalSeedBaseMap = 0; |
||
310 | } |
||
311 | |||
391 | werner | 312 | // ************ Kernel ************** |
1102 | werner | 313 | void SeedDispersal::createKernel(Grid<float> &kernel, const double max_seed) |
391 | werner | 314 | { |
415 | werner | 315 | |
391 | werner | 316 | double max_dist = treemig_distanceTo(0.0001); |
317 | double cell_size = mSeedMap.cellsize(); |
||
415 | werner | 318 | int max_radius = int(max_dist / cell_size); |
391 | werner | 319 | // e.g.: cell_size: regeneration grid (e.g. 400qm), px-size: light-grid (4qm) |
445 | werner | 320 | double occupation = cell_size*cell_size / (cPxSize*cPxSize * mTM_occupancy); |
391 | werner | 321 | |
415 | werner | 322 | kernel.clear(); |
391 | werner | 323 | |
415 | werner | 324 | kernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
325 | int kernel_offset = max_radius; |
||
326 | |||
391 | werner | 327 | // filling of the kernel.... use the treemig |
415 | werner | 328 | QPoint center = QPoint(kernel_offset, kernel_offset); |
329 | const float *sk_end = kernel.end(); |
||
330 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
||
331 | double d = kernel.distance(center, kernel.indexOf(p)); |
||
1102 | werner | 332 | *p = d<=max_dist?static_cast<float>(treemig(d)):0.f; |
391 | werner | 333 | } |
334 | |||
335 | // normalize |
||
1102 | werner | 336 | float sum = kernel.sum(); |
391 | werner | 337 | if (sum==0. || occupation==0.) |
338 | throw IException("create seed kernel: sum of probabilities = 0!"); |
||
339 | |||
340 | // the sum of all kernel cells has to equal 1 |
||
1102 | werner | 341 | kernel.multiply(1.f/sum); |
391 | werner | 342 | // probabilities are derived in multiplying by seed number, and dividing by occupancy criterion |
1102 | werner | 343 | kernel.multiply( static_cast<float>( max_seed / occupation)); |
391 | werner | 344 | // all cells that get more seeds than the occupancy criterion are considered to have no seed limitation for regeneration |
415 | werner | 345 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
391 | werner | 346 | *p = qMin(*p, 1.f); |
347 | } |
||
348 | // set the parent cell to 1 |
||
415 | werner | 349 | kernel.valueAtIndex(kernel_offset, kernel_offset)=1.f; |
350 | |||
351 | |||
391 | werner | 352 | // some final statistics.... |
550 | werner | 353 | if (logLevelInfo()) qDebug() << "kernel setup.Species:"<< mSpecies->id() << "kernel-size: " << kernel.sizeX() << "x" << kernel.sizeY() << "pixels, sum: " << kernel.sum(); |
391 | werner | 354 | } |
355 | |||
356 | /* R-Code: |
||
930 | werner | 357 | treemig=function(as1,as2,ks,d) # two-part exponential function, cf. Lischke & Loeffler (2006), Annex |
391 | werner | 358 | { |
359 | p1=(1-ks)*exp(-d/as1)/as1 |
||
360 | if(as2>0){p2=ks*exp(-d/as2)/as2}else{p2=0} |
||
361 | p1+p2 |
||
362 | } |
||
363 | */ |
||
670 | werner | 364 | |
365 | /// the used kernel function |
||
366 | /// see also Appendix B of iland paper II (note the different variable names) |
||
391 | werner | 367 | double SeedDispersal::treemig(const double &distance) |
368 | { |
||
369 | double p1 = (1.-mTM_ks)*exp(-distance/mTM_as1)/mTM_as1; |
||
370 | double p2 = 0.; |
||
371 | if (mTM_as2>0.) |
||
372 | p2 = mTM_ks*exp(-distance/mTM_as2)/mTM_as2; |
||
373 | return p1 + p2; |
||
374 | } |
||
375 | |||
376 | /// calculate the distance where the probability falls below 'value' |
||
377 | double SeedDispersal::treemig_distanceTo(const double value) |
||
378 | { |
||
379 | double dist = 0.; |
||
380 | while (treemig(dist)>value && dist<10000.) |
||
381 | dist+=10; |
||
382 | return dist; |
||
383 | } |
||
384 | |||
764 | werner | 385 | void SeedDispersal::setupExternalSeedsForSpecies(Species *species) |
386 | { |
||
387 | if (!mExtSeedData.contains(species->id())) |
||
388 | return; // nothing to do |
||
389 | qDebug() << "setting up external seed map for" << species->id(); |
||
390 | QVector<double> &pcts = mExtSeedData[species->id()]; |
||
391 | mExternalSeedMap.setup(mSeedMap); |
||
392 | mExternalSeedMap.initialize(0.f); |
||
393 | for (int sector_x=0; sector_x<mExtSeedSizeX; ++sector_x) |
||
394 | for (int sector_y=0; sector_y<mExtSeedSizeY; ++sector_y) { |
||
395 | int xmin,xmax,ymin,ymax; |
||
396 | int fx = mExternalSeedMap.sizeX() / mExtSeedSizeX; // number of cells per sector |
||
397 | xmin = sector_x*fx; |
||
398 | xmax = (sector_x+1)*fx; |
||
399 | fx = mExternalSeedMap.sizeY() / mExtSeedSizeY; // number of cells per sector |
||
400 | ymin = sector_y*fx; |
||
401 | ymax = (sector_y+1)*fx; |
||
402 | // now loop over the whole sector |
||
403 | int index = sector_y*mExtSeedSizeX + sector_x; |
||
404 | double p = pcts[index]; |
||
405 | for (int y=ymin;y<ymax;++y) |
||
406 | for (int x=xmin;x<xmax;++x) { |
||
407 | // check |
||
408 | if (mExternalSeedBaseMap->valueAtIndex(x,y)==2.f) |
||
409 | if (drandom()<p) |
||
410 | mExternalSeedMap.valueAtIndex(x,y) = 1.f; // flag |
||
411 | } |
||
391 | werner | 412 | |
764 | werner | 413 | } |
414 | } |
||
415 | |||
416 | |||
391 | werner | 417 | // ************ Dispersal ************** |
418 | |||
419 | |||
375 | werner | 420 | /// debug function: loads a image of arbirtrary size... |
421 | void SeedDispersal::loadFromImage(const QString &fileName) |
||
422 | { |
||
423 | mSeedMap.clear(); |
||
424 | loadGridFromImage(fileName, mSeedMap); |
||
425 | for (float* p=mSeedMap.begin();p!=mSeedMap.end();++p) |
||
426 | *p = *p>0.8?1.f:0.f; |
||
427 | |||
428 | } |
||
429 | |||
430 | void SeedDispersal::clear() |
||
431 | { |
||
764 | werner | 432 | if (!mExternalSeedMap.isEmpty()) { |
433 | // we have a preprocessed initial value for the external seed map (see setupExternalSeeds() et al) |
||
434 | mSeedMap.copy(mExternalSeedMap); |
||
435 | return; |
||
436 | } |
||
437 | // clear the map |
||
1102 | werner | 438 | float background_value = static_cast<float>(mExternalSeedBackgroundInput); // there is potentitally a background probability <>0 for all pixels. |
836 | werner | 439 | mSeedMap.initialize(background_value); |
472 | werner | 440 | if (mHasExternalSeedInput) { |
441 | // if external seed input is enabled, the buffer area of the seed maps is |
||
442 | // "turned on", i.e. set to 1. |
||
1102 | werner | 443 | int buf_size = GlobalSettings::instance()->settings().valueInt("model.world.buffer",0.) / mSeedMap.cellsize(); |
491 | werner | 444 | // if a special buffer is defined, reduce the size of the input |
445 | if (mExternalSeedBuffer>0) |
||
446 | buf_size -= mExternalSeedBuffer; |
||
472 | werner | 447 | if (buf_size>0) { |
448 | int ix,iy; |
||
449 | for (iy=0;iy<mSeedMap.sizeY();++iy) |
||
450 | for (ix=0;ix<mSeedMap.sizeX(); ++ix) |
||
491 | werner | 451 | if (iy<buf_size || iy>=mSeedMap.sizeY()-buf_size || ix<buf_size || ix>=mSeedMap.sizeX()-buf_size) { |
452 | if (mExternalSeedDirection==0) { |
||
453 | // seeds from all directions |
||
454 | mSeedMap.valueAtIndex(ix,iy)=1.f; |
||
455 | } else { |
||
456 | // seeds only from specific directions |
||
457 | float value = 0.f; |
||
458 | if (isBitSet(mExternalSeedDirection,1) && ix>=mSeedMap.sizeX()-buf_size) value = 1; // north |
||
459 | if (isBitSet(mExternalSeedDirection,2) && iy<buf_size) value = 1; // east |
||
460 | if (isBitSet(mExternalSeedDirection,3) && ix<buf_size) value = 1; // south |
||
461 | if (isBitSet(mExternalSeedDirection,4) && iy>=mSeedMap.sizeY()-buf_size) value = 1; // west |
||
462 | mSeedMap.valueAtIndex(ix,iy)=value; |
||
463 | } |
||
464 | } |
||
472 | werner | 465 | } else { |
466 | qDebug() << "external seed input: Error: invalid buffer size???"; |
||
467 | } |
||
468 | } |
||
375 | werner | 469 | } |
470 | |||
471 | void SeedDispersal::execute() |
||
472 | { |
||
991 | werner | 473 | #ifdef ILAND_GUI |
474 | int year = GlobalSettings::instance()->currentYear(); |
||
475 | QString path; |
||
472 | werner | 476 | if (mDumpSeedMaps) { |
1102 | werner | 477 | path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 478 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_before_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
618 | werner | 479 | qDebug() << "saved seed map image to" << path; |
992 | werner | 480 | } |
989 | werner | 481 | #else |
992 | werner | 482 | if (mDumpSeedMaps) |
989 | werner | 483 | qDebug() << "saving of seedmaps only supported in the iLand GUI."; |
484 | #endif |
||
992 | werner | 485 | |
1102 | werner | 486 | DebugTimer t("seed dispersal", true); |
391 | werner | 487 | { |
375 | werner | 488 | // (1) detect edges |
619 | werner | 489 | if (edgeDetection()) { |
1106 | werner | 490 | |
491 | #ifdef ILAND_GUI |
||
492 | if (mDumpSeedMaps) { |
||
493 | gridToImage(seedMap(), true, -1., 1.).save(QString("%1/seed_edge_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
494 | } |
||
495 | #endif |
||
496 | |||
497 | // (2) distribute seed probabilites from edges |
||
619 | werner | 498 | distribute(); |
391 | werner | 499 | } |
619 | werner | 500 | } |
989 | werner | 501 | #ifdef ILAND_GUI |
1102 | werner | 502 | if (mDumpSeedMaps) { |
503 | qDebug() << "finished seed dispersal for species" << mSpecies->id() << t.elapsed(); |
||
472 | werner | 504 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_after_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
1102 | werner | 505 | } |
1064 | werner | 506 | |
507 | if (!mDumpNextYearFileName.isEmpty()) { |
||
508 | Helper::saveToTextFile(GlobalSettings::instance()->path(mDumpNextYearFileName), gridToESRIRaster(seedMap())); |
||
509 | qDebug() << "saved seed map for " << species()->id() << "to" << GlobalSettings::instance()->path(mDumpNextYearFileName); |
||
510 | mDumpNextYearFileName = QString(); |
||
511 | } |
||
512 | |||
989 | werner | 513 | #endif |
375 | werner | 514 | } |
515 | |||
373 | werner | 516 | /** scans the seed image and detects "edges". |
517 | edges are then subsequently marked (set to -1). This is pass 1 of the seed distribution process. |
||
518 | */ |
||
619 | werner | 519 | bool SeedDispersal::edgeDetection() |
373 | werner | 520 | { |
521 | float *p_above, *p, *p_below; |
||
522 | int dy = mSeedMap.sizeY(); |
||
523 | int dx = mSeedMap.sizeX(); |
||
524 | int x,y; |
||
619 | werner | 525 | bool found = false; |
1106 | werner | 526 | |
527 | // fill mini-gaps |
||
373 | werner | 528 | for (y=1;y<dy-1;++y){ |
529 | p = mSeedMap.ptr(1,y); |
||
530 | p_above = p - dx; // one line above |
||
531 | p_below = p + dx; // one line below |
||
532 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
1106 | werner | 533 | if (*p < 0.999f) { |
534 | |||
535 | if ((*(p_above-1)==1.f) + (*p_above==1.f) + (*(p_above+1)==1.f) + |
||
536 | (*(p-1)==1.f) + (*(p+1)==1.f) + |
||
537 | (*(p_below-1)==1.f) + (*p_below==1.f) + (*(p_below+1)==1.f) > 3) |
||
538 | *p=0.999f; // if more than 3 neighbors are active pixels, the value is high |
||
539 | } |
||
540 | |||
541 | } |
||
542 | } |
||
543 | |||
544 | |||
545 | // now detect the edges |
||
546 | int n_edges=0 ; |
||
547 | for (y=1;y<dy-1;++y){ |
||
548 | p = mSeedMap.ptr(1,y); |
||
549 | p_above = p - dx; // one line above |
||
550 | p_below = p + dx; // one line below |
||
551 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
552 | if (*p == 1.f) { |
||
619 | werner | 553 | found = true; |
1106 | werner | 554 | if ( (*(p_above-1)<0.999f && *(p_above-1)>=0.f) |
555 | || (*p_above<0.999f && *p_above>=0.f) |
||
556 | || (*(p_above+1)<0.999f && *(p_above+1)>=0.f) |
||
557 | || (*(p-1)<0.999f && *(p-1)>=0.f) |
||
558 | || (*(p+1)<0.999f && (*p+1)>=0.f) |
||
559 | || (*(p_below-1)<0.999f && *(p_below-1)>=0.f) |
||
560 | || (*p_below<0.999f && *p_below>=0.f) |
||
561 | || (*(p_below+1)<0.999f && *(p_below+1)>=0.f ) ) { |
||
562 | *p=-1.f; // if any surrounding pixel is >=0 & <0.999: -> mark as edge |
||
563 | ++n_edges; |
||
564 | } |
||
373 | werner | 565 | } |
566 | |||
567 | } |
||
568 | } |
||
1106 | werner | 569 | //qDebug() << "species:" << mSpecies->id() << "# of edge-pixels:" << n_edges; |
619 | werner | 570 | return found; |
373 | werner | 571 | } |
572 | |||
573 | /** do the seed probability distribution. |
||
574 | This is phase 2. Apply the seed kernel for each "edge" point identified in phase 1. |
||
575 | */ |
||
576 | void SeedDispersal::distribute() |
||
577 | { |
||
578 | int x,y; |
||
579 | |||
580 | float *end = mSeedMap.end(); |
||
581 | float *p = mSeedMap.begin(); |
||
415 | werner | 582 | // choose the kernel depending whether there is a seed year for the current species or not |
583 | Grid<float> &kernel = species()->isSeedYear()?mKernelSeedYear:mKernelNonSeedYear; |
||
584 | int offset = kernel.sizeX() / 2; // offset is the index of the center pixel |
||
373 | werner | 585 | for(;p!=end;++p) { |
375 | werner | 586 | if (*p==-1.f) { |
373 | werner | 587 | // edge pixel found. Now apply the kernel.... |
588 | QPoint pt=mSeedMap.indexOf(p); |
||
415 | werner | 589 | for (y=-offset;y<=offset;++y) |
590 | for (x=-offset;x<=offset;++x) |
||
373 | werner | 591 | if (mSeedMap.isIndexValid(pt.x()+x, pt.y()+y)) { |
592 | float &val = mSeedMap.valueAtIndex(pt.x()+x, pt.y()+y); |
||
1106 | werner | 593 | if (val!=-1.f) |
415 | werner | 594 | val = qMin(1.f - (1.f - val)*(1.f-kernel.valueAtIndex(x+offset, y+offset)),1.f ); |
373 | werner | 595 | } |
375 | werner | 596 | *p=1.f; // mark as processed |
373 | werner | 597 | } // *p==1 |
598 | } // for() |
||
599 | } |