Subversion Repositories public iLand

Rev

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/grasscover.cpp':
1
Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/grasscover.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
#include "grasscover.h"
20
#include "grasscover.h"
21
21
22
#include "globalsettings.h"
22
#include "globalsettings.h"
23
#include "debugtimer.h"
23
#include "debugtimer.h"
24
#include "xmlhelper.h"
24
#include "xmlhelper.h"
25
#include "model.h"
25
#include "model.h"
26
#include "modelcontroller.h"
26
#include "modelcontroller.h"
27
27
28
const int GrassCover::GRASSCOVERSTEPS;
28
const int GrassCover::GRASSCOVERSTEPS;
29
29
30
GrassCover::GrassCover()
30
GrassCover::GrassCover()
31
{
31
{
32
    mLayers = new GrassCoverLayers();
32
    mLayers = new GrassCoverLayers();
33
    mLayers->setGrid(mGrid, this);
33
    mLayers->setGrid(mGrid, this);
34
    mEnabled = false;
34
    mEnabled = false;
35
    mType = Invalid;
35
    mType = Invalid;
36
}
36
}
37
37
38
GrassCover::~GrassCover()
38
GrassCover::~GrassCover()
39
{
39
{
40
    delete mLayers;
40
    delete mLayers;
41
}
41
}
42
42
43
void GrassCover::setup()
43
void GrassCover::setup()
44
{
44
{
45
    XmlHelper xml=GlobalSettings::instance()->settings();
45
    XmlHelper xml=GlobalSettings::instance()->settings();
46
    if (!xml.valueBool("model.settings.grass.enabled")) {
46
    if (!xml.valueBool("model.settings.grass.enabled")) {
47
        // clear grid
47
        // clear grid
48
        mGrid.clear();
48
        mGrid.clear();
49
        GlobalSettings::instance()->controller()->removeLayers(mLayers);
49
        GlobalSettings::instance()->controller()->removeLayers(mLayers);
50
        mEnabled=false;
50
        mEnabled=false;
51
        qDebug() << "grass module not enabled";
51
        qDebug() << "grass module not enabled";
52
        return;
52
        return;
53
    }
53
    }
54
    // create the grid
54
    // create the grid
55
    mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize());
55
    mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize());
56
    mGrid.wipe();
56
    mGrid.wipe();
57
    // mask out out-of-project areas
57
    // mask out out-of-project areas
58
    HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();
58
    HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();
59
    for (int i=0;i<mGrid.count();++i)
59
    for (int i=0;i<mGrid.count();++i)
60
        if (!hg->valueAtIndex(mGrid.index5(i)).isValid())
60
        if (!hg->valueAtIndex(mGrid.index5(i)).isValid())
61
            mGrid[i] = -1;
61
            mGrid[i] = -1;
62
62
63
    mType = Invalid;
63
    mType = Invalid;
64
    QString type = xml.value("model.settings.grass.type");
64
    QString type = xml.value("model.settings.grass.type");
65
    if (type == QStringLiteral("pixel"))
65
    if (type == QStringLiteral("pixel"))
66
        mType = Pixel;
66
        mType = Pixel;
67
67
68
    if (type==QStringLiteral("continuous"))
68
    if (type==QStringLiteral("continuous"))
69
        mType = Continuous;
69
        mType = Continuous;
70
70
71
    if (mType == Invalid)
71
    if (mType == Invalid)
72
        throw IException("GrassCover::setup: invalid 'grass.type'. Allowed: 'continous' and 'pixel'.");
72
        throw IException("GrassCover::setup: invalid 'grass.type'. Allowed: 'continous' and 'pixel'.");
