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