WebSVN - public iLand - Diff - Rev 265 and 277


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
/** @class ResourceUnit
  ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
  (Climate, Soil, Water, ...).

  */
#include <QtCore>
#include "global.h"

#include "resourceunit.h"
#include "resourceunitspecies.h"
#include "speciesset.h"
#include "species.h"
#include "production3pg.h"
#include "model.h"
#include "climate.h"
#include "watercycle.h"
#include "helper.h"

ResourceUnit::~ResourceUnit()
{
    if (mWater)
        delete mWater;
}

ResourceUnit::ResourceUnit(const int index)
{
    mSpeciesSet = 0;
    mClimate = 0;
    mIndex = index;
    mWater = new WaterCycle();

    mTrees.reserve(100); // start with space for 100 trees.
}

void ResourceUnit::setup()
{
    mWater->setup(this);
}

/// set species and setup the species-per-RU-data
void ResourceUnit::setSpeciesSet(SpeciesSet *set)
{
    mSpeciesSet = set;
    mRUSpecies.clear();
    mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
    for (int i=0;i<set->count();i++) {
        Species *s = const_cast<Species*>(mSpeciesSet->species(i));
        if (!s)
            throw IException("ResourceUnit::setSpeciesSet: invalid index!");

        /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
           If the container memory is relocated (QVector), the pointer gets invalid!!!
           Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
        mRUSpecies[i].setup(s,this); // setup this element

    }
}

ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
{
    return mRUSpecies[species->index()];
}

Tree &ResourceUnit::newTree()
{
    // start simple: just append to the vector...
    mTrees.append(Tree());
    return mTrees.back();
}

/// remove dead trees from tree list
/// reduce size of vector if lots of space is free
/// tests showed that this way of cleanup is very fast,
/// because no memory allocations are performed (simple memmove())
/// when trees are moved.
void ResourceUnit::cleanTreeList()
{
    QVector<Tree>::iterator last=mTrees.end()-1;
    QVector<Tree>::iterator current = mTrees.begin();
    while (last>=current && (*last).isDead())
        --last;

    while (current<last) {
        if ((*current).isDead()) {
            *current = *last; // copy data!
            --last; //
            while (last>=current && (*last).isDead())
                --last;
        }
        ++current;
    }
    ++last; // last points now to the first dead tree

    // free ressources
    mTrees.erase(last, mTrees.end());
    if (mTrees.capacity()>100) {
        if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
            int target_size = mTrees.count()*2;
            qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
            mTrees.reserve(qMax(target_size, 100));
        }
    }
}

void ResourceUnit::newYear()
{
    mAggregatedWLA = 0.;
    mAggregatedLA = 0.;
    mAggregatedLR = 0.;
    mEffectiveArea = 0.;
    mPixelCount = mStockedPixelCount = 0;
    // clear statistics global and per species...
}



/** production() is the "stand-level" part of the biomass production (3PG).
    - The amount of radiation intercepted by the stand is calculated
    - The 3PG production for each species and ressource unit is invoked  */
void ResourceUnit::production()
{

    if (mAggregatedWLA==0 || mPixelCount==0) {
        // nothing to do...
        mStatistics.clear();
        return;
    }

    // the pixel counters are filled during the height-grid-calculations
    mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)

    // calculate the leaf area index (LAI)
    double LAI = mAggregatedLA / mStockedArea;
    // calculate the intercepted radiation fraction using the law of Beer Lambert
    const double k = Model::settings().lightExtinctionCoefficient;
    double interception_fraction = 1. - exp(-k * LAI);
    mEffectiveArea = mStockedArea * interception_fraction; // m2

    // calculate the total weighted leaf area on this RU:
    mLRI_modification = interception_fraction *  mStockedArea / mAggregatedWLA;
    if (mLRI_modification == 0.)
        qDebug() << "lri modifaction==0!";


    DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
            .arg(LAI)
            .arg(interception_fraction)
            .arg(mLRI_modification)
            .arg(mStockedArea);
    );
    // soil water model - this determines soil water contents needed for response calculations
    {
    DebugTimer tw("water:run");
    mWater->run();
    }

    // invoke species specific calculation (3PG)
    //QVector<ResourceUnitSpecies>::iterator i;
    ResourceUnitSpecies *i;
    QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();
    mStatistics.clear();
    for (i=mRUSpecies.begin(); i!=iend; ++i) {
        i->calculate();
        i->statistics().clear();
        i->statisticsDead().clear();
        qDebug() << "species" << (*i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea();
    }
}

void ResourceUnit::calculateInterceptedArea()
{
    if (mAggregatedLR==0) {
        mEffectiveArea_perWLA = 0.;
        return;
    }
    Q_ASSERT(mAggregatedLR>0.);
    mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
    qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR  << "eff.area./wla:" << mEffectiveArea_perWLA;
}