73
73
74
    if (mType == Pixel) {
74
    if (mType == Pixel) {
75
        // setup of pixel based / discrete approach
75
        // setup of pixel based / discrete approach
76
        QString formula = xml.value("model.settings.grass.grassDuration");
76
        QString formula = xml.value("model.settings.grass.grassDuration");
77
        if (formula.isEmpty())
77
        if (formula.isEmpty())
78
            throw IException("GrassCover::setup(): missing equation for 'grassDuration'.");
78
            throw IException("GrassCover::setup(): missing equation for 'grassDuration'.");
79
        mPDF.setup(formula, 0., 100.);
79
        mPDF.setup(formula, 0., 100.);
80
        //mGrassEffect.setExpression(formula);
80
        //mGrassEffect.setExpression(formula);
81
81
82
        mGrassLIFThreshold = static_cast<float>( xml.valueDouble("model.settings.grass.LIFThreshold", 0.2) );
82
        mGrassLIFThreshold = static_cast<float>( xml.valueDouble("model.settings.grass.LIFThreshold", 0.2) );
83
83
84
        // clear array
84
        // clear array
85
        for (int i=0;i<GRASSCOVERSTEPS;++i) {
85
        for (int i=0;i<GRASSCOVERSTEPS;++i) {
86
            mEffect[i] = 0.;
86
            mEffect[i] = 0.;
87
        }
87
        }
88
88
89
    } else {
89
    } else {
90
90
91
        // setup of continuous grass concept
91
        // setup of continuous grass concept
92
92
93
        QString formula = xml.value("model.settings.grass.grassPotential");
93
        QString formula = xml.value("model.settings.grass.grassPotential");
94
        if (formula.isEmpty())
94
        if (formula.isEmpty())
95
            throw IException("setup of 'grass': required expression 'grassPotential' is missing.");
95
            throw IException("setup of 'grass': required expression 'grassPotential' is missing.");
96
        mGrassPotential.setExpression(formula);
96
        mGrassPotential.setExpression(formula);
97
        mGrassPotential.linearize(0.,1., qMin(GRASSCOVERSTEPS, 1000));
97
        mGrassPotential.linearize(0.,1., qMin(GRASSCOVERSTEPS, 1000));
98
98
99
        formula = xml.value("model.settings.grass.grassEffect");
99
        formula = xml.value("model.settings.grass.grassEffect");
100
        if (formula.isEmpty())
100
        if (formula.isEmpty())
101
            throw IException("setup of 'grass': required expression 'grassEffect' is missing.");
101
            throw IException("setup of 'grass': required expression 'grassEffect' is missing.");
102
        mGrassEffect.setExpression(formula);
102
        mGrassEffect.setExpression(formula);
103
        mMaxTimeLag = static_cast<int>( xml.valueDouble("model.settings.grass.maxTimeLag") );
103
        mMaxTimeLag = static_cast<int>( xml.valueDouble("model.settings.grass.maxTimeLag") );
104
        if (mMaxTimeLag==0)
104
        if (mMaxTimeLag==0)
105
            throw IException("setup of 'grass': value of 'maxTimeLag' is invalid or missing.");
105
            throw IException("setup of 'grass': value of 'maxTimeLag' is invalid or missing.");
106
        mGrowthRate = GRASSCOVERSTEPS / mMaxTimeLag;
106
        mGrowthRate = GRASSCOVERSTEPS / mMaxTimeLag;
107
107
108
        // set up the effect on regeneration in NSTEPS steps
108
        // set up the effect on regeneration in NSTEPS steps
109
        for (int i=0;i<GRASSCOVERSTEPS;++i) {
109
        for (int i=0;i<GRASSCOVERSTEPS;++i) {
110
            double effect = mGrassEffect.calculate(i/double(GRASSCOVERSTEPS-1));
110
            double effect = mGrassEffect.calculate(i/double(GRASSCOVERSTEPS-1));
111
            mEffect[i] = limit(effect, 0., 1.);
111
            mEffect[i] = limit(effect, 0., 1.);
112
        }
112
        }
113
113
114
        mMaxState = static_cast<qint16>( limit(mGrassPotential.calculate(1.f), 0., 1.)*(GRASSCOVERSTEPS-1)  ); // the max value of the potential function
114
        mMaxState = static_cast<qint16>( limit(mGrassPotential.calculate(1.f), 0., 1.)*(GRASSCOVERSTEPS-1)  ); // the max value of the potential function
115
    }
115
    }
116
116
117
    GlobalSettings::instance()->controller()->addLayers(mLayers, QStringLiteral("grass cover"));
117
    GlobalSettings::instance()->controller()->addLayers(mLayers, QStringLiteral("grass cover"));
118
    mEnabled = true;
118
    mEnabled = true;
119
    qDebug() << "setup of grass cover complete.";
119
    qDebug() << "setup of grass cover complete.";
