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