Rev 1178 | Rev 1180 | 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 |
1176 | werner | 82 | // copy values for the species parameters: |
83 | mSpecies->treeMigKernel(mTM_as1, mTM_as2, mTM_ks); |
||
445 | werner | 84 | mTM_fecundity_cell = mSpecies->fecundity_m2() * seedmap_size*seedmap_size * mTM_occupancy; // scale to production for the whole cell |
85 | mNonSeedYearFraction = mSpecies->nonSeedYearFraction(); |
||
1176 | werner | 86 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.settings.seedDispersal")); |
87 | mKernelThresholdArea = xml.valueDouble(".longDistanceDispersal.thresholdArea", 0.0001); |
||
88 | mKernelThresholdLDD = xml.valueDouble(".longDistanceDispersal.thresholdLDD", 0.0001); |
||
89 | mLDDProbability = xml.valueDouble(".longDistanceDispersal.minProbability", 0.0001); |
||
90 | mLDDRings = xml.valueInt(".longDistanceDispersal.rings", 4); |
||
415 | werner | 91 | |
1176 | werner | 92 | mLDDProbability = qMax(mLDDProbability, static_cast<float>(mKernelThresholdArea)); |
415 | werner | 93 | |
1176 | werner | 94 | |
95 | mFecundityFactor = createKernel(mKernelSeedYear, mTM_fecundity_cell); |
||
96 | |||
415 | werner | 97 | // the kernel for non seed years looks similar, but is simply linearly scaled down |
98 | // using the species parameter NonSeedYearFraction. |
||
99 | // the central pixel still gets the value of 1 (i.e. 100% probability) |
||
445 | werner | 100 | createKernel(mKernelNonSeedYear, mTM_fecundity_cell*mNonSeedYearFraction); |
415 | werner | 101 | |
1167 | werner | 102 | if (mSpecies->fecunditySerotiny()>0.) { |
103 | // an extra seed map is used for storing information related to post-fire seed rain |
||
1168 | werner | 104 | mSeedMapSerotiny.clear(); |
105 | mSeedMapSerotiny.setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
||
106 | mSeedMapSerotiny.initialize(0.); |
||
1167 | werner | 107 | |
108 | // set up the special seed kernel for post fire seed rain |
||
109 | createKernel(mKernelSerotiny, mTM_fecundity_cell * mSpecies->fecunditySerotiny()); |
||
110 | qDebug() << "created extra seed map and serotiny seed kernel for species" << mSpecies->name() << "with fecundity factor" << mSpecies->fecunditySerotiny(); |
||
111 | } |
||
112 | mHasPendingSerotiny = false; |
||
113 | |||
472 | werner | 114 | // debug info |
115 | mDumpSeedMaps = GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false); |
||
116 | if (mDumpSeedMaps) { |
||
1102 | werner | 117 | QString path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 118 | Helper::saveToTextFile(QString("%1/seedkernelYes_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelSeedYear)); |
119 | Helper::saveToTextFile(QString("%1/seedkernelNo_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelNonSeedYear)); |
||
1168 | werner | 120 | if (!mKernelSerotiny.isEmpty()) |
121 | Helper::saveToTextFile(QString("%1/seedkernelSerotiny_%2.csv").arg(path).arg(mSpecies->id()),gridToString(mKernelSerotiny)); |
||
417 | werner | 122 | } |
1176 | werner | 123 | // long distance dispersal |
124 | setupLDD(); |
||
125 | |||
472 | werner | 126 | // external seeds |
481 | werner | 127 | mHasExternalSeedInput = false; |
491 | werner | 128 | mExternalSeedBuffer = 0; |
129 | mExternalSeedDirection = 0; |
||
836 | werner | 130 | mExternalSeedBackgroundInput = 0.; |
472 | werner | 131 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.externalSeedEnabled",false)) { |
764 | werner | 132 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) { |
133 | // external seed input specified by sectors and around the project area (seedbelt) |
||
134 | setupExternalSeedsForSpecies(mSpecies); |
||
135 | } else { |
||
136 | // external seeds specified fixedly per cardinal direction |
||
137 | // current species in list?? |
||
138 | mHasExternalSeedInput = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSpecies").contains(mSpecies->id()); |
||
139 | QString dir = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedSource").toLower(); |
||
140 | // encode cardinal positions as bits: e.g: "e,w" -> 6 |
||
141 | mExternalSeedDirection += dir.contains("n")?1:0; |
||
142 | mExternalSeedDirection += dir.contains("e")?2:0; |
||
143 | mExternalSeedDirection += dir.contains("s")?4:0; |
||
144 | mExternalSeedDirection += dir.contains("w")?8:0; |
||
837 | werner | 145 | QStringList buffer_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBuffer").split(QRegExp("([^\\.\\w]+)")); |
764 | werner | 146 | int index = buffer_list.indexOf(mSpecies->id()); |
147 | if (index>=0) { |
||
148 | mExternalSeedBuffer = buffer_list[index+1].toInt(); |
||
149 | qDebug() << "enabled special buffer for species" <<mSpecies->id() << ": distance of" << mExternalSeedBuffer << "pixels = " << mExternalSeedBuffer*20. << "m"; |
||
150 | } |
||
836 | werner | 151 | |
152 | // background seed rain (i.e. for the full landscape), use regexp |
||
837 | werner | 153 | QStringList background_input_list = GlobalSettings::instance()->settings().value("model.settings.seedDispersal.externalSeedBackgroundInput").split(QRegExp("([^\\.\\w]+)")); |
836 | werner | 154 | index = background_input_list.indexOf(mSpecies->id()); |
155 | if (index>=0) { |
||
156 | mExternalSeedBackgroundInput = background_input_list[index+1].toDouble(); |
||
157 | qDebug() << "enabled background seed input (for full area) for species" <<mSpecies->id() << ": p=" << mExternalSeedBackgroundInput; |
||
158 | } |
||
159 | |||
764 | werner | 160 | if (mHasExternalSeedInput) |
161 | qDebug() << "External seed input enabled for" << mSpecies->id(); |
||
491 | werner | 162 | } |
472 | werner | 163 | } |
415 | werner | 164 | |
373 | werner | 165 | // setup of seed kernel |
391 | werner | 166 | // const int max_radius = 15; // pixels |
167 | // |
||
168 | // mSeedKernel.clear(); |
||
169 | // mSeedKernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
||
170 | // mKernelOffset = max_radius; |
||
171 | // // filling of the kernel.... for simplicity: a linear kernel |
||
172 | // QPoint center = QPoint(mKernelOffset, mKernelOffset); |
||
173 | // const double max_dist = max_radius * seedmap_size; |
||
174 | // for (float *p=mSeedKernel.begin(); p!=mSeedKernel.end();++p) { |
||
175 | // double d = mSeedKernel.distance(center, mSeedKernel.indexOf(p)); |
||
176 | // *p = qMax( 1. - d / max_dist, 0.); |
||
177 | // } |
||
373 | werner | 178 | |
179 | |||
180 | // randomize seed map.... set 1/3 to "filled" |
||
375 | werner | 181 | //for (int i=0;i<mSeedMap.count(); i++) |
182 | // mSeedMap.valueAtIndex(mSeedMap.randomPosition()) = 1.; |
||
373 | werner | 183 | |
184 | |||
375 | werner | 185 | // QImage img = gridToImage(mSeedMap, true, -1., 1.); |
186 | // img.save("seedmap.png"); |
||
187 | |||
188 | // img = gridToImage(mSeedMap, true, -1., 1.); |
||
764 | werner | 189 | // img.save("seedmap_e.png"); |
373 | werner | 190 | } |
191 | |||
764 | werner | 192 | void SeedDispersal::setupExternalSeeds() |
193 | { |
||
194 | mExternalSeedBaseMap = 0; |
||
195 | if (!GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.seedBelt.enabled",false)) |
||
196 | return; |
||
197 | |||
978 | werner | 198 | DebugTimer t("setup of external seed maps."); |
764 | werner | 199 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.settings.seedDispersal.seedBelt")); |
1102 | werner | 200 | int seedbelt_width =xml.valueInt(".width",10); |
764 | werner | 201 | // setup of sectors |
202 | // setup of base map |
||
203 | const float seedmap_size = 20.f; |
||
204 | mExternalSeedBaseMap = new Grid<float>; |
||
205 | mExternalSeedBaseMap->setup(GlobalSettings::instance()->model()->heightGrid()->metricRect(), seedmap_size ); |
||
206 | mExternalSeedBaseMap->initialize(0.); |
||
207 | if (mExternalSeedBaseMap->count()*4 != GlobalSettings::instance()->model()->heightGrid()->count()) |
||
947 | werner | 208 | 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 | 209 | // make a copy of the 10m height grid in lower resolution and mark pixels that are forested and outside of |
210 | // the project area. |
||
211 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) |
||
212 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
765 | werner | 213 | bool val = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isForestOutside(); |
764 | werner | 214 | mExternalSeedBaseMap->valueAtIndex(x,y) = val?1.f:0.f; |
215 | if(GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(x*2,y*2).isValid()) |
||
216 | mExternalSeedBaseMap->valueAtIndex(x,y) = -1.f; |
||
217 | } |
||
836 | werner | 218 | QString path = GlobalSettings::instance()->path(GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath")); |
765 | werner | 219 | |
220 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
||
989 | werner | 221 | #ifdef ILAND_GUI |
765 | werner | 222 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
223 | img.save(path + "/seedbeltmap_before.png"); |
||
989 | werner | 224 | #endif |
765 | werner | 225 | } |
764 | werner | 226 | // img.save("seedmap.png"); |
227 | // now scan the pixels of the belt: paint all pixels that are close to the project area |
||
228 | // we do this 4 times (for all cardinal direcitons) |
||
229 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
230 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
231 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
232 | continue; |
||
765 | werner | 233 | int look_forward = std::min(x + seedbelt_width, mExternalSeedBaseMap->sizeX()-1); |
1106 | werner | 234 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 235 | // fill pixels |
236 | for(; x<look_forward;++x) { |
||
237 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
238 | if (v==1.f) v=2.f; |
||
239 | } |
||
240 | } |
||
241 | } |
||
242 | } |
||
243 | // right to left |
||
244 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
245 | for (int x=mExternalSeedBaseMap->sizeX();x>=0;--x) { |
||
246 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
247 | continue; |
||
248 | int look_forward = std::max(x - seedbelt_width, 0); |
||
1106 | werner | 249 | if (mExternalSeedBaseMap->valueAtIndex(look_forward, y)==-1.f) { |
764 | werner | 250 | // fill pixels |
251 | for(; x>look_forward;--x) { |
||
252 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
253 | if (v==1.f) v=2.f; |
||
254 | } |
||
255 | } |
||
256 | } |
||
257 | } |
||
258 | // up and down *** |
||
259 | // from top to bottom |
||
260 | for (int x=0;x<mExternalSeedBaseMap->sizeX();x++) { |
||
261 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
262 | |||
263 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
264 | continue; |
||
765 | werner | 265 | int look_forward = std::min(y + seedbelt_width, mExternalSeedBaseMap->sizeY()-1); |
764 | werner | 266 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
267 | // fill pixels |
||
268 | for(; y<look_forward;++y) { |
||
269 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
270 | if (v==1.f) v=2.f; |
||
271 | } |
||
272 | } |
||
273 | } |
||
274 | } |
||
275 | // bottom to top *** |
||
276 | for (int y=0;y<mExternalSeedBaseMap->sizeY();y++) { |
||
277 | for (int x=mExternalSeedBaseMap->sizeX();x>=0;--x) { |
||
278 | if (mExternalSeedBaseMap->valueAtIndex(x, y)!=1.) |
||
279 | continue; |
||
280 | int look_forward = std::max(y - seedbelt_width, 0); |
||
281 | if (mExternalSeedBaseMap->valueAtIndex(x, look_forward)==-1.) { |
||
282 | // fill pixels |
||
283 | for(; y>look_forward;--y) { |
||
284 | float &v = mExternalSeedBaseMap->valueAtIndex(x, y); |
||
285 | if (v==1.f) v=2.f; |
||
286 | } |
||
287 | } |
||
288 | } |
||
289 | } |
||
765 | werner | 290 | if (GlobalSettings::instance()->settings().valueBool("model.settings.seedDispersal.dumpSeedMapsEnabled",false)) { |
989 | werner | 291 | #ifdef ILAND_GUI |
765 | werner | 292 | QImage img = gridToImage(*mExternalSeedBaseMap, true, -1., 2.); |
972 | werner | 293 | img.save(path + "/seedbeltmap_after.png"); |
989 | werner | 294 | #endif |
765 | werner | 295 | } |
764 | werner | 296 | mExtSeedData.clear(); |
1102 | werner | 297 | int sectors_x = xml.valueInt("sizeX",0); |
298 | int sectors_y = xml.valueInt("sizeY",0); |
||
764 | werner | 299 | if(sectors_x<1 || sectors_y<1) |
300 | throw IException(QString("setup of external seed dispersal: invalid number of sectors x=%1 y=%3").arg(sectors_x).arg(sectors_y)); |
||
301 | QDomElement elem = xml.node("."); |
||
302 | for(QDomNode n = elem.firstChild(); !n.isNull(); n = n.nextSibling()) { |
||
303 | if (n.nodeName().startsWith("species")) { |
||
304 | QStringList coords = n.nodeName().split("_"); |
||
305 | if (coords.count()!=3) |
||
306 | throw IException("external seed species definition is not valid: " + n.nodeName()); |
||
307 | int x = coords[1].toInt(); |
||
308 | int y = coords[2].toInt(); |
||
309 | if (x<0 || x>=sectors_x || y<0 || y>=sectors_y) |
||
310 | throw IException(QString("invalid sector for specifiing external seed input (x y): %1 %2 ").arg(x).arg(y) ); |
||
311 | int index = y*sectors_x + x; |
||
312 | |||
313 | QString text = xml.value("." + n.nodeName()); |
||
314 | qDebug() << "processing element " << n.nodeName() << "x,y:" << x << y << text; |
||
315 | // we assume pairs of name and fraction |
||
316 | QStringList species_list = text.split(" "); |
||
317 | for (int i=0;i<species_list.count();++i) { |
||
318 | QVector<double> &space = mExtSeedData[species_list[i]]; |
||
319 | if (space.isEmpty()) |
||
320 | space.resize(sectors_x*sectors_y); // are initialized to 0s |
||
321 | double fraction = species_list[++i].toDouble(); |
||
322 | space[index] = fraction; |
||
323 | } |
||
324 | } |
||
325 | } |
||
326 | mExtSeedSizeX = sectors_x; |
||
327 | mExtSeedSizeY = sectors_y; |
||
328 | qDebug() << "setting up of external seed maps finished"; |
||
329 | } |
||
330 | |||
331 | void SeedDispersal::finalizeExternalSeeds() |
||
332 | { |
||
333 | if (mExternalSeedBaseMap) |
||
334 | delete mExternalSeedBaseMap; |
||
335 | mExternalSeedBaseMap = 0; |
||
336 | } |
||
337 | |||
1167 | werner | 338 | void SeedDispersal::seedProductionSerotiny(const QPoint &position_index) |
339 | { |
||
1168 | werner | 340 | if (mSeedMapSerotiny.isEmpty()) |
1167 | werner | 341 | throw IException("Invalid use seedProductionSerotiny(): tried to set a seed source for a non-serotinous species!"); |
1168 | werner | 342 | mSeedMapSerotiny.valueAtIndex(position_index.x()/mIndexFactor, position_index.y()/mIndexFactor)=1.f; |
1167 | werner | 343 | mHasPendingSerotiny = true; |
344 | } |
||
345 | |||
391 | werner | 346 | // ************ Kernel ************** |
1176 | werner | 347 | float SeedDispersal::createKernel(Grid<float> &kernel, const double max_seed) |
391 | werner | 348 | { |
415 | werner | 349 | |
1176 | werner | 350 | double max_dist = treemig_distanceTo(mKernelThresholdArea); |
391 | werner | 351 | double cell_size = mSeedMap.cellsize(); |
415 | werner | 352 | int max_radius = int(max_dist / cell_size); |
391 | werner | 353 | // e.g.: cell_size: regeneration grid (e.g. 400qm), px-size: light-grid (4qm) |
445 | werner | 354 | double occupation = cell_size*cell_size / (cPxSize*cPxSize * mTM_occupancy); |
391 | werner | 355 | |
415 | werner | 356 | kernel.clear(); |
391 | werner | 357 | |
415 | werner | 358 | kernel.setup(mSeedMap.cellsize(), 2*max_radius + 1 , 2*max_radius + 1); |
359 | int kernel_offset = max_radius; |
||
360 | |||
391 | werner | 361 | // filling of the kernel.... use the treemig |
1179 | werner | 362 | double dist_center_cell = sqrt(cell_size*cell_size/M_PI); |
415 | werner | 363 | QPoint center = QPoint(kernel_offset, kernel_offset); |
364 | const float *sk_end = kernel.end(); |
||
365 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
||
366 | double d = kernel.distance(center, kernel.indexOf(p)); |
||
1178 | werner | 367 | if (d==0.) |
1179 | werner | 368 | *p = treemig_centercell(dist_center_cell); // r is the radius of a circle with the same area as a cell |
1178 | werner | 369 | else |
1179 | werner | 370 | *p = d<=max_dist?static_cast<float>(( treemig(d+dist_center_cell) + treemig(d+dist_center_cell))/2. ):0.f; |
391 | werner | 371 | } |
372 | |||
373 | // normalize |
||
1102 | werner | 374 | float sum = kernel.sum(); |
391 | werner | 375 | if (sum==0. || occupation==0.) |
376 | throw IException("create seed kernel: sum of probabilities = 0!"); |
||
377 | |||
378 | // the sum of all kernel cells has to equal 1 |
||
1178 | werner | 379 | // kernel.multiply(1.f/sum); |
380 | |||
1176 | werner | 381 | float fecundity_factor = static_cast<float>( max_seed / occupation); |
391 | werner | 382 | // probabilities are derived in multiplying by seed number, and dividing by occupancy criterion |
1176 | werner | 383 | kernel.multiply( fecundity_factor ); |
391 | werner | 384 | // all cells that get more seeds than the occupancy criterion are considered to have no seed limitation for regeneration |
415 | werner | 385 | for (float *p=kernel.begin(); p!=sk_end;++p) { |
391 | werner | 386 | *p = qMin(*p, 1.f); |
387 | } |
||
388 | // set the parent cell to 1 |
||
1178 | werner | 389 | //kernel.valueAtIndex(kernel_offset, kernel_offset)=1.f; |
415 | werner | 390 | |
391 | |||
391 | werner | 392 | // some final statistics.... |
550 | werner | 393 | if (logLevelInfo()) qDebug() << "kernel setup.Species:"<< mSpecies->id() << "kernel-size: " << kernel.sizeX() << "x" << kernel.sizeY() << "pixels, sum: " << kernel.sum(); |
1176 | werner | 394 | |
395 | return fecundity_factor; |
||
391 | werner | 396 | } |
397 | |||
1176 | werner | 398 | void SeedDispersal::setupLDD() |
399 | { |
||
400 | mLDDDensity.clear(); mLDDDistance.clear(); |
||
401 | if (mKernelThresholdLDD >= mKernelThresholdArea) { |
||
402 | // no long distance dispersal |
||
403 | return; |
||
404 | |||
405 | } |
||
406 | double r_min = treemig_distanceTo(mKernelThresholdArea); |
||
407 | double r_max = treemig_distanceTo(mKernelThresholdLDD); |
||
408 | |||
409 | |||
410 | mLDDDistance.push_back(r_min); |
||
411 | for (int i=0;i<mLDDRings;++i) { |
||
412 | double r_in = mLDDDistance.last(); |
||
413 | mLDDDistance.push_back(mLDDDistance.last() + (r_max-r_min)/static_cast<float>(mLDDRings)); |
||
414 | double r_out = mLDDDistance.last(); |
||
415 | // calculate the value of the kernel for the middle of the ring |
||
416 | double ring_in = treemig(r_in); // kernel value at the inner border of the ring |
||
417 | double ring_out = treemig(r_out); // kernel value at the outer border of the ring |
||
1178 | werner | 418 | double ring_val = (ring_in + ring_out)/2.; // this is the average p |
1176 | werner | 419 | ring_val *= mFecundityFactor; // include fecundity (along the lines of the kernel calculations) |
420 | // calculate the area of the ring |
||
421 | double ring_area = (r_out*r_out - r_in*r_in)*M_PI; // in square meters |
||
422 | double ring_px = ring_area / (mSeedMap.cellsize()*mSeedMap.cellsize()); // # of seed pixels on the ring |
||
423 | double n_px = ring_val * ring_px / mLDDProbability; |
||
424 | |||
425 | mLDDDensity.push_back(n_px); |
||
426 | } |
||
427 | qDebug() << "Setup LDD for" << species()->name() << ", using probability: "<< mLDDProbability<< ": Distances:" << mLDDDistance << ", seed pixels:" << mLDDDensity; |
||
428 | |||
429 | |||
430 | } |
||
431 | |||
391 | werner | 432 | /* R-Code: |
930 | werner | 433 | treemig=function(as1,as2,ks,d) # two-part exponential function, cf. Lischke & Loeffler (2006), Annex |
391 | werner | 434 | { |
435 | p1=(1-ks)*exp(-d/as1)/as1 |
||
436 | if(as2>0){p2=ks*exp(-d/as2)/as2}else{p2=0} |
||
437 | p1+p2 |
||
438 | } |
||
439 | */ |
||
670 | werner | 440 | |
441 | /// the used kernel function |
||
442 | /// see also Appendix B of iland paper II (note the different variable names) |
||
1178 | werner | 443 | /// the function returns the seed density at a point with distance 'distance'. |
391 | werner | 444 | double SeedDispersal::treemig(const double &distance) |
445 | { |
||
446 | double p1 = (1.-mTM_ks)*exp(-distance/mTM_as1)/mTM_as1; |
||
447 | double p2 = 0.; |
||
448 | if (mTM_as2>0.) |
||
449 | p2 = mTM_ks*exp(-distance/mTM_as2)/mTM_as2; |
||
1178 | werner | 450 | double s = p1 + p2; |
451 | // 's' is the density for radius 'distance' - not for specific point with that distance. |
||
452 | // (i.e. the integral over the one-dimensional treemig function is 1, but if applied for 2d cells, the |
||
453 | // sum would be much larger as all seeds arriving at 'distance' would be arriving somewhere at the circle with radius 'distance') |
||
454 | // convert that to a density at a point, by dividing with the circumference at the circle with radius 'distance' |
||
455 | s = s / (2.*std::max(distance, 0.01)*M_PI); |
||
456 | |||
457 | return s; |
||
391 | werner | 458 | } |
459 | |||
1178 | werner | 460 | double SeedDispersal::treemig_centercell(const double &max_distance) |
461 | { |
||
462 | // use 100 steps and calculate dispersal kernel for consecutive rings |
||
463 | double sum = 0.; |
||
464 | for (int i=0;i<100;i++) { |
||
465 | double r_in = i*max_distance/100.; |
||
466 | double r_out = (i+1)*max_distance/100.; |
||
467 | double ring_area = (r_out*r_out-r_in*r_in)*M_PI; |
||
468 | // the value of each ring is: treemig(r) * area of the ring |
||
469 | sum += treemig((r_out+r_in)/2.)*ring_area; |
||
470 | } |
||
471 | return sum; |
||
472 | } |
||
473 | |||
391 | werner | 474 | /// calculate the distance where the probability falls below 'value' |
475 | double SeedDispersal::treemig_distanceTo(const double value) |
||
476 | { |
||
477 | double dist = 0.; |
||
478 | while (treemig(dist)>value && dist<10000.) |
||
479 | dist+=10; |
||
480 | return dist; |
||
481 | } |
||
482 | |||
764 | werner | 483 | void SeedDispersal::setupExternalSeedsForSpecies(Species *species) |
484 | { |
||
485 | if (!mExtSeedData.contains(species->id())) |
||
486 | return; // nothing to do |
||
487 | qDebug() << "setting up external seed map for" << species->id(); |
||
488 | QVector<double> &pcts = mExtSeedData[species->id()]; |
||
489 | mExternalSeedMap.setup(mSeedMap); |
||
490 | mExternalSeedMap.initialize(0.f); |
||
491 | for (int sector_x=0; sector_x<mExtSeedSizeX; ++sector_x) |
||
492 | for (int sector_y=0; sector_y<mExtSeedSizeY; ++sector_y) { |
||
493 | int xmin,xmax,ymin,ymax; |
||
494 | int fx = mExternalSeedMap.sizeX() / mExtSeedSizeX; // number of cells per sector |
||
495 | xmin = sector_x*fx; |
||
496 | xmax = (sector_x+1)*fx; |
||
497 | fx = mExternalSeedMap.sizeY() / mExtSeedSizeY; // number of cells per sector |
||
498 | ymin = sector_y*fx; |
||
499 | ymax = (sector_y+1)*fx; |
||
500 | // now loop over the whole sector |
||
501 | int index = sector_y*mExtSeedSizeX + sector_x; |
||
502 | double p = pcts[index]; |
||
503 | for (int y=ymin;y<ymax;++y) |
||
504 | for (int x=xmin;x<xmax;++x) { |
||
505 | // check |
||
506 | if (mExternalSeedBaseMap->valueAtIndex(x,y)==2.f) |
||
507 | if (drandom()<p) |
||
508 | mExternalSeedMap.valueAtIndex(x,y) = 1.f; // flag |
||
509 | } |
||
391 | werner | 510 | |
764 | werner | 511 | } |
512 | } |
||
513 | |||
514 | |||
391 | werner | 515 | // ************ Dispersal ************** |
516 | |||
517 | |||
375 | werner | 518 | /// debug function: loads a image of arbirtrary size... |
519 | void SeedDispersal::loadFromImage(const QString &fileName) |
||
520 | { |
||
521 | mSeedMap.clear(); |
||
522 | loadGridFromImage(fileName, mSeedMap); |
||
523 | for (float* p=mSeedMap.begin();p!=mSeedMap.end();++p) |
||
524 | *p = *p>0.8?1.f:0.f; |
||
525 | |||
526 | } |
||
527 | |||
528 | void SeedDispersal::clear() |
||
529 | { |
||
764 | werner | 530 | if (!mExternalSeedMap.isEmpty()) { |
531 | // we have a preprocessed initial value for the external seed map (see setupExternalSeeds() et al) |
||
532 | mSeedMap.copy(mExternalSeedMap); |
||
533 | return; |
||
534 | } |
||
535 | // clear the map |
||
1102 | werner | 536 | float background_value = static_cast<float>(mExternalSeedBackgroundInput); // there is potentitally a background probability <>0 for all pixels. |
836 | werner | 537 | mSeedMap.initialize(background_value); |
472 | werner | 538 | if (mHasExternalSeedInput) { |
539 | // if external seed input is enabled, the buffer area of the seed maps is |
||
540 | // "turned on", i.e. set to 1. |
||
1167 | werner | 541 | int buf_size = GlobalSettings::instance()->settings().valueInt("model.world.buffer",0.) / static_cast<int>(mSeedMap.cellsize()); |
491 | werner | 542 | // if a special buffer is defined, reduce the size of the input |
543 | if (mExternalSeedBuffer>0) |
||
544 | buf_size -= mExternalSeedBuffer; |
||
472 | werner | 545 | if (buf_size>0) { |
546 | int ix,iy; |
||
547 | for (iy=0;iy<mSeedMap.sizeY();++iy) |
||
548 | for (ix=0;ix<mSeedMap.sizeX(); ++ix) |
||
491 | werner | 549 | if (iy<buf_size || iy>=mSeedMap.sizeY()-buf_size || ix<buf_size || ix>=mSeedMap.sizeX()-buf_size) { |
550 | if (mExternalSeedDirection==0) { |
||
551 | // seeds from all directions |
||
552 | mSeedMap.valueAtIndex(ix,iy)=1.f; |
||
553 | } else { |
||
554 | // seeds only from specific directions |
||
555 | float value = 0.f; |
||
556 | if (isBitSet(mExternalSeedDirection,1) && ix>=mSeedMap.sizeX()-buf_size) value = 1; // north |
||
557 | if (isBitSet(mExternalSeedDirection,2) && iy<buf_size) value = 1; // east |
||
558 | if (isBitSet(mExternalSeedDirection,3) && ix<buf_size) value = 1; // south |
||
559 | if (isBitSet(mExternalSeedDirection,4) && iy>=mSeedMap.sizeY()-buf_size) value = 1; // west |
||
560 | mSeedMap.valueAtIndex(ix,iy)=value; |
||
561 | } |
||
562 | } |
||
472 | werner | 563 | } else { |
564 | qDebug() << "external seed input: Error: invalid buffer size???"; |
||
565 | } |
||
566 | } |
||
375 | werner | 567 | } |
568 | |||
1176 | werner | 569 | static int _debug_ldd=0; |
375 | werner | 570 | void SeedDispersal::execute() |
571 | { |
||
991 | werner | 572 | #ifdef ILAND_GUI |
573 | int year = GlobalSettings::instance()->currentYear(); |
||
574 | QString path; |
||
472 | werner | 575 | if (mDumpSeedMaps) { |
1102 | werner | 576 | path = GlobalSettings::instance()->path( GlobalSettings::instance()->settings().value("model.settings.seedDispersal.dumpSeedMapsPath") ); |
472 | werner | 577 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_before_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
618 | werner | 578 | qDebug() << "saved seed map image to" << path; |
992 | werner | 579 | } |
989 | werner | 580 | #else |
992 | werner | 581 | if (mDumpSeedMaps) |
989 | werner | 582 | qDebug() << "saving of seedmaps only supported in the iLand GUI."; |
583 | #endif |
||
992 | werner | 584 | |
1102 | werner | 585 | DebugTimer t("seed dispersal", true); |
391 | werner | 586 | { |
375 | werner | 587 | // (1) detect edges |
619 | werner | 588 | if (edgeDetection()) { |
1106 | werner | 589 | |
590 | #ifdef ILAND_GUI |
||
591 | if (mDumpSeedMaps) { |
||
592 | gridToImage(seedMap(), true, -1., 1.).save(QString("%1/seed_edge_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
593 | } |
||
594 | #endif |
||
595 | |||
596 | // (2) distribute seed probabilites from edges |
||
619 | werner | 597 | distribute(); |
391 | werner | 598 | } |
1167 | werner | 599 | |
600 | // special case serotiny |
||
601 | if (mHasPendingSerotiny) { |
||
1168 | werner | 602 | qDebug() << "calculating extra seed rain (serotiny)...."; |
603 | #ifdef ILAND_GUI |
||
604 | if (mDumpSeedMaps) { |
||
605 | gridToImage(mSeedMapSerotiny, true, 0., 1.).save(QString("%1/seed_serotiny_before_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
606 | } |
||
607 | #endif |
||
608 | if (edgeDetection(&mSeedMapSerotiny)) |
||
609 | distribute(&mSeedMapSerotiny); |
||
1167 | werner | 610 | // copy back data |
1168 | werner | 611 | float *sero=mSeedMapSerotiny.begin(); |
1167 | werner | 612 | for (float* p=mSeedMap.begin();p!=mSeedMap.end();++p, ++sero) |
613 | *p = std::max(*p, *sero); |
||
614 | |||
1168 | werner | 615 | float total = mSeedMapSerotiny.sum(); |
616 | #ifdef ILAND_GUI |
||
617 | if (mDumpSeedMaps) { |
||
618 | gridToImage(mSeedMapSerotiny, true, 0., 1.).save(QString("%1/seed_serotiny_after_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
||
619 | } |
||
620 | #endif |
||
621 | mSeedMapSerotiny.initialize(0.f); // clear |
||
1167 | werner | 622 | mHasPendingSerotiny = false; |
623 | qDebug() << "serotiny event: extra seed input" << total << "(total sum of seed probability over all pixels of the serotiny seed map) of species" << mSpecies->name(); |
||
619 | werner | 624 | } |
1167 | werner | 625 | } // debugtimer |
989 | werner | 626 | #ifdef ILAND_GUI |
1102 | werner | 627 | if (mDumpSeedMaps) { |
1168 | werner | 628 | qDebug() << "finished seed dispersal for species. time: " << mSpecies->id() << t.elapsed(); |
472 | werner | 629 | gridToImage(seedMap(), true, 0., 1.).save(QString("%1/seed_after_%2_%3.png").arg(path).arg(mSpecies->id()).arg(year)); |
1102 | werner | 630 | } |
1064 | werner | 631 | |
632 | if (!mDumpNextYearFileName.isEmpty()) { |
||
633 | Helper::saveToTextFile(GlobalSettings::instance()->path(mDumpNextYearFileName), gridToESRIRaster(seedMap())); |
||
634 | qDebug() << "saved seed map for " << species()->id() << "to" << GlobalSettings::instance()->path(mDumpNextYearFileName); |
||
635 | mDumpNextYearFileName = QString(); |
||
636 | } |
||
1176 | werner | 637 | qDebug() << "LDD-count:" << _debug_ldd; |
1064 | werner | 638 | |
989 | werner | 639 | #endif |
375 | werner | 640 | } |
641 | |||
373 | werner | 642 | /** scans the seed image and detects "edges". |
643 | edges are then subsequently marked (set to -1). This is pass 1 of the seed distribution process. |
||
644 | */ |
||
1167 | werner | 645 | bool SeedDispersal::edgeDetection(Grid<float> *seed_map) |
373 | werner | 646 | { |
647 | float *p_above, *p, *p_below; |
||
1167 | werner | 648 | Grid<float> &seedmap = seed_map ? *seed_map : mSeedMap; // switch to extra seed map if provided |
649 | int dy = seedmap.sizeY(); |
||
650 | int dx = seedmap.sizeX(); |
||
373 | werner | 651 | int x,y; |
619 | werner | 652 | bool found = false; |
1106 | werner | 653 | |
654 | // fill mini-gaps |
||
1168 | werner | 655 | int n_gaps_filled=0; |
373 | werner | 656 | for (y=1;y<dy-1;++y){ |
1167 | werner | 657 | p = seedmap.ptr(1,y); |
373 | werner | 658 | p_above = p - dx; // one line above |
659 | p_below = p + dx; // one line below |
||
660 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
1106 | werner | 661 | if (*p < 0.999f) { |
662 | |||
663 | if ((*(p_above-1)==1.f) + (*p_above==1.f) + (*(p_above+1)==1.f) + |
||
664 | (*(p-1)==1.f) + (*(p+1)==1.f) + |
||
1168 | werner | 665 | (*(p_below-1)==1.f) + (*p_below==1.f) + (*(p_below+1)==1.f) > 3) { |
1106 | werner | 666 | *p=0.999f; // if more than 3 neighbors are active pixels, the value is high |
1168 | werner | 667 | ++n_gaps_filled; |
668 | } |
||
1106 | werner | 669 | } |
670 | |||
671 | } |
||
672 | } |
||
673 | |||
674 | |||
675 | // now detect the edges |
||
676 | int n_edges=0 ; |
||
677 | for (y=1;y<dy-1;++y){ |
||
1168 | werner | 678 | p = seedmap.ptr(1,y); |
1106 | werner | 679 | p_above = p - dx; // one line above |
680 | p_below = p + dx; // one line below |
||
681 | for (x=1;x<dx-1;++x,++p,++p_below, ++p_above) { |
||
682 | if (*p == 1.f) { |
||
619 | werner | 683 | found = true; |
1106 | werner | 684 | if ( (*(p_above-1)<0.999f && *(p_above-1)>=0.f) |
685 | || (*p_above<0.999f && *p_above>=0.f) |
||
686 | || (*(p_above+1)<0.999f && *(p_above+1)>=0.f) |
||
687 | || (*(p-1)<0.999f && *(p-1)>=0.f) |
||
688 | || (*(p+1)<0.999f && (*p+1)>=0.f) |
||
689 | || (*(p_below-1)<0.999f && *(p_below-1)>=0.f) |
||
690 | || (*p_below<0.999f && *p_below>=0.f) |
||
691 | || (*(p_below+1)<0.999f && *(p_below+1)>=0.f ) ) { |
||
692 | *p=-1.f; // if any surrounding pixel is >=0 & <0.999: -> mark as edge |
||
693 | ++n_edges; |
||
694 | } |
||
373 | werner | 695 | } |
696 | |||
697 | } |
||
698 | } |
||
1168 | werner | 699 | if (mDumpSeedMaps) |
700 | qDebug() << "species:" << mSpecies->id() << "# of gaps filled: " << n_gaps_filled << "# of edge-pixels:" << n_edges; |
||
619 | werner | 701 | return found; |
373 | werner | 702 | } |
703 | |||
704 | /** do the seed probability distribution. |
||
705 | This is phase 2. Apply the seed kernel for each "edge" point identified in phase 1. |
||
706 | */ |
||
1167 | werner | 707 | void SeedDispersal::distribute(Grid<float> *seed_map) |
373 | werner | 708 | { |
709 | int x,y; |
||
1167 | werner | 710 | Grid<float> &seedmap = seed_map ? *seed_map : mSeedMap; // switch to extra seed map if provided |
711 | float *end = seedmap.end(); |
||
712 | float *p = seedmap.begin(); |
||
415 | werner | 713 | // choose the kernel depending whether there is a seed year for the current species or not |
1168 | werner | 714 | Grid<float> *kernel = species()->isSeedYear()? &mKernelSeedYear : &mKernelNonSeedYear; |
1167 | werner | 715 | // extra case: serotiny |
716 | if (seed_map) |
||
1168 | werner | 717 | kernel = &mKernelSerotiny; |
1167 | werner | 718 | |
1168 | werner | 719 | int offset = kernel->sizeX() / 2; // offset is the index of the center pixel |
373 | werner | 720 | for(;p!=end;++p) { |
375 | werner | 721 | if (*p==-1.f) { |
373 | werner | 722 | // edge pixel found. Now apply the kernel.... |
1167 | werner | 723 | QPoint pt=seedmap.indexOf(p); |
1168 | werner | 724 | for (y=-offset;y<=offset;++y) { |
725 | for (x=-offset;x<=offset;++x) { |
||
726 | float &kernel_value = kernel->valueAtIndex(x+offset, y+offset); |
||
727 | if (kernel_value>0.f && seedmap.isIndexValid(pt.x()+x, pt.y()+y)) { |
||
1167 | werner | 728 | float &val = seedmap.valueAtIndex(pt.x()+x, pt.y()+y); |
1106 | werner | 729 | if (val!=-1.f) |
1168 | werner | 730 | val = qMin(1.f - (1.f - val)*(1.f-kernel_value),1.f ); |
373 | werner | 731 | } |
1168 | werner | 732 | } |
733 | } |
||
1176 | werner | 734 | // long distance dispersal |
735 | if (!mLDDDensity.isEmpty()) { |
||
736 | double m = species()->isSeedYear() ? 1. : mNonSeedYearFraction; |
||
737 | for (int r=0;r<mLDDDensity.size(); ++r) { |
||
738 | float ldd_val = mLDDProbability; // pixels will have this probability |
||
739 | int n = round( mLDDDensity[r]*m ); // number of pixels to activate |
||
740 | for (int i=0;i<n;++i) { |
||
741 | // distance and direction: |
||
742 | double radius = nrandom(mLDDDistance[r], mLDDDistance[r+1]) / seedmap.cellsize(); // choose a random distance (in pixels) |
||
743 | double phi = drandom()*2.*M_PI; // choose a random direction |
||
744 | QPoint ldd(pt.x() + radius*cos(phi), pt.y() + radius*sin(phi)); |
||
745 | if (seedmap.isIndexValid(ldd)) { |
||
746 | float &val = seedmap.valueAtIndex(ldd); |
||
747 | _debug_ldd++; |
||
748 | // use the same adding of probabilities |
||
749 | if (val!=-1.f) |
||
750 | val = qMin(1.f - (1.f - val)*(1.f-ldd_val), 1.f); |
||
751 | } |
||
752 | } |
||
753 | } |
||
754 | } |
||
375 | werner | 755 | *p=1.f; // mark as processed |
373 | werner | 756 | } // *p==1 |
757 | } // for() |
||
758 | } |