Rev 1217 | Rev 1220 | 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; |
1180 | werner | 64 | mProbMode = false; |
391 | werner | 65 | |
66 | const float seedmap_size = 20.f; |
||
373 | werner | 67 | // setup of seed map |
68 | mSeedMap.clear(); |
||
391 | werner | 69 | mSeedMap.setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
373 | werner | 70 | mSeedMap.initialize(0.); |
1180 | werner | 71 | if (!mProbMode) { |
72 | mSourceMap.setup(mSeedMap); |
||
73 | mSourceMap.initialize(0.); |
||
74 | } |
||
764 | werner | 75 | mExternalSeedMap.clear(); |
391 | werner | 76 | mIndexFactor = int(seedmap_size) / cPxSize; // ratio seed grid / lip-grid: |
550 | werner | 77 | if (logLevelInfo()) qDebug() << "Seed map setup. Species:"<< mSpecies->id() << "kernel-size: " << mSeedMap.sizeX() << "x" << mSeedMap.sizeY() << "pixels."; |
373 | werner | 78 | |
445 | werner | 79 | if (mSpecies==0) |
80 | throw IException("Setup of SeedDispersal: Species not defined."); |
||
81 | |||
802 | werner | 82 | if (fmod(GlobalSettings::instance()->settings().valueDouble("model.world.buffer",0),seedmap_size) != 0.) |
83 | 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,...))."); |
||
84 | |||
415 | werner | 85 | // settings |
445 | werner | 86 | mTM_occupancy = 1.; // is currently constant |
1176 | werner | 87 | // copy values for the species parameters: |
88 | mSpecies->treeMigKernel(mTM_as1, mTM_as2, mTM_ks); |
||
445 | werner | 89 | mTM_fecundity_cell = mSpecies->fecundity_m2() * seedmap_size*seedmap_size * mTM_occupancy; // scale to production for the whole cell |
90 | mNonSeedYearFraction = mSpecies->nonSeedYearFraction(); |
||
1176 | werner | 91 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.settings.seedDispersal")); |
92 | mKernelThresholdArea = xml.valueDouble(".longDistanceDispersal.thresholdArea", 0.0001); |
||
93 | mKernelThresholdLDD = xml.valueDouble(".longDistanceDispersal.thresholdLDD", 0.0001); |
||
1180 | werner | 94 | mLDDSeedlings = xml.valueDouble(".longDistanceDispersal.LDDSeedlings", 0.0001); |
1176 | werner | 95 | mLDDRings = xml.valueInt(".longDistanceDispersal.rings", 4); |
415 | werner | 96 | |
1180 | werner | 97 | mLDDSeedlings = qMax(mLDDSeedlings, static_cast<float>(mKernelThresholdArea)); |
415 | werner | 98 | |
1180 | werner | 99 | // long distance dispersal |
100 | double ldd_area = setupLDD(); |
||
1176 | werner | 101 | |
1180 | werner | 102 | createKernel(mKernelSeedYear, mTM_fecundity_cell, 1. - ldd_area); |
1176 | werner | 103 | |
415 | werner | 104 | // the kernel for non seed years looks similar, but is simply linearly scaled down |
105 | // using the species parameter NonSeedYearFraction. |
||
106 | // the central pixel still gets the value of 1 (i.e. 100% probability) |
||
1180 | werner | 107 | createKernel(mKernelNonSeedYear, mTM_fecundity_cell*mNonSeedYearFraction, 1. - ldd_area); |
415 | werner | 108 | |
1167 | werner | 109 | if (mSpecies->fecunditySerotiny()>0.) { |
110 | // an extra seed map is used for storing information related to post-fire seed rain |
||
1168 | werner | 111 | mSeedMapSerotiny.clear(); |
112 | mSeedMapSerotiny.setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
||
113 | mSeedMapSerotiny.initialize(0.); |
||
1167 | werner | 114 | |
115 | // set up the special seed kernel for post fire seed rain |
||
1180 | werner | 116 | createKernel(mKernelSerotiny, mTM_fecundity_cell * mSpecies->fecunditySerotiny(),1.); |
1167 | werner | 117 | qDebug() << "created extra seed map and serotiny seed kernel for species" << mSpecies->name() << "with fecundity factor" << mSpecies->fecunditySerotiny(); |
118 | } |
||
119 | mHasPendingSerotiny = false; |
||
120 | |||
472 | werner | 121 | // debug info |
122 | mDumpSeedMaps = GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false); |
||
123 | if (mDumpSeedMaps) { |
||
1102 | werner | 124 | QString path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 125 | Helper::saveToTextFile(QString("%1/seedkernelYes_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelSeedYear)); |
126 | Helper::saveToTextFile(QString("%1/seedkernelNo_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelNonSeedYear)); |
||
1168 | werner | 127 | if (!mKernelSerotiny.isEmpty()) |
128 | Helper::saveToTextFile(QString("%1/seedkernelSerotiny_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelSerotiny)); |
||
417 | werner | 129 | } |
1176 | werner | 130 | |
1180 | werner | 131 | |
472 | werner | 132 | // external seeds |
481 | werner | 133 | mHasExternalSeedInput = false; |
491 | werner | 134 | mExternalSeedBuffer = 0; |
135 | mExternalSeedDirection = 0; |
||
836 | werner | 136 | mExternalSeedBackgroundInput = 0.; |
472 | werner | 137 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.externalSeedEnabled",false)) { |
764 | werner | 138 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) { |
139 | // external seed input specified by sectors and around the project area (seedbelt) |
||
140 | setupExternalSeedsForSpecies(mSpecies); |
||
141 | } else { |
||
142 | // external seeds specified fixedly per cardinal direction |
||
143 | // current species in list?? |
||
144 | mHasExternalSeedInput = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSpecies").contains(mSpecies->id()); |
||
145 | QString dir = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSource").toLower(); |
||
146 | // encode cardinal positions as bits: e.g: "e,w" -> 6 |
||
147 | mExternalSeedDirection += dir.contains("n")?1:0; |
||
148 | mExternalSeedDirection += dir.contains("e")?2:0; |
||
149 | mExternalSeedDirection += dir.contains("s")?4:0; |
||
150 | mExternalSeedDirection += dir.contains("w")?8:0; |
||
837 | werner | 151 | QStringList buffer_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBuffer").split(QRegExp("([^\\.\\w]+)")); |
764 | werner | 152 | int index = buffer_list.indexOf(mSpecies->id()); |
153 | if (index>=0) { |
||
154 | mExternalSeedBuffer = buffer_list[index+1].toInt(); |
||
155 | qDebug() << "enabled special buffer for species" <<mSpecies->id() << ": distance of" << mExternalSeedBuffer << "pixels = " << mExternalSeedBuffer*20. << "m"; |
||
156 | } |
||
836 | werner | 157 | |
158 | // background seed rain (i.e. for the full landscape), use regexp |
||
837 | werner | 159 | QStringList background_input_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBackgroundInput").split(QRegExp("([^\\.\\w]+)")); |
836 | werner | 160 | index = background_input_list.indexOf(mSpecies->id()); |
161 | if (index>=0) { |
||
162 | mExternalSeedBackgroundInput = background_input_list[index+1].toDouble(); |
||
163 | qDebug() << "enabled background seed input (for full area) for species" <<mSpecies->id() << ": p=" << mExternalSeedBackgroundInput; |
||
164 | } |
||
165 | |||
764 | werner | 166 | if (mHasExternalSeedInput) |
167 | qDebug() << "External seed input enabled for" << mSpecies->id(); |
||
491 | werner | 168 | } |
472 | werner | 169 | } |
415 | werner | 170 | |
373 | werner | 171 | // setup of seed kernel |
391 | werner | 172 | // const int max_radius = 15; // pixels |
173 | // |
||
174 | // mSeedKernel.clear(); |
||
175 | // mSeedKernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
||
176 | // mKernelOffset = max_radius; |
||
177 | // // filling of the kernel.... for simplicity: a linear kernel |
||
178 | // QPoint center = QPoint(mKernelOffset, mKernelOffset); |
||
179 | // const double max_dist = max_radius * seedmap_size; |
||
180 | // for (float *p=mSeedKernel.begin(); p!=mSeedKernel.end();++p) { |
||
181 | // double d = mSeedKernel.distance(center, mSeedKernel.indexOf(p)); |
||
182 | // *p = qMax( 1. - d / max_dist, 0.); |
||
183 | // } |
||
373 | werner | 184 | |
185 | |||
186 | // randomize seed map.... set 1/3 to "filled" |
||
375 | werner | 187 | //for (int i=0;i<mSeedMap.count(); i++) |
188 | // mSeedMap.valueAtIndex(mSeedMap.randomPosition()) = 1.; |
||
373 | werner | 189 | |
190 | |||
375 | werner | 191 | // QImage img = gridToImage(mSeedMap, true, -1., 1.); |
192 | // img.save("seedmap.png"); |
||
193 | |||
194 | // img = gridToImage(mSeedMap, true, -1., 1.); |
||
764 | werner | 195 | // img.save("seedmap_e.png"); |
373 | werner | 196 | } |
197 | |||
764 | werner | 198 | void SeedDispersal::setupExternalSeeds() |
199 | { |
||
200 | mExternalSeedBaseMap = 0; |
||
201 | if (!GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) |
||
202 | return; |
||
203 | |||
978 | werner | 204 | DebugTimer t("setup of external seed maps."); |
764 | werner | 205 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.settings.seedDispersal.seedBelt")); |
1102 | werner | 206 | int seedbelt_width =xml.valueInt(".width",10); |
764 | werner | 207 | // setup of sectors |
208 | // setup of base map |
||
209 | const float seedmap_size = 20.f; |
||
210 | mExternalSeedBaseMap = new Grid<float>; |
||
211 | mExternalSeedBaseMap->setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
||
212 | mExternalSeedBaseMap->initialize(0.); |
||
213 | if (mExternalSeedBaseMap->count()*4 != GlobalSettings::instance()->model()->heightGrid()->count()) |
||
947 | werner | 214 | 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 | 215 | // make a copy of the 10m height grid in lower resolution and mark pixels that are forested and outside of |
216 | // the project area. |
||
217 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) |
||
218 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
765 | werner | 219 | bool val = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isForestOutside(); |
764 | werner | 220 | mExternalSeedBaseMap->valueAtIndex(x,y) = val?1.f:0.f; |
221 | if(GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isValid()) |
||
222 | mExternalSeedBaseMap->valueAtIndex(x,y) = -1.f; |
||
223 | } |
||
836 | werner | 224 | QString path = GlobalSettings::instance()->path(GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath")); |
765 | werner | 225 | |
226 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
||
989 | werner | 227 | #ifdef ILAND_GUI |
765 | werner | 228 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
229 | img.save(path + "/seedbeltmap_before.png"); |
||
989 | werner | 230 | #endif |
765 | werner | 231 | } |
764 | werner | 232 | // img.save("seedmap.png"); |
233 | // now scan the pixels of the belt: paint all pixels that are close to the project area |
||
234 | // we do this 4 times (for all cardinal direcitons) |
||
235 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
236 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
237 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
238 | continue; |
||
765 | werner | 239 | int look_forward = std::min(x + seedbelt_width, mExternalSeedBaseMap->sizeX()-1); |
1106 | werner | 240 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 241 | // fill pixels |
242 | for(; x<look_forward;++x) { |
||
243 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
244 | if (v==1.f) v=2.f; |
||
245 | } |
||
246 | } |
||
247 | } |
||
248 | } |
||
249 | // right to left |
||
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(x - seedbelt_width, 0); |
||
1106 | werner | 255 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 256 | // fill pixels |
257 | for(; x>look_forward;--x) { |
||
258 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
259 | if (v==1.f) v=2.f; |
||
260 | } |
||
261 | } |
||
262 | } |
||
263 | } |
||
264 | // up and down *** |
||
265 | // from top to bottom |
||
266 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
267 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
268 | |||
269 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
270 | continue; |
||
765 | werner | 271 | int look_forward = std::min(y + seedbelt_width, mExternalSeedBaseMap->sizeY()-1); |
764 | werner | 272 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
273 | // fill pixels |
||
274 | for(; y<look_forward;++y) { |
||
275 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
276 | if (v==1.f) v=2.f; |
||
277 | } |
||
278 | } |
||
279 | } |
||
280 | } |
||
281 | // bottom to top *** |
||
282 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
283 | for (int x=mExternalSeedBaseMap->sizeX();x>=0;--x) { |
||
284 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
285 | continue; |
||
286 | int look_forward = std::max(y - seedbelt_width, 0); |
||
287 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
||
288 | // fill pixels |
||
289 | for(; y>look_forward;--y) { |
||
290 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
291 | if (v==1.f) v=2.f; |
||
292 | } |
||
293 | } |
||
294 | } |
||
295 | } |
||
1180 | werner | 296 | |
765 | werner | 297 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
989 | werner | 298 | #ifdef ILAND_GUI |
765 | werner | 299 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
972 | werner | 300 | img.save(path + "/seedbeltmap_after.png"); |
989 | werner | 301 | #endif |
765 | werner | 302 | } |
764 | werner | 303 | mExtSeedData.clear(); |
1102 | werner | 304 | int sectors_x = xml.valueInt("sizeX",0); |
305 | int sectors_y = xml.valueInt("sizeY",0); |
||
764 | werner | 306 | if(sectors_x<1 || sectors_y<1) |
307 | throw IException(QString("setup of external seed dispersal: invalid number of sectors x=%1 y=%3").arg(sectors_x).arg(sectors_y)); |
||
308 | QDomElement elem = xml.node("."); |
||
309 | for(QDomNode n = elem.firstChild(); !n.isNull(); n = n.nextSibling()) { |
||
310 | if (n.nodeName().startsWith("species")) { |
||
311 | QStringList coords = n.nodeName().split("_"); |
||
312 | if (coords.count()!=3) |
||
313 | throw IException("external seed species definition is not valid: " + n.nodeName()); |
||
314 | int x = coords[1].toInt(); |
||
315 | int y = coords[2].toInt(); |
||
316 | if (x<0 || x>=sectors_x || y<0 || y>=sectors_y) |
||
317 | throw IException(QString("invalid sector for specifiing external seed input (x y): %1 %2 ").arg(x).arg(y) ); |
||
318 | int index = y*sectors_x + x; |
||
319 | |||
320 | QString text = xml.value("." + n.nodeName()); |
||
321 | qDebug() << "processing element " << n.nodeName() << "x,y:" << x << y << text; |
||
322 | // we assume pairs of name and fraction |
||
323 | QStringList species_list = text.split(" "); |
||
324 | for (int i=0;i<species_list.count();++i) { |
||
325 | QVector<double> &space = mExtSeedData[species_list[i]]; |
||
326 | if (space.isEmpty()) |
||
327 | space.resize(sectors_x*sectors_y); // are initialized to 0s |
||
328 | double fraction = species_list[++i].toDouble(); |
||
329 | space[index] = fraction; |
||
330 | } |
||
331 | } |
||
332 | } |
||
333 | mExtSeedSizeX = sectors_x; |
||
334 | mExtSeedSizeY = sectors_y; |
||
335 | qDebug() << "setting up of external seed maps finished"; |
||
336 | } |
||
337 | |||
338 | void SeedDispersal::finalizeExternalSeeds() |
||
339 | { |
||
340 | if (mExternalSeedBaseMap) |
||
341 | delete mExternalSeedBaseMap; |
||
342 | mExternalSeedBaseMap = 0; |
||
343 | } |
||
344 | |||
1167 | werner | 345 | void SeedDispersal::seedProductionSerotiny(const QPoint &position_index) |
346 | { |
||
1168 | werner | 347 | if (mSeedMapSerotiny.isEmpty()) |
1167 | werner | 348 | throw IException("Invalid use seedProductionSerotiny(): tried to set a seed source for a non-serotinous species!"); |
1168 | werner | 349 | mSeedMapSerotiny.valueAtIndex(position_index.x()/mIndexFactor, position_index.y()/mIndexFactor)=1.f; |
1167 | werner | 350 | mHasPendingSerotiny = true; |
351 | } |
||
352 | |||
391 | werner | 353 | // ************ Kernel ************** |
1180 | werner | 354 | void SeedDispersal::createKernel(Grid<float> &kernel, const double max_seed, const double scale_area) |
391 | werner | 355 | { |
415 | werner | 356 | |
1180 | werner | 357 | double max_dist = treemig_distanceTo(mKernelThresholdArea / species()->fecundity_m2()); |
391 | werner | 358 | double cell_size = mSeedMap.cellsize(); |
415 | werner | 359 | int max_radius = int(max_dist / cell_size); |
391 | werner | 360 | // e.g.: cell_size: regeneration grid (e.g. 400qm), px-size: light-grid (4qm) |
445 | werner | 361 | double occupation = cell_size*cell_size / (cPxSize*cPxSize * mTM_occupancy); |
391 | werner | 362 | |
415 | werner | 363 | kernel.clear(); |
391 | werner | 364 | |
415 | werner | 365 | kernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
366 | int kernel_offset = max_radius; |
||
367 | |||
1180 | werner | 368 | // filling of the kernel.... use the treemig density function |
1179 | werner | 369 | double dist_center_cell = sqrt(cell_size*cell_size/M_PI); |
415 | werner | 370 | QPoint center = QPoint(kernel_offset, kernel_offset); |
371 | const float *sk_end = kernel.end(); |
||
372 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
||
373 | double d = kernel.distance(center, kernel.indexOf(p)); |
||
1178 | werner | 374 | if (d==0.) |
1179 | werner | 375 | *p = treemig_centercell(dist_center_cell); // r is the radius of a circle with the same area as a cell |
1178 | werner | 376 | else |
1180 | werner | 377 | *p = d<=max_dist?static_cast<float>(( treemig(d+dist_center_cell) + treemig(d-dist_center_cell))/2.f * cell_size*cell_size ):0.f; |
391 | werner | 378 | } |
379 | |||
380 | // normalize |
||
1102 | werner | 381 | float sum = kernel.sum(); |
391 | werner | 382 | if (sum==0. || occupation==0.) |
383 | throw IException("create seed kernel: sum of probabilities = 0!"); |
||
384 | |||
1180 | werner | 385 | // the sum of all kernel cells has to equal 1 (- long distance dispersal) |
386 | kernel.multiply(scale_area/sum); |
||
1178 | werner | 387 | |
1180 | werner | 388 | |
389 | if (mProbMode) { |
||
390 | // probabilities are derived in multiplying by seed number, and dividing by occupancy criterion |
||
391 | float fecundity_factor = static_cast<float>( max_seed / occupation); |
||
392 | kernel.multiply( fecundity_factor ); |
||
393 | // all cells that get more seeds than the occupancy criterion are considered to have no seed limitation for regeneration |
||
394 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
||
395 | *p = qMin(*p, 1.f); |
||
396 | } |
||
391 | werner | 397 | } |
398 | // set the parent cell to 1 |
||
1178 | werner | 399 | //kernel.valueAtIndex(kernel_offset, kernel_offset)=1.f; |
415 | werner | 400 | |
401 | |||
391 | werner | 402 | // some final statistics.... |
1180 | werner | 403 | if (logLevelInfo()) |
404 | qDebug() << "kernel setup. Species:"<< mSpecies->id() << "kernel-size: " << kernel.sizeX() << "x" << kernel.sizeY() << "pixels, sum (after scaling): " << kernel.sum(); |
||
1176 | werner | 405 | |
1180 | werner | 406 | |
391 | werner | 407 | } |
408 | |||
1180 | werner | 409 | double SeedDispersal::setupLDD() |
1176 | werner | 410 | { |
411 | mLDDDensity.clear(); mLDDDistance.clear(); |
||
412 | if (mKernelThresholdLDD >= mKernelThresholdArea) { |
||
413 | // no long distance dispersal |
||
1180 | werner | 414 | return 0.; |
1176 | werner | 415 | |
416 | } |
||
1180 | werner | 417 | double r_min = treemig_distanceTo(mKernelThresholdArea / species()->fecundity_m2()); |
418 | double r_max = treemig_distanceTo(mKernelThresholdLDD / species()->fecundity_m2()); |
||
1176 | werner | 419 | |
420 | |||
421 | mLDDDistance.push_back(r_min); |
||
1180 | werner | 422 | double ldd_sum = 0.; |
1176 | werner | 423 | for (int i=0;i<mLDDRings;++i) { |
424 | double r_in = mLDDDistance.last(); |
||
425 | mLDDDistance.push_back(mLDDDistance.last() + (r_max-r_min)/static_cast<float>(mLDDRings)); |
||
426 | double r_out = mLDDDistance.last(); |
||
427 | // calculate the value of the kernel for the middle of the ring |
||
428 | double ring_in = treemig(r_in); // kernel value at the inner border of the ring |
||
429 | double ring_out = treemig(r_out); // kernel value at the outer border of the ring |
||
1180 | werner | 430 | double ring_val = ring_in*0.4 + ring_out*0.6; // this is the average p -- 0.4/0.6 better estimate the nonlinear behavior (fits very well for medium to large kernels, e.g. piab) |
431 | // |
||
1176 | werner | 432 | // calculate the area of the ring |
433 | double ring_area = (r_out*r_out - r_in*r_in)*M_PI; // in square meters |
||
1180 | werner | 434 | // the number of px considers the fecundity |
435 | double n_px = ring_val * ring_area * species()->fecundity_m2() / mLDDSeedlings; |
||
436 | ldd_sum += ring_val * ring_area; // this fraction of the full kernel (=1) is distributed in theis ring |
||
1176 | werner | 437 | |
438 | mLDDDensity.push_back(n_px); |
||
439 | } |
||
1204 | werner | 440 | if (logLevelInfo()) |
441 | qDebug() << "Setup LDD for" << species()->name() << ", using probability: "<< mLDDSeedlings<< ": Distances:" << mLDDDistance << ", seed pixels:" << mLDDDensity << "covered prob:" << ldd_sum; |
||
1176 | werner | 442 | |
1180 | werner | 443 | return ldd_sum; |
1176 | werner | 444 | } |
445 | |||
391 | werner | 446 | /* R-Code: |
930 | werner | 447 | treemig=function(as1,as2,ks,d) # two-part exponential function, cf. Lischke & Loeffler (2006), Annex |
391 | werner | 448 | { |
449 | p1=(1-ks)*exp(-d/as1)/as1 |
||
450 | if(as2>0){p2=ks*exp(-d/as2)/as2}else{p2=0} |
||
451 | p1+p2 |
||
452 | } |
||
453 | */ |
||
670 | werner | 454 | |
455 | /// the used kernel function |
||
456 | /// see also Appendix B of iland paper II (note the different variable names) |
||
1178 | werner | 457 | /// the function returns the seed density at a point with distance 'distance'. |
391 | werner | 458 | double SeedDispersal::treemig(const double &distance) |
459 | { |
||
460 | double p1 = (1.-mTM_ks)*exp(-distance/mTM_as1)/mTM_as1; |
||
461 | double p2 = 0.; |
||
462 | if (mTM_as2>0.) |
||
463 | p2 = mTM_ks*exp(-distance/mTM_as2)/mTM_as2; |
||
1178 | werner | 464 | double s = p1 + p2; |
465 | // 's' is the density for radius 'distance' - not for specific point with that distance. |
||
466 | // (i.e. the integral over the one-dimensional treemig function is 1, but if applied for 2d cells, the |
||
467 | // sum would be much larger as all seeds arriving at 'distance' would be arriving somewhere at the circle with radius 'distance') |
||
468 | // convert that to a density at a point, by dividing with the circumference at the circle with radius 'distance' |
||
469 | s = s / (2.*std::max(distance, 0.01)*M_PI); |
||
470 | |||
471 | return s; |
||
391 | werner | 472 | } |
473 | |||
1178 | werner | 474 | double SeedDispersal::treemig_centercell(const double &max_distance) |
475 | { |
||
476 | // use 100 steps and calculate dispersal kernel for consecutive rings |
||
477 | double sum = 0.; |
||
478 | for (int i=0;i<100;i++) { |
||
479 | double r_in = i*max_distance/100.; |
||
480 | double r_out = (i+1)*max_distance/100.; |
||
481 | double ring_area = (r_out*r_out-r_in*r_in)*M_PI; |
||
482 | // the value of each ring is: treemig(r) * area of the ring |
||
483 | sum += treemig((r_out+r_in)/2.)*ring_area; |
||
484 | } |
||
485 | return sum; |
||
486 | } |
||
487 | |||
391 | werner | 488 | /// calculate the distance where the probability falls below 'value' |
489 | double SeedDispersal::treemig_distanceTo(const double value) |
||
490 | { |
||
491 | double dist = 0.; |
||
492 | while (treemig(dist)>value && dist<10000.) |
||
493 | dist+=10; |
||
494 | return dist; |
||
495 | } |
||
496 | |||
764 | werner | 497 | void SeedDispersal::setupExternalSeedsForSpecies(Species *species) |
498 | { |
||
499 | if (!mExtSeedData.contains(species->id())) |
||
500 | return; // nothing to do |
||
501 | qDebug() << "setting up external seed map for" << species->id(); |
||
502 | QVector<double> &pcts = mExtSeedData[species->id()]; |
||
503 | mExternalSeedMap.setup(mSeedMap); |
||
504 | mExternalSeedMap.initialize(0.f); |
||
505 | for (int sector_x=0; sector_x<mExtSeedSizeX; ++sector_x) |
||
506 | for (int sector_y=0; sector_y<mExtSeedSizeY; ++sector_y) { |
||
507 | int xmin,xmax,ymin,ymax; |
||
508 | int fx = mExternalSeedMap.sizeX() / mExtSeedSizeX; // number of cells per sector |
||
509 | xmin = sector_x*fx; |
||
510 | xmax = (sector_x+1)*fx; |
||
511 | fx = mExternalSeedMap.sizeY() / mExtSeedSizeY; // number of cells per sector |
||
512 | ymin = sector_y*fx; |
||
513 | ymax = (sector_y+1)*fx; |
||
514 | // now loop over the whole sector |
||
515 | int index = sector_y*mExtSeedSizeX + sector_x; |
||
516 | double p = pcts[index]; |
||
517 | for (int y=ymin;y<ymax;++y) |
||
518 | for (int x=xmin;x<xmax;++x) { |
||
519 | // check |
||
520 | if (mExternalSeedBaseMap->valueAtIndex(x,y)==2.f) |
||
521 | if (drandom()<p) |
||
522 | mExternalSeedMap.valueAtIndex(x,y) = 1.f; // flag |
||
523 | } |
||
391 | werner | 524 | |
764 | werner | 525 | } |
1180 | werner | 526 | if (!mProbMode) { |
527 | // scale external seed values to have pixels with LAI=3 |
||
528 | for (float *p=mExternalSeedMap.begin(); p!=mExternalSeedMap.end(); ++p) |
||
529 | *p *= 3.f * mExternalSeedMap.cellsize()*mExternalSeedMap.cellsize(); |
||
530 | } |
||
764 | werner | 531 | } |
532 | |||
533 | |||
391 | werner | 534 | // ************ Dispersal ************** |
535 | |||
536 | |||
375 | werner | 537 | /// debug function: loads a image of arbirtrary size... |
538 | void SeedDispersal::loadFromImage(const QString &fileName) |
||
539 | { |
||
540 | mSeedMap.clear(); |
||
541 | loadGridFromImage(fileName, mSeedMap); |
||
542 | for (float* p=mSeedMap.begin();p!=mSeedMap.end();++p) |
||
543 | *p = *p>0.8?1.f:0.f; |
||
544 | |||
545 | } |
||
546 | |||
547 | void SeedDispersal::clear() |
||
548 | { |
||
1180 | werner | 549 | Grid<float> *seed_map = &mSeedMap; |
550 | if (!mProbMode) { |
||
551 | seed_map = &mSourceMap; |
||
552 | mSeedMap.initialize(0.f); |
||
553 | } |
||
764 | werner | 554 | if (!mExternalSeedMap.isEmpty()) { |
555 | // we have a preprocessed initial value for the external seed map (see setupExternalSeeds() et al) |
||
1180 | werner | 556 | seed_map->copy(mExternalSeedMap); |
764 | werner | 557 | return; |
558 | } |
||
559 | // clear the map |
||
1102 | werner | 560 | float background_value = static_cast<float>(mExternalSeedBackgroundInput); // there is potentitally a background probability <>0 for all pixels. |
1180 | werner | 561 | seed_map->initialize(background_value); |
472 | werner | 562 | if (mHasExternalSeedInput) { |
563 | // if external seed input is enabled, the buffer area of the seed maps is |
||
564 | // "turned on", i.e. set to 1. |
||
1180 | werner | 565 | int buf_size = GlobalSettings::instance()->settings().valueInt("model.world.buffer",0.) / static_cast<int>(seed_map->cellsize()); |
491 | werner | 566 | // if a special buffer is defined, reduce the size of the input |
567 | if (mExternalSeedBuffer>0) |
||
568 | buf_size -= mExternalSeedBuffer; |
||
472 | werner | 569 | if (buf_size>0) { |
570 | int ix,iy; |
||
1180 | werner | 571 | for (iy=0;iy<seed_map->sizeY();++iy) |
572 | for (ix=0;ix<seed_map->sizeX(); ++ix) |
||
573 | if (iy<buf_size || iy>=seed_map->sizeY()-buf_size || ix<buf_size || ix>=seed_map->sizeX()-buf_size) { |
||
491 | werner | 574 | if (mExternalSeedDirection==0) { |
575 | // seeds from all directions |
||
1180 | werner | 576 | seed_map->valueAtIndex(ix,iy)=1.f; |
491 | werner | 577 | } else { |
578 | // seeds only from specific directions |
||
579 | float value = 0.f; |
||
1180 | werner | 580 | if (isBitSet(mExternalSeedDirection,1) && ix>=seed_map->sizeX()-buf_size) value = 1; // north |
491 | werner | 581 | if (isBitSet(mExternalSeedDirection,2) && iy<buf_size) value = 1; // east |
582 | if (isBitSet(mExternalSeedDirection,3) && ix<buf_size) value = 1; // south |
||
1180 | werner | 583 | if (isBitSet(mExternalSeedDirection,4) && iy>=seed_map->sizeY()-buf_size) value = 1; // west |
584 | seed_map->valueAtIndex(ix,iy)=value; |
||
491 | werner | 585 | } |
586 | } |
||
472 | werner | 587 | } else { |
588 | qDebug() << "external seed input: Error: invalid buffer size???"; |
||
589 | } |
||
590 | } |
||
375 | werner | 591 | } |
592 | |||
1176 | werner | 593 | static int _debug_ldd=0; |
375 | werner | 594 | void SeedDispersal::execute() |
595 | { |
||
991 | werner | 596 | #ifdef ILAND_GUI |
597 | int year = GlobalSettings::instance()->currentYear(); |
||
598 | QString path; |
||
472 | werner | 599 | if (mDumpSeedMaps) { |
1102 | werner | 600 | path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 601 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_before_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
618 | werner | 602 | qDebug() << "saved seed map image to" << path; |
992 | werner | 603 | } |
989 | werner | 604 | #else |
992 | werner | 605 | if (mDumpSeedMaps) |
989 | werner | 606 | qDebug() << "saving of seedmaps only supported in the iLand GUI."; |
607 | #endif |
||
1180 | werner | 608 | if (mProbMode) { |
992 | werner | 609 | |
1180 | werner | 610 | DebugTimer t("seed dispersal", true); |
1106 | werner | 611 | |
1180 | werner | 612 | // (1) detect edges |
613 | if (edgeDetection()) { |
||
614 | |||
1106 | werner | 615 | #ifdef ILAND_GUI |
1180 | werner | 616 | if (mDumpSeedMaps) { |
617 | gridToImage(seedMap(), true, -1., 1.).save(QString("%1/seed_edge_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
618 | } |
||
1106 | werner | 619 | #endif |
620 | |||
1180 | werner | 621 | // (2) distribute seed probabilites from edges |
622 | distribute(); |
||
623 | } |
||
1167 | werner | 624 | |
1180 | werner | 625 | // special case serotiny |
626 | if (mHasPendingSerotiny) { |
||
627 | qDebug() << "calculating extra seed rain (serotiny)...."; |
||
1168 | werner | 628 | #ifdef ILAND_GUI |
1180 | werner | 629 | if (mDumpSeedMaps) { |
630 | gridToImage(mSeedMapSerotiny, true, 0., 1.).save(QString("%1/seed_serotiny_before_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
631 | } |
||
1168 | werner | 632 | #endif |
1180 | werner | 633 | if (edgeDetection(&mSeedMapSerotiny)) |
634 | distribute(&mSeedMapSerotiny); |
||
635 | // copy back data |
||
636 | float *sero=mSeedMapSerotiny.begin(); |
||
637 | for (float* p=mSeedMap.begin();p!=mSeedMap.end();++p, ++sero) |
||
638 | *p = std::max(*p, *sero); |
||
1167 | werner | 639 | |
1180 | werner | 640 | float total = mSeedMapSerotiny.sum(); |
1168 | werner | 641 | #ifdef ILAND_GUI |
1180 | werner | 642 | if (mDumpSeedMaps) { |
643 | gridToImage(mSeedMapSerotiny, true, 0., 1.).save(QString("%1/seed_serotiny_after_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
644 | } |
||
645 | #endif |
||
646 | mSeedMapSerotiny.initialize(0.f); // clear |
||
647 | mHasPendingSerotiny = false; |
||
648 | qDebug() << "serotiny event: extra seed input" << total << "(total sum of seed probability over all pixels of the serotiny seed map) of species" << mSpecies->name(); |
||
1168 | werner | 649 | } |
1180 | werner | 650 | |
651 | } else { |
||
652 | // distribute actual values |
||
653 | DebugTimer t("seed dispersal", true); |
||
654 | // fill seed map from source map |
||
655 | distributeSeeds(); |
||
656 | |||
619 | werner | 657 | } |
989 | werner | 658 | #ifdef ILAND_GUI |
1102 | werner | 659 | if (mDumpSeedMaps) { |
1180 | werner | 660 | //qDebug() << "finished seed dispersal for species. time: " << mSpecies->id() << t.elapsed(); |
472 | werner | 661 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_after_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
1102 | werner | 662 | } |
1064 | werner | 663 | |
664 | if (!mDumpNextYearFileName.isEmpty()) { |
||
665 | Helper::saveToTextFile(GlobalSettings::instance()->path(mDumpNextYearFileName), gridToESRIRaster(seedMap())); |
||
666 | qDebug() << "saved seed map for " << species()->id() << "to" << GlobalSettings::instance()->path(mDumpNextYearFileName); |
||
667 | mDumpNextYearFileName = QString(); |
||
668 | } |
||
1176 | werner | 669 | qDebug() << "LDD-count:" << _debug_ldd; |
1064 | werner | 670 | |
989 | werner | 671 | #endif |
375 | werner | 672 | } |
673 | |||
373 | werner | 674 | /** scans the seed image and detects "edges". |
675 | edges are then subsequently marked (set to -1). This is pass 1 of the seed distribution process. |
||
676 | */ |
||
1167 | werner | 677 | bool SeedDispersal::edgeDetection(Grid<float> *seed_map) |
373 | werner | 678 | { |
679 | float *p_above, *p, *p_below; |
||
1167 | werner | 680 | Grid<float> &seedmap = seed_map ? *seed_map : mSeedMap; // switch to extra seed map if provided |
681 | int dy = seedmap.sizeY(); |
||
682 | int dx = seedmap.sizeX(); |
||
373 | werner | 683 | int x,y; |
619 | werner | 684 | bool found = false; |
1106 | werner | 685 | |
686 | // fill mini-gaps |
||
1168 | werner | 687 | int n_gaps_filled=0; |
373 | werner | 688 | for (y=1;y<dy-1;++y){ |
1167 | werner | 689 | p = seedmap.ptr(1,y); |
373 | werner | 690 | p_above = p - dx; // one line above |
691 | p_below = p + dx; // one line below |
||
692 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
1106 | werner | 693 | if (*p < 0.999f) { |
694 | |||
695 | if ((*(p_above-1)==1.f) + (*p_above==1.f) + (*(p_above+1)==1.f) + |
||
696 | (*(p-1)==1.f) + (*(p+1)==1.f) + |
||
1168 | werner | 697 | (*(p_below-1)==1.f) + (*p_below==1.f) + (*(p_below+1)==1.f) > 3) { |
1106 | werner | 698 | *p=0.999f; // if more than 3 neighbors are active pixels, the value is high |
1168 | werner | 699 | ++n_gaps_filled; |
700 | } |
||
1106 | werner | 701 | } |
702 | |||
703 | } |
||
704 | } |
||
705 | |||
706 | |||
707 | // now detect the edges |
||
708 | int n_edges=0 ; |
||
709 | for (y=1;y<dy-1;++y){ |
||
1168 | werner | 710 | p = seedmap.ptr(1,y); |
1106 | werner | 711 | p_above = p - dx; // one line above |
712 | p_below = p + dx; // one line below |
||
713 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
714 | if (*p == 1.f) { |
||
619 | werner | 715 | found = true; |
1106 | werner | 716 | if ( (*(p_above-1)<0.999f && *(p_above-1)>=0.f) |
717 | || (*p_above<0.999f && *p_above>=0.f) |
||
718 | || (*(p_above+1)<0.999f && *(p_above+1)>=0.f) |
||
719 | || (*(p-1)<0.999f && *(p-1)>=0.f) |
||
720 | || (*(p+1)<0.999f && (*p+1)>=0.f) |
||
721 | || (*(p_below-1)<0.999f && *(p_below-1)>=0.f) |
||
722 | || (*p_below<0.999f && *p_below>=0.f) |
||
723 | || (*(p_below+1)<0.999f && *(p_below+1)>=0.f ) ) { |
||
724 | *p=-1.f; // if any surrounding pixel is >=0 & <0.999: -> mark as edge |
||
725 | ++n_edges; |
||
726 | } |
||
373 | werner | 727 | } |
728 | |||
729 | } |
||
730 | } |
||
1168 | werner | 731 | if (mDumpSeedMaps) |
732 | qDebug() << "species:" << mSpecies->id() << "# of gaps filled: " << n_gaps_filled << "# of edge-pixels:" << n_edges; |
||
619 | werner | 733 | return found; |
373 | werner | 734 | } |
735 | |||
736 | /** do the seed probability distribution. |
||
737 | This is phase 2. Apply the seed kernel for each "edge" point identified in phase 1. |
||
738 | */ |
||
1167 | werner | 739 | void SeedDispersal::distribute(Grid<float> *seed_map) |
373 | werner | 740 | { |
741 | int x,y; |
||
1167 | werner | 742 | Grid<float> &seedmap = seed_map ? *seed_map : mSeedMap; // switch to extra seed map if provided |
743 | float *end = seedmap.end(); |
||
744 | float *p = seedmap.begin(); |
||
415 | werner | 745 | // choose the kernel depending whether there is a seed year for the current species or not |
1168 | werner | 746 | Grid<float> *kernel = species()->isSeedYear()? &mKernelSeedYear : &mKernelNonSeedYear; |
1167 | werner | 747 | // extra case: serotiny |
748 | if (seed_map) |
||
1168 | werner | 749 | kernel = &mKernelSerotiny; |
1167 | werner | 750 | |
1168 | werner | 751 | int offset = kernel->sizeX() / 2; // offset is the index of the center pixel |
373 | werner | 752 | for(;p!=end;++p) { |
375 | werner | 753 | if (*p==-1.f) { |
373 | werner | 754 | // edge pixel found. Now apply the kernel.... |
1167 | werner | 755 | QPoint pt=seedmap.indexOf(p); |
1168 | werner | 756 | for (y=-offset;y<=offset;++y) { |
757 | for (x=-offset;x<=offset;++x) { |
||
758 | float &kernel_value = kernel->valueAtIndex(x+offset, y+offset); |
||
759 | if (kernel_value>0.f && seedmap.isIndexValid(pt.x()+x, pt.y()+y)) { |
||
1167 | werner | 760 | float &val = seedmap.valueAtIndex(pt.x()+x, pt.y()+y); |
1106 | werner | 761 | if (val!=-1.f) |
1168 | werner | 762 | val = qMin(1.f - (1.f - val)*(1.f-kernel_value),1.f ); |
373 | werner | 763 | } |
1168 | werner | 764 | } |
765 | } |
||
1176 | werner | 766 | // long distance dispersal |
767 | if (!mLDDDensity.isEmpty()) { |
||
768 | double m = species()->isSeedYear() ? 1. : mNonSeedYearFraction; |
||
769 | for (int r=0;r<mLDDDensity.size(); ++r) { |
||
1180 | werner | 770 | float ldd_val = mLDDSeedlings; // pixels will have this probability |
1176 | werner | 771 | int n = round( mLDDDensity[r]*m ); // number of pixels to activate |
772 | for (int i=0;i<n;++i) { |
||
773 | // distance and direction: |
||
774 | double radius = nrandom(mLDDDistance[r], mLDDDistance[r+1]) / seedmap.cellsize(); // choose a random distance (in pixels) |
||
775 | double phi = drandom()*2.*M_PI; // choose a random direction |
||
776 | QPoint ldd(pt.x() + radius*cos(phi), pt.y() + radius*sin(phi)); |
||
777 | if (seedmap.isIndexValid(ldd)) { |
||
778 | float &val = seedmap.valueAtIndex(ldd); |
||
779 | _debug_ldd++; |
||
780 | // use the same adding of probabilities |
||
781 | if (val!=-1.f) |
||
782 | val = qMin(1.f - (1.f - val)*(1.f-ldd_val), 1.f); |
||
783 | } |
||
784 | } |
||
785 | } |
||
786 | } |
||
375 | werner | 787 | *p=1.f; // mark as processed |
373 | werner | 788 | } // *p==1 |
789 | } // for() |
||
790 | } |
||
1180 | werner | 791 | |
1182 | werner | 792 | // because C modulo operation gives negative numbers for negative values, here a fix |
793 | // that always returns positive numbers: http://www.lemoda.net/c/modulo-operator/ |
||
794 | #define MOD(a,b) ((((a)%(b))+(b))%(b)) |
||
795 | |||
1180 | werner | 796 | void SeedDispersal::distributeSeeds(Grid<float> *seed_map) |
797 | { |
||
798 | Grid<float> &sourcemap = seed_map ? *seed_map : mSourceMap; // switch to extra seed map if provided |
||
799 | Grid<float> &kernel = mKernelSeedYear; |
||
1182 | werner | 800 | |
801 | // *** estimate seed production (based on leaf area) *** |
||
1180 | werner | 802 | // calculate number of seeds; the source map holds now m2 leaf area on 20x20m pixels |
803 | // after this step, each source cell has a value between 0 (no source) and 1 (fully covered cell) |
||
804 | float fec = species()->fecundity_m2(); |
||
805 | if (!species()->isSeedYear()) |
||
806 | fec *= mNonSeedYearFraction; |
||
807 | for (float *p=sourcemap.begin(); p!=sourcemap.end(); ++p){ |
||
808 | if (*p) { |
||
809 | // if LAI >3, then full potential is assumed, below LAI=3 a linear ramp is used |
||
810 | *p = std::min(*p / (sourcemap.cellsize()*sourcemap.cellsize()) /3.f, 3.f); |
||
811 | } |
||
812 | } |
||
813 | |||
814 | // sink mode |
||
815 | |||
1182 | werner | 816 | // // now look for each pixel in the targetmap and sum up seeds*kernel |
817 | // int idx=0; |
||
818 | // int offset = kernel.sizeX() / 2; // offset is the index of the center pixel |
||
819 | // //const Grid<ResourceUnit*> &ru_map = GlobalSettings::instance()->model()->RUgrid(); |
||
820 | // DebugTimer tsink("seed_sink"); { |
||
821 | // for (float *t=mSeedMap.begin(); t!=mSeedMap.end(); ++t, ++idx) { |
||
822 | // // skip out-of-project areas |
||
823 | // //if (!ru_map.constValueAtIndex(mSeedMap.index5(idx))) |
||
824 | // // continue; |
||
825 | // // apply the kernel |
||
826 | // QPoint sm=mSeedMap.indexOf(t)-QPoint(offset, offset); |
||
827 | // for (int iy=0;iy<kernel.sizeY();++iy) { |
||
828 | // for (int ix=0;ix<kernel.sizeX();++ix) { |
||
829 | // if (sourcemap.isIndexValid(sm.x()+ix, sm.y()+iy)) |
||
830 | // *t+=sourcemap(sm.x()+ix, sm.y()+iy) * kernel(ix, iy); |
||
831 | // } |
||
832 | // } |
||
833 | // } |
||
834 | // } // debugtimer |
||
835 | // mSeedMap.initialize(0.f); // just for debugging... |
||
1180 | werner | 836 | |
837 | int offset = kernel.sizeX() / 2; // offset is the index of the center pixel |
||
838 | // source mode |
||
839 | |||
1182 | werner | 840 | // *** seed distribution (Kernel + long distance dispersal) *** |
841 | if (GlobalSettings::instance()->model()->settings().torusMode==false) { |
||
842 | // ** standard case (no torus) ** |
||
843 | for (float *src=sourcemap.begin(); src!=sourcemap.end(); ++src) { |
||
844 | if (*src>0.f) { |
||
845 | QPoint sm=sourcemap.indexOf(src)-QPoint(offset, offset); |
||
846 | int sx = sm.x(), sy=sm.y(); |
||
847 | for (int iy=0;iy<kernel.sizeY();++iy) { |
||
848 | for (int ix=0;ix<kernel.sizeX();++ix) { |
||
849 | if (mSeedMap.isIndexValid(sx+ix, sy+iy)) |
||
850 | mSeedMap.valueAtIndex(sx+ix, sy+iy)+= *src * kernel(ix, iy); |
||
851 | } |
||
852 | } |
||
853 | // long distance dispersal |
||
854 | if (!mLDDDensity.isEmpty()) { |
||
855 | QPoint pt=sourcemap.indexOf(src); |
||
1180 | werner | 856 | |
1182 | werner | 857 | for (int r=0;r<mLDDDensity.size(); ++r) { |
858 | float ldd_val = mLDDSeedlings / fec; // pixels will have this probability [note: fecundity will be multiplied below] |
||
859 | int n; |
||
860 | if (mLDDDensity[r]<1) |
||
861 | n = drandom()<mLDDDensity[r] ? 1 : 0; |
||
862 | else |
||
863 | n = round( mLDDDensity[r] ); // number of pixels to activate |
||
864 | for (int i=0;i<n;++i) { |
||
865 | // distance and direction: |
||
866 | double radius = nrandom(mLDDDistance[r], mLDDDistance[r+1]) / mSeedMap.cellsize(); // choose a random distance (in pixels) |
||
867 | double phi = drandom()*2.*M_PI; // choose a random direction |
||
868 | QPoint ldd(pt.x() + radius*cos(phi), pt.y() + radius*sin(phi)); |
||
869 | if (mSeedMap.isIndexValid(ldd)) { |
||
870 | float &val = mSeedMap.valueAtIndex(ldd); |
||
871 | _debug_ldd++; |
||
872 | val += ldd_val; |
||
873 | } |
||
874 | } |
||
875 | } |
||
1180 | werner | 876 | } |
1182 | werner | 877 | |
1180 | werner | 878 | } |
1182 | werner | 879 | } |
880 | } else { |
||
881 | // **** seed distribution in torus mode *** |
||
882 | int seedmap_offset = sourcemap.indexAt(QPointF(0., 0.)).x(); // the seed maps have x extra rows/columns |
||
883 | QPoint torus_pos; |
||
884 | int seedpx_per_ru = static_cast<int>((cRUSize/sourcemap.cellsize())); |
||
885 | for (float *src=sourcemap.begin(); src!=sourcemap.end(); ++src) { |
||
886 | if (*src>0.f) { |
||
887 | QPoint sm=sourcemap.indexOf(src); |
||
888 | // get the origin of the resource unit *on* the seedmap in *seedmap-coords*: |
||
889 | QPoint offset_ru( ((sm.x()-seedmap_offset) / seedpx_per_ru) * seedpx_per_ru + seedmap_offset, |
||
890 | ((sm.y()-seedmap_offset) / seedpx_per_ru) * seedpx_per_ru + seedmap_offset); // coords RU origin |
||
1180 | werner | 891 | |
1182 | werner | 892 | QPoint offset_in_ru((sm.x()-seedmap_offset) % seedpx_per_ru, (sm.y()-seedmap_offset) % seedpx_per_ru ); // offset of current point within the RU |
893 | |||
894 | //QPoint sm=sourcemap.indexOf(src)-QPoint(offset, offset); |
||
895 | for (int iy=0;iy<kernel.sizeY();++iy) { |
||
896 | for (int ix=0;ix<kernel.sizeX();++ix) { |
||
897 | torus_pos = offset_ru + QPoint(MOD((offset_in_ru.x() - offset + ix), seedpx_per_ru), MOD((offset_in_ru.y() - offset + iy), seedpx_per_ru)); |
||
898 | |||
899 | if (mSeedMap.isIndexValid(torus_pos)) |
||
900 | mSeedMap.valueAtIndex(torus_pos)+= *src * kernel(ix, iy); |
||
901 | } |
||
902 | } |
||
903 | // long distance dispersal |
||
904 | if (!mLDDDensity.isEmpty()) { |
||
905 | |||
906 | for (int r=0;r<mLDDDensity.size(); ++r) { |
||
907 | float ldd_val = mLDDSeedlings / fec; // pixels will have this probability [note: fecundity will be multiplied below] |
||
908 | int n; |
||
909 | if (mLDDDensity[r]<1) |
||
910 | n = drandom()<mLDDDensity[r] ? 1 : 0; |
||
911 | else |
||
912 | n = round( mLDDDensity[r] ); // number of pixels to activate |
||
913 | for (int i=0;i<n;++i) { |
||
914 | // distance and direction: |
||
915 | double radius = nrandom(mLDDDistance[r], mLDDDistance[r+1]) / mSeedMap.cellsize(); // choose a random distance (in pixels) |
||
916 | double phi = drandom()*2.*M_PI; // choose a random direction |
||
917 | QPoint ldd( radius*cos(phi), + radius*sin(phi)); // destination (offset) |
||
918 | torus_pos = offset_ru + QPoint(MOD((offset_in_ru.x()+ldd.x()),seedpx_per_ru), MOD((offset_in_ru.y()+ldd.y()),seedpx_per_ru) ); |
||
919 | |||
920 | if (mSeedMap.isIndexValid(torus_pos)) { |
||
921 | float &val = mSeedMap.valueAtIndex(torus_pos); |
||
922 | _debug_ldd++; |
||
923 | val += ldd_val; |
||
924 | } |
||
1180 | werner | 925 | } |
926 | } |
||
927 | } |
||
1182 | werner | 928 | |
1180 | werner | 929 | } |
930 | } |
||
1182 | werner | 931 | } // torus |
1180 | werner | 932 | |
933 | |||
934 | |||
935 | // now the seed sources (0..1) are spatially distributed by the kernel (and LDD) without altering the magnitude; |
||
936 | // now we include the fecundity (=seedling potential per m2 crown area), and convert to the establishment probability p_seed. |
||
937 | // The number of (potential) seedlings per m2 on each cell is: cell * fecundity[m2] |
||
938 | // We assume that the availability of 10 potential seedlings/m2 is enough for unconstrained establishment; |
||
1181 | werner | 939 | const float n_unlimited = 100.f; |
1180 | werner | 940 | for (float *p=mSeedMap.begin(); p!=mSeedMap.end(); ++p){ |
941 | if (*p>0.f) { |
||
942 | *p = std::min(*p*fec / n_unlimited, 1.f); |
||
943 | } |
||
944 | } |
||
945 | } |