120
120
121
}
121
}
122
122
123
void GrassCover::setInitialValues(const QVector<float *> &LIFpixels, const int percent)
123
void GrassCover::setInitialValues(const QVector<float *> &LIFpixels, const int percent)
124
{
124
{
125
    if (!enabled())
125
    if (!enabled())
126
        return;
126
        return;
127
    if (mType == Continuous) {
127
    if (mType == Continuous) {
128
        grass_grid_type cval = static_cast<grass_grid_type>( limit(percent / 100., 0., 1.)*(GRASSCOVERSTEPS-1) );
128
        grass_grid_type cval = static_cast<grass_grid_type>( limit(percent / 100., 0., 1.)*(GRASSCOVERSTEPS-1) );
129
        if (cval > mMaxState)
129
        if (cval > mMaxState)
130
            cval = mMaxState;
130
            cval = mMaxState;
131
131
132
        Grid<float> *lif_grid = GlobalSettings::instance()->model()->grid();
132
        Grid<float> *lif_grid = GlobalSettings::instance()->model()->grid();
133
        for (QVector<float *>::const_iterator it = LIFpixels.constBegin(); it!=LIFpixels.constEnd(); ++it)
133
        for (QVector<float *>::const_iterator it = LIFpixels.constBegin(); it!=LIFpixels.constEnd(); ++it)
134
            mGrid.valueAtIndex(lif_grid->indexOf(*it)) = cval;
134
            mGrid.valueAtIndex(lif_grid->indexOf(*it)) = cval;
135
    } else {
135
    } else {
136
        // mType == Pixel
136
        // mType == Pixel
137
        Grid<float> *lif_grid = GlobalSettings::instance()->model()->grid();
137
        Grid<float> *lif_grid = GlobalSettings::instance()->model()->grid();
138
        for (QVector<float *>::const_iterator it = LIFpixels.constBegin(); it!=LIFpixels.constEnd(); ++it) {
138
        for (QVector<float *>::const_iterator it = LIFpixels.constBegin(); it!=LIFpixels.constEnd(); ++it) {
139
            if (percent > irandom(0,100))
139
            if (percent > irandom(0,100))
140
                mGrid.valueAtIndex(lif_grid->indexOf(*it)) = static_cast<qint16>( mPDF.get() );
140
                mGrid.valueAtIndex(lif_grid->indexOf(*it)) = static_cast<qint16>( mPDF.get() );
141
            else
141
            else
142
                mGrid.valueAtIndex(lif_grid->indexOf(*it)) = 0;
142
                mGrid.valueAtIndex(lif_grid->indexOf(*it)) = 0;
143
        }
143
        }
144
144
145
    }
145
    }
146
}
146
}
147
147
148
void GrassCover::execute()
148
void GrassCover::execute()
149
{
149
{
150
    if (!enabled())
150
    if (!enabled())
151
        return;
151
        return;
152
152
153
    DebugTimer t("GrassCover");
153
    DebugTimer t("GrassCover");
154
154
155
    // Main function of the grass submodule
155
    // Main function of the grass submodule
156
    float *lif = GlobalSettings::instance()->model()->grid()->begin();
156
    float *lif = GlobalSettings::instance()->model()->grid()->begin();
157
    float *end_lif = GlobalSettings::instance()->model()->grid()->end();
157
    float *end_lif = GlobalSettings::instance()->model()->grid()->end();
158
    grass_grid_type *gr = mGrid.begin();
158
    grass_grid_type *gr = mGrid.begin();
159
159
160
    if (mType == Continuous) {
160
    if (mType == Continuous) {
161
        // loop over every LIF pixel
161
        // loop over every LIF pixel
162
        int skipped=0;
162
        int skipped=0;
163
        for (; lif!=end_lif;++lif, ++gr) {
163
        for (; lif!=end_lif;++lif, ++gr) {
164
            // calculate potential grass vegetation cover
164
            // calculate potential grass vegetation cover
165
            if (*lif == 1.f && *gr==mMaxState) {
165
            if (*lif == 1.f && *gr==mMaxState) {
166
                ++skipped;
166
                ++skipped;
167
                continue;
167
                continue;
168
            }
168
            }
169
169
170
            int potential = static_cast<int>( limit(mGrassPotential.calculate(*lif), 0., 1.)*(GRASSCOVERSTEPS-1) );
170
            int potential = static_cast<int>( limit(mGrassPotential.calculate(*lif), 0., 1.)*(GRASSCOVERSTEPS-1) );
171
            *gr = static_cast<qint16>( qMin( int(*gr) + mGrowthRate, potential) );
171
            *gr = static_cast<qint16>( qMin( int(*gr) + mGrowthRate, potential) );
172
172
173
        }
173
        }
174
        //qDebug() << "skipped" << skipped;
174
        //qDebug() << "skipped" << skipped;
175
    } else {
175
    } else {
176
        // type = Pixel
176
        // type = Pixel
177
        for (; lif!=end_lif;++lif, ++gr) {
177
        for (; lif!=end_lif;++lif, ++gr) {
178
            if (*gr<0)
178
            if (*gr<0)
179
                continue;
179
                continue;
180
            if (*gr>1)
180
            if (*gr>1)
181
                (*gr)--; // count down the years (until gr=1)
181
                (*gr)--; // count down the years (until gr=1)
182
182
183
            if (*gr==0 && *lif>mGrassLIFThreshold) {
183
            if (*gr==0 && *lif>mGrassLIFThreshold) {
184
                // enable grass cover
184
                // enable grass cover
185
                qint16 v = static_cast<qint16>( qMax(mPDF.get(), 0.) );
185
                qint16 v = static_cast<qint16>( qMax(mPDF.get(), 0.) );
186
                *gr = v + 1; // switch on...
186
                *gr = v + 1; // switch on...
187
            }
187
            }
188
            if (*gr==1 && *lif<mGrassLIFThreshold) {
188
            if (*gr==1 && *lif<mGrassLIFThreshold) {
189
                // now LIF is below the threshold - this enables the pixel get grassy again
189
                // now LIF is below the threshold - this enables the pixel get grassy again
190
                *gr = 0;
190
                *gr = 0;
191
            }
191
            }
192
        }
192
        }
193
    }
193
    }
194
194
195
}
195
}
196
196
197
197
198
198
199
double GrassCoverLayers::value(const grass_grid_type &data, const int index) const
199
double GrassCoverLayers::value(const grass_grid_type &data, const int index) const
200
{
200
{
201
    if (!mGrassCover->enabled()) return 0.;
201
    if (!mGrassCover->enabled()) return 0.;
202
    switch(index){
202
    switch(index){
203
    case 0: return mGrassCover->effect(data); //effect
203
    case 0: return mGrassCover->effect(data); //effect
204
    case 1: return mGrassCover->cover(data); // cover
204
    case 1: return mGrassCover->cover(data); // cover
205
    default: throw IException(QString("invalid variable index for a GrassCoverLayers: %1").arg(index));
205
    default: throw IException(QString("invalid variable index for a GrassCoverLayers: %1").arg(index));
206
    }
206
    }
207
}
207
}
208
208
209
const QVector<LayeredGridBase::LayerElement> &GrassCoverLayers::names()
209
const QVector<LayeredGridBase::LayerElement> &GrassCoverLayers::names()
210
{
210
{
211
    if (mNames.isEmpty())
211
    if (mNames.isEmpty())
212
        mNames = QVector<LayeredGridBase::LayerElement>()
212
        mNames = QVector<LayeredGridBase::LayerElement>()
213
                << LayeredGridBase::LayerElement(QLatin1Literal("effect"), QLatin1Literal("prohibiting effect on regeneration [0..1]"), GridViewGreens)
213
                << LayeredGridBase::LayerElement(QLatin1Literal("effect"), QLatin1Literal("prohibiting effect on regeneration [0..1]"), GridViewGreens)
214
                << LayeredGridBase::LayerElement(QLatin1Literal("cover"), QLatin1Literal("current grass cover on pixels [0..1 for continuous, or #(years+2) for pixel mode]"), GridViewGreens);
214
                << LayeredGridBase::LayerElement(QLatin1Literal("cover"), QLatin1Literal("current grass cover on pixels [0..1 for continuous, or #(years+2) for pixel mode]"), GridViewGreens);
215
    return mNames;
215
    return mNames;
216
216
217
}
217
}
218
 
218