[GRASS-SVN] r72769 - grass-addons/grass7/raster/r.spread.sod
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 4 18:17:32 PDT 2018
Author: wenzeslaus
Date: 2018-06-04 18:17:32 -0700 (Mon, 04 Jun 2018)
New Revision: 72769
Modified:
grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
grass-addons/grass7/raster/r.spread.sod/Img.cpp
grass-addons/grass7/raster/r.spread.sod/Img.h
grass-addons/grass7/raster/r.spread.sod/Spore.cpp
grass-addons/grass7/raster/r.spread.sod/Spore.h
grass-addons/grass7/raster/r.spread.sod/main.cpp
Log:
r.spread.sod: add mortality of infected trees
Now using branch mortality from upstream
(https://github.com/ncsu-landscape-dynamics/SOD-modeling-cpp).
Modified: grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md 2018-06-05 01:17:32 UTC (rev 72769)
@@ -4,8 +4,19 @@
The format is based on [Keep a Changelog](http://keepachangelog.com/).
-## 2017-09-05 - March 2018 update
+## 2018-03-26 - mortality update
+### Added
+
+- Mortality (Vaclav Petras)
+ - Enabled by a flag, mandatory mortality rate, optional start time
+ - Optional output series of accumulated dead tree counts per cell
+- Image class constructor taking another image using its dimensions
+ and a provided value (Vaclav Petras)
+- Multiply for image class is commutative (Vaclav Petras)
+
+## 2018-03-26 - March 2018 update
+
### Changed
- Changes for automatic compilation in GRASS GIS Addons (Vaclav Petras)
Modified: grass-addons/grass7/raster/r.spread.sod/Img.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.cpp 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/Img.cpp 2018-06-05 01:17:32 UTC (rev 72769)
@@ -69,6 +69,15 @@
std::copy(other.data, other.data + (width * height), data);
}
+Img::Img(const Img& other, int value)
+{
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = new int[width * height]{value};
+}
+
Img::Img(Img&& other)
{
width = other.width;
@@ -269,17 +278,12 @@
return out;
}
-Img Img::operator*(double factor) const
+Img Img::operator*(double value) const
{
- auto re_width = this->width;
- auto re_height = this->height;
- auto out = Img(re_width, re_height, this->w_e_res, this->n_s_res);
+ auto out = Img(width, height, w_e_res, n_s_res);
- for (int i = 0; i < re_height; i++) {
- for (int j = 0; j < re_width; j++) {
- out.data[i * width + j] = this->data[i * width + j] * factor;
- }
- }
+ std::transform(data, data + (width * height), out.data,
+ [&value](const int& a) { return a * value; });
return out;
}
@@ -348,6 +352,11 @@
return *this;
}
+Img operator*(double factor, const Img& image)
+{
+ return image * factor;
+}
+
Img::~Img()
{
if (data) {
Modified: grass-addons/grass7/raster/r.spread.sod/Img.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.h 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/Img.h 2018-06-05 01:17:32 UTC (rev 72769)
@@ -44,6 +44,8 @@
Img();
Img(Img&& other);
Img(const Img& other);
+ // Use dimensions of an existing object, but supply a new value
+ Img(const Img& other, int value);
//Img(int width,int height);
Img(const char *fileName);
Img(int width, int height, int w_e_res, int n_s_res);
@@ -101,7 +103,7 @@
Img operator-(const Img& image) const;
Img operator*(const Img& image) const;
Img operator/(const Img& image) const;
- Img operator*(double factor) const;
+ Img operator*(double value) const;
Img operator/(double value) const;
Img& operator+=(int value);
Img& operator-=(int value);
@@ -111,6 +113,9 @@
Img& operator-=(const Img& image);
Img& operator*=(const Img& image);
Img& operator/=(const Img& image);
+
+ friend Img operator*(double factor, const Img& image);
+
~Img();
void toGrassRaster(const char *name);
Modified: grass-addons/grass7/raster/r.spread.sod/Spore.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.cpp 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.cpp 2018-06-05 01:17:32 UTC (rev 72769)
@@ -121,7 +121,7 @@
}
}
-void Sporulation::SporeSpreadDisp_singleSpecies(Img& S, Img& I,
+void Sporulation::SporeSpreadDisp_singleSpecies(Img& S, Img& I, Img& I2,
const Img& lvtree_rast,
std::vector<std::tuple<int, int> >& outside_spores,
Rtype rtype, const double *weather,
@@ -193,6 +193,7 @@
prob_S *= weather_value;
if (U < prob_S) {
I(row, col) += 1;
+ I2(row, col) += 1;
S(row, col) -= 1;
}
}
Modified: grass-addons/grass7/raster/r.spread.sod/Spore.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.h 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.h 2018-06-05 01:17:32 UTC (rev 72769)
@@ -44,7 +44,7 @@
Sporulation(unsigned random_seed, const Img &size);
void SporeGen(const Img& I, const double *weather,
double weather_value, double rate);
- void SporeSpreadDisp_singleSpecies(Img& S, Img& I,
+ void SporeSpreadDisp_singleSpecies(Img& S, Img& I, Img& I2,
const Img& lvtree_rast, std::vector<std::tuple<int, int> > &outside_spores, Rtype rtype,
const double *weather, double weather_value,
double scale1, double kappa = 2,
Modified: grass-addons/grass7/raster/r.spread.sod/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/main.cpp 2018-06-04 20:01:24 UTC (rev 72768)
+++ grass-addons/grass7/raster/r.spread.sod/main.cpp 2018-06-05 01:17:32 UTC (rev 72769)
@@ -245,6 +245,8 @@
struct Option *start_time, *end_time, *seasonality;
struct Option *spore_rate, *wind;
struct Option *radial_type, *scale_1, *scale_2, *kappa, *gamma;
+ struct Option *infected_to_dead_rate, *first_year_to_die;
+ struct Option *dead_series;
struct Option *seed, *runs, *threads;
struct Option *output, *output_series;
struct Option *stddev, *stddev_series;
@@ -253,6 +255,7 @@
struct SodFlags
{
+ struct Flag *mortality;
struct Flag *generate_seed;
struct Flag *series_as_single_run;
};
@@ -452,6 +455,46 @@
opt.gamma->options = "0-1";
opt.gamma->guisection = _("Spores");
+ opt.infected_to_dead_rate = G_define_option();
+ opt.infected_to_dead_rate->type = TYPE_DOUBLE;
+ opt.infected_to_dead_rate->key = "mortality_rate";
+ opt.infected_to_dead_rate->label =
+ _("Mortality rate of infected trees");
+ opt.infected_to_dead_rate->description =
+ _("Percentage of infected trees that die in a given year"
+ " (trees are removed from the infected pool)");
+ opt.infected_to_dead_rate->options = "0-1";
+ opt.infected_to_dead_rate->guisection = _("Mortality");
+
+ opt.first_year_to_die = G_define_option();
+ opt.first_year_to_die->type = TYPE_INTEGER;
+ opt.first_year_to_die->key = "mortality_start";
+ opt.first_year_to_die->label =
+ _("Year of simulation when mortality is first applied");
+ opt.first_year_to_die->description =
+ _("How many years it takes for an infected tree to die"
+ " (value 1 for trees dying at the end of the first year)");
+ opt.first_year_to_die->guisection = _("Mortality");
+
+ opt.dead_series = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
+ opt.dead_series->key = "dead_series";
+ opt.dead_series->label =
+ _("Basename for series of number of dead trees");
+ opt.dead_series->description =
+ _("Basename for output series of number of dead trees"
+ " (requires mortality to be activated)");
+ opt.dead_series->required = NO;
+ opt.dead_series->guisection = _("Mortality");
+
+ flg.mortality = G_define_flag();
+ flg.mortality->key = 'm';
+ flg.mortality->label =
+ _("Apply mortality");
+ flg.mortality->description =
+ _("After a given amount of time, start removing dead trees"
+ " from the infected pool with a given rate");
+ flg.mortality->guisection = _("Mortality");
+
opt.seed = G_define_option();
opt.seed->key = "random_seed";
opt.seed->type = TYPE_INTEGER;
@@ -503,6 +546,15 @@
opt.weather_file, opt.weather_value, NULL);
G_option_collective(opt.moisture_file, opt.temperature_file, NULL);
+ // mortality
+ // flag and rate required always
+ // for simplicity of the code outputs allowed only with output
+ // for single run (avgs from runs likely not needed)
+ G_option_requires(flg.mortality, opt.infected_to_dead_rate, NULL);
+ G_option_requires(opt.first_year_to_die, flg.mortality, NULL);
+ G_option_requires_all(opt.dead_series, flg.mortality,
+ flg.series_as_single_run, NULL);
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -554,7 +606,29 @@
}
Date dd_start(start_time, 01, 01);
Date dd_end(end_time, 12, 31);
+ // difference in years (in dates) but including both years
+ auto num_years = dd_end.getYear() - dd_start.getYear() + 1;
+ // mortality
+ bool mortality = false;
+ unsigned first_year_to_die = 1; // starts at 1 (same as the opt)
+ double infected_to_dead_rate = 0.0;
+ if (flg.mortality->answer) {
+ mortality = true;
+ if (opt.first_year_to_die->answer) {
+ first_year_to_die = std::stoi(opt.first_year_to_die->answer);
+ if (first_year_to_die > num_years) {
+ G_fatal_error(
+ _("%s is too large (%d). It must be smaller or "
+ " equal than number of simulation years (%d)."),
+ opt.first_year_to_die->key,
+ first_year_to_die, num_years);
+ }
+ }
+ if (opt.infected_to_dead_rate->answer)
+ infected_to_dead_rate = std::stod(opt.infected_to_dead_rate->answer);
+ }
+
unsigned seed_value;
if (opt.seed->answer) {
seed_value = std::stoul(opt.seed->answer);
@@ -656,6 +730,19 @@
std::vector<Sporulation> sporulations;
std::vector<Img> sus_species_rasts(num_runs, S_species_rast);
std::vector<Img> inf_species_rasts(num_runs, I_species_rast);
+
+ // infected cohort for each year (index is cohort age)
+ // age starts with 0 (in year 1), 0 is oldest
+ std::vector<std::vector<Img> > inf_species_cohort_rasts(
+ num_runs, std::vector<Img>(num_years, Img(S_species_rast, 0)));
+
+ // we are using only the first dead img for visualization, but for
+ // parallelization we need all allocated anyway
+ std::vector<Img> dead_in_current_year(num_runs, Img(S_species_rast, 0));
+ // dead trees accumulated over years
+ // TODO: allow only when series as single run
+ Img accumulated_dead(Img(S_species_rast, 0));
+
sporulations.reserve(num_runs);
for (unsigned i = 0; i < num_runs; ++i)
sporulations.emplace_back(seed_value++, I_species_rast);
@@ -664,20 +751,22 @@
std::vector<unsigned> unresolved_weeks;
unresolved_weeks.reserve(max_weeks_in_year);
+ Date dd_current(dd_start);
+
// main simulation loop (weekly steps)
- for (int current_week = 0; ; current_week++, dd_start.increasedByWeek()) {
- if (dd_start < dd_end)
- if (!ss || !(dd_start.getMonth() > 9))
+ for (int current_week = 0; ; current_week++, dd_current.increasedByWeek()) {
+ if (dd_current < dd_end)
+ if (!ss || !(dd_current.getMonth() > 9))
unresolved_weeks.push_back(current_week);
// if all the oaks are infected, then exit
if (all_infected(S_species_rast)) {
- cerr << "In the " << dd_start << " all suspectible oaks are infected!" << endl;
+ cerr << "In the " << dd_current << " all suspectible oaks are infected!" << endl;
break;
}
// check whether the spore occurs in the month
- if (dd_start.isYearEnd() || dd_start >= dd_end) {
+ if (dd_current.isYearEnd() || dd_current >= dd_end) {
if (!unresolved_weeks.empty()) {
unsigned week_in_chunk = 0;
@@ -714,7 +803,9 @@
}
sporulations[run].SporeGen(inf_species_rasts[run], week_weather, weather_value, spore_rate);
+ auto current_age = dd_current.getYear() - dd_start.getYear();
sporulations[run].SporeSpreadDisp_singleSpecies(sus_species_rasts[run], inf_species_rasts[run],
+ inf_species_cohort_rasts[run][current_age],
lvtree_rast, outside_spores[run],
rtype, week_weather,
weather_value, scale1, kappa, pwdir, scale2,
@@ -724,6 +815,35 @@
}
unresolved_weeks.clear();
}
+ if (mortality && (dd_current.getYear() <= dd_end.getYear())) {
+ // to avoid problem with Jan 1 of the following year
+ // we explicitely check if we are in a valid year range
+ unsigned simulation_year = dd_current.getYear() - dd_start.getYear();
+ // only run to the current year of simulation
+ // (first year is 0):
+ // max index == sim year
+ // reduced by first time when trees start dying
+ // (counted from 1: first year == 1)
+ // e.g. for sim year 3, year dying 4, max index is 0
+ // max index = sim year - (dying year - 1)
+ // index is negative before we reach the year
+ // (so we can skip these years)
+ // sim year - (dying year - 1) < 0
+ // sim year < dying year - 1
+ if (simulation_year >= first_year_to_die - 1) {
+ auto max_index = simulation_year - (first_year_to_die - 1);
+ #pragma omp parallel for num_threads(threads)
+ for (unsigned run = 0; run < num_runs; run++) {
+ dead_in_current_year[run].zero();
+ for (unsigned age = 0; age <= max_index; age++) {
+ Img dead_in_cohort = infected_to_dead_rate * inf_species_cohort_rasts[run][age];
+ inf_species_cohort_rasts[run][age] -= dead_in_cohort;
+ dead_in_current_year[run] += dead_in_cohort;
+ }
+ inf_species_rasts[run] -= dead_in_current_year[run];
+ }
+ }
+ }
if ((opt.output_series->answer && !flg.series_as_single_run->answer)
|| opt.stddev_series->answer) {
// aggregate in the series
@@ -735,7 +855,7 @@
if (opt.output_series->answer) {
// write result
// date is always end of the year, even for seasonal spread
- string name = generate_name(opt.output_series->answer, dd_start);
+ string name = generate_name(opt.output_series->answer, dd_current);
if (flg.series_as_single_run->answer)
inf_species_rasts[0].toGrassRaster(name.c_str());
else
@@ -750,12 +870,19 @@
}
stddev /= num_runs;
stddev.for_each([](int& a){a = std::sqrt(a);});
- string name = generate_name(opt.stddev_series->answer, dd_start);
+ string name = generate_name(opt.stddev_series->answer, dd_current);
stddev.toGrassRaster(name.c_str());
}
+ if (mortality && opt.dead_series->answer) {
+ accumulated_dead += dead_in_current_year[0];
+ if (opt.dead_series->answer) {
+ string name = generate_name(opt.dead_series->answer, dd_current);
+ accumulated_dead.toGrassRaster(name.c_str());
+ }
+ }
}
- if (dd_start >= dd_end)
+ if (dd_current >= dd_end)
break;
}
More information about the grass-commit
mailing list