void ResourceUnit::yearEnd()
{
    // calculate statistics for all tree species of the ressource unit
    int c = mRUSpecies.count();
    for (int i=0;i<c; i++) {
 
 
        mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees
        mRUSpecies[i].updateGWL(); // get sum od dead trees
        mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl)
        mStatistics.add(mRUSpecies[i].statistics());
    }
    mStatistics.calculate(); // aggreagte on stand level
}

/// refresh of tree based statistics.
void ResourceUnit::createStandStatistics()
{
    // clear statistics (ru-level and ru-species level)
    mStatistics.clear();
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i].statistics().clear();
        mRUSpecies[i].statisticsDead().clear();
    }

    // add all trees to the statistics objects of the species
    foreach(const Tree &t, mTrees) {
        if (!t.isDead())
            resourceUnitSpecies(t.species()).statistics().add(&t, 0);
    }
    // summarize statistics for the whole resource unit
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i].statistics().calculate();
        mStatistics.add(mRUSpecies[i].statistics());
    }
}
 

Rev 265	Rev 277
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
2	/** @class ResourceUnit	2	/** @class ResourceUnit
3	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components	3	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
4	(Climate, Soil, Water, ...).	4	(Climate, Soil, Water, ...).
5		5
6	*/	6	*/
7	#include <QtCore>	7	#include <QtCore>
8	#include "global.h"	8	#include "global.h"
9		9
10	#include "resourceunit.h"	10	#include "resourceunit.h"
11	#include "resourceunitspecies.h"	11	#include "resourceunitspecies.h"
12	#include "speciesset.h"	12	#include "speciesset.h"
13	#include "species.h"	13	#include "species.h"
14	#include "production3pg.h"	14	#include "production3pg.h"
15	#include "model.h"	15	#include "model.h"
16	#include "climate.h"	16	#include "climate.h"
17	#include "watercycle.h"	17	#include "watercycle.h"
18	#include "helper.h"	18	#include "helper.h"
19		19
20	ResourceUnit::~ResourceUnit()	20	ResourceUnit::~ResourceUnit()
21	{	21	{
22	if (mWater)	22	if (mWater)
23	delete mWater;	23	delete mWater;
24	}	24	}
25		25
26	ResourceUnit::ResourceUnit(const int index)	26	ResourceUnit::ResourceUnit(const int index)
27	{	27	{
28	mSpeciesSet = 0;	28	mSpeciesSet = 0;
29	mClimate = 0;	29	mClimate = 0;
30	mIndex = index;	30	mIndex = index;
31	mWater = new WaterCycle();	31	mWater = new WaterCycle();
32		32
33	mTrees.reserve(100); // start with space for 100 trees.	33	mTrees.reserve(100); // start with space for 100 trees.
34	}	34	}
35		35
36	void ResourceUnit::setup()	36	void ResourceUnit::setup()
37	{	37	{
38	mWater->setup(this);	38	mWater->setup(this);
39	}	39	}
40		40
41	/// set species and setup the species-per-RU-data	41	/// set species and setup the species-per-RU-data
42	void ResourceUnit::setSpeciesSet(SpeciesSet *set)	42	void ResourceUnit::setSpeciesSet(SpeciesSet *set)
43	{	43	{
44	mSpeciesSet = set;	44	mSpeciesSet = set;
45	mRUSpecies.clear();	45	mRUSpecies.clear();
46	mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated	46	mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
47	for (int i=0;i<set->count();i++) {	47	for (int i=0;i<set->count();i++) {
48	Species s = const_cast<Species>(mSpeciesSet->species(i));	48	Species s = const_cast<Species>(mSpeciesSet->species(i));
49	if (!s)	49	if (!s)
50	throw IException("ResourceUnit::setSpeciesSet: invalid index!");	50	throw IException("ResourceUnit::setSpeciesSet: invalid index!");
51		51
52	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.	52	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
53	If the container memory is relocated (QVector), the pointer gets invalid!!!	53	If the container memory is relocated (QVector), the pointer gets invalid!!!
54	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */	54	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
55	mRUSpecies[i].setup(s,this); // setup this element	55	mRUSpecies[i].setup(s,this); // setup this element
-		56
56	}	57	}
57	}	58	}
58		59
59	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)	60	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
60	{	61	{
61	return mRUSpecies[species->index()];	62	return mRUSpecies[species->index()];
62	}	63	}
63		64
64	Tree &ResourceUnit::newTree()	65	Tree &ResourceUnit::newTree()
65	{	66	{
66	// start simple: just append to the vector...	67	// start simple: just append to the vector...
67	mTrees.append(Tree());	68	mTrees.append(Tree());
68	return mTrees.back();	69	return mTrees.back();
69	}	70	}
70		71
71	/// remove dead trees from tree list	72	/// remove dead trees from tree list
72	/// reduce size of vector if lots of space is free	73	/// reduce size of vector if lots of space is free
73	/// tests showed that this way of cleanup is very fast,	74	/// tests showed that this way of cleanup is very fast,
74	/// because no memory allocations are performed (simple memmove())	75	/// because no memory allocations are performed (simple memmove())
75	/// when trees are moved.	76	/// when trees are moved.
76	void ResourceUnit::cleanTreeList()	77	void ResourceUnit::cleanTreeList()
77	{	78	{
78	QVector<Tree>::iterator last=mTrees.end()-1;	79	QVector<Tree>::iterator last=mTrees.end()-1;
79	QVector<Tree>::iterator current = mTrees.begin();	80	QVector<Tree>::iterator current = mTrees.begin();
80	while (last>=current && (*last).isDead())	81	while (last>=current && (*last).isDead())
81	--last;	82	--last;
82		83
83	while (current<last) {	84	while (current<last) {
84	if ((*current).isDead()) {	85	if ((*current).isDead()) {
85	current = last; // copy data!	86	current = last; // copy data!
86	--last; //	87	--last; //
87	while (last>=current && (*last).isDead())	88	while (last>=current && (*last).isDead())
88	--last;	89	--last;
89	}	90	}
90	++current;	91	++current;
91	}	92	}
92	++last; // last points now to the first dead tree	93	++last; // last points now to the first dead tree
93		94
94	// free ressources	95	// free ressources
95	mTrees.erase(last, mTrees.end());	96	mTrees.erase(last, mTrees.end());
96	if (mTrees.capacity()>100) {	97	if (mTrees.capacity()>100) {
97	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {	98	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
98	int target_size = mTrees.count()*2;	99	int target_size = mTrees.count()*2;
99	qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;	100	qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
100	mTrees.reserve(qMax(target_size, 100));	101	mTrees.reserve(qMax(target_size, 100));
101	}	102	}
102	}	103	}
103	}	104	}
104		105
105	void ResourceUnit::newYear()	106	void ResourceUnit::newYear()
106	{	107	{
107	mAggregatedWLA = 0.;	108	mAggregatedWLA = 0.;
108	mAggregatedLA = 0.;	109	mAggregatedLA = 0.;
109	mAggregatedLR = 0.;	110	mAggregatedLR = 0.;
110	mEffectiveArea = 0.;	111	mEffectiveArea = 0.;
111	mPixelCount = mStockedPixelCount = 0;	112	mPixelCount = mStockedPixelCount = 0;
112	// clear statistics global and per species...	113	// clear statistics global and per species...
113	}	114	}
114		115
115		116
116		117
117	/** production() is the "stand-level" part of the biomass production (3PG).	118	/** production() is the "stand-level" part of the biomass production (3PG).
118	- The amount of radiation intercepted by the stand is calculated	119	- The amount of radiation intercepted by the stand is calculated
119	- The 3PG production for each species and ressource unit is invoked */	120	- The 3PG production for each species and ressource unit is invoked */
120	void ResourceUnit::production()	121	void ResourceUnit::production()
121	{	122	{
122		123
123	if (mAggregatedWLA==0 \|\| mPixelCount==0) {	124	if (mAggregatedWLA==0 \|\| mPixelCount==0) {
124	// nothing to do...	125	// nothing to do...
125	mStatistics.clear();	126	mStatistics.clear();
126	return;	127	return;
127	}	128	}
128		129
129	// the pixel counters are filled during the height-grid-calculations	130	// the pixel counters are filled during the height-grid-calculations
130	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)	131	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)
131		132
132	// calculate the leaf area index (LAI)	133	// calculate the leaf area index (LAI)
133	double LAI = mAggregatedLA / mStockedArea;	134	double LAI = mAggregatedLA / mStockedArea;
134	// calculate the intercepted radiation fraction using the law of Beer Lambert	135	// calculate the intercepted radiation fraction using the law of Beer Lambert
135	const double k = Model::settings().lightExtinctionCoefficient;	136	const double k = Model::settings().lightExtinctionCoefficient;
136	double interception_fraction = 1. - exp(-k * LAI);	137	double interception_fraction = 1. - exp(-k * LAI);
137	mEffectiveArea = mStockedArea * interception_fraction; // m2	138	mEffectiveArea = mStockedArea * interception_fraction; // m2
138		139
139	// calculate the total weighted leaf area on this RU:	140	// calculate the total weighted leaf area on this RU:
140	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA;	141	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA;
141	if (mLRI_modification == 0.)	142	if (mLRI_modification == 0.)
142	qDebug() << "lri modifaction==0!";	143	qDebug() << "lri modifaction==0!";
143		144
144		145
145	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")	146	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
146	.arg(LAI)	147	.arg(LAI)
147	.arg(interception_fraction)	148	.arg(interception_fraction)
148	.arg(mLRI_modification)	149	.arg(mLRI_modification)
149	.arg(mStockedArea);	150	.arg(mStockedArea);
150	);	151	);
151	// soil water model - this determines soil water contents needed for response calculations	152	// soil water model - this determines soil water contents needed for response calculations
152	{	153	{
153	DebugTimer tw("water:run");	154	DebugTimer tw("water:run");
154	mWater->run();	155	mWater->run();
155	}	156	}
156		157
157	// invoke species specific calculation (3PG)	158	// invoke species specific calculation (3PG)
158	//QVector<ResourceUnitSpecies>::iterator i;	159	//QVector<ResourceUnitSpecies>::iterator i;
159	ResourceUnitSpecies *i;	160	ResourceUnitSpecies *i;
160	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();	161	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();
161	mStatistics.clear();	162	mStatistics.clear();
162	for (i=mRUSpecies.begin(); i!=iend; ++i) {	163	for (i=mRUSpecies.begin(); i!=iend; ++i) {
163	i->calculate();	164	i->calculate();
164	i->statistics().clear();	165	i->statistics().clear();
165	i->statisticsDead().clear();	166	i->statisticsDead().clear();
166	qDebug() << "species" << (*i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea();	167	qDebug() << "species" << (*i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea();
167	}	168	}
168	}	169	}
169		170
170	void ResourceUnit::calculateInterceptedArea()	171	void ResourceUnit::calculateInterceptedArea()
171	{	172	{
172	if (mAggregatedLR==0) {	173	if (mAggregatedLR==0) {
173	mEffectiveArea_perWLA = 0.;	174	mEffectiveArea_perWLA = 0.;
174	return;	175	return;
175	}	176	}
176	Q_ASSERT(mAggregatedLR>0.);	177	Q_ASSERT(mAggregatedLR>0.);
177	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;	178	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
178	qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;	179	qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;
179	}	180	}
180		181
181	void ResourceUnit::yearEnd()	182	void ResourceUnit::yearEnd()
182	{	183	{
183	// calculate statistics for all tree species of the ressource unit	184	// calculate statistics for all tree species of the ressource unit
184	int c = mRUSpecies.count();	185	int c = mRUSpecies.count();
185	for (int i=0;i<c; i++) {	186	for (int i=0;i<c; i++) {
186	mRUSpecies[i].statistics().calculate();	-
187	mRUSpecies[i].statisticsDead().calculate();	-
-		187	mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees
-		188	mRUSpecies[i].updateGWL(); // get sum od dead trees
-		189	mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl)
188	mStatistics.add(mRUSpecies[i].statistics());	190	mStatistics.add(mRUSpecies[i].statistics());
189	}	191	}
190	mStatistics.calculate(); // aggreagte on stand level	192	mStatistics.calculate(); // aggreagte on stand level
191	}	193	}
192		194
193	/// refresh of tree based statistics.	195	/// refresh of tree based statistics.
194	void ResourceUnit::createStandStatistics()	196	void ResourceUnit::createStandStatistics()
195	{	197	{
196	// clear statistics (ru-level and ru-species level)	198	// clear statistics (ru-level and ru-species level)
197	mStatistics.clear();	199	mStatistics.clear();
198	for (int i=0;i<mRUSpecies.count();i++) {	200	for (int i=0;i<mRUSpecies.count();i++) {
199	mRUSpecies[i].statistics().clear();	201	mRUSpecies[i].statistics().clear();
200	mRUSpecies[i].statisticsDead().clear();	202	mRUSpecies[i].statisticsDead().clear();
201	}	203	}
202		204
203	// add all trees to the statistics objects of the species	205	// add all trees to the statistics objects of the species
204	foreach(const Tree &t, mTrees) {	206	foreach(const Tree &t, mTrees) {
205	if (!t.isDead())	207	if (!t.isDead())
206	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	208	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
207	}	209	}
208	// summarize statistics for the whole resource unit	210	// summarize statistics for the whole resource unit
209	for (int i=0;i<mRUSpecies.count();i++) {	211	for (int i=0;i<mRUSpecies.count();i++) {
210	mRUSpecies[i].statistics().calculate();	212	mRUSpecies[i].statistics().calculate();
211	mStatistics.add(mRUSpecies[i].statistics());	213	mStatistics.add(mRUSpecies[i].statistics());
212	}	214	}
213	}	215	}
214		216

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 265 → 277