[GRASS-SVN] r72867 - grass-addons/grass7/raster/r.spread.sod
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 21 13:09:01 PDT 2018
Author: wenzeslaus
Date: 2018-06-21 13:09:01 -0700 (Thu, 21 Jun 2018)
New Revision: 72867
Added:
grass-addons/grass7/raster/r.spread.sod/raster.h
Removed:
grass-addons/grass7/raster/r.spread.sod/Img.cpp
grass-addons/grass7/raster/r.spread.sod/Img.h
Modified:
grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
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/date.h
grass-addons/grass7/raster/r.spread.sod/main.cpp
Log:
r.spread.sod: add critical temperature
Now using branch lanternfly from upstream
(https://github.com/ncsu-landscape-dynamics/SOD-modeling-cpp).
This also significantly improves the raster (image) handling class
which is now a more general template class.
Modified: grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md 2018-06-21 20:09:01 UTC (rev 72867)
@@ -4,10 +4,50 @@
The format is based on [Keep a Changelog](http://keepachangelog.com/).
-## 2018-03-26 - mortality update
+## 2018-06-21 - Critical Temperature
### Added
+- Time-series of temperature raster maps specified as a text file.
+- Temperature rasters are used at a specified month to check against
+ a provided critical temperature and if the condition is met,
+ infected trees become susceptible again.
+
+### Changed
+
+- Img was replaced by a generalized Raster template class which can
+ handle both integers and floating point numbers. (Vaclav Petras)
+
+## 2018-06-13 - Spotted Lanternfly
+
+### Added
+
+- The date class now supports also month increments. (Vaclav Petras)
+- A new step option allows user to choose between weekly and monthly
+ increments in simulation. (Vaclav Petras)
+- A custom season can now be selected by the user. (Vaclav Petras)
+- A new test executable for the date class added and available in an
+ alternative Makefile. (Vaclav Petras)
+
+### Changed
+
+- The season option is no longer yes or no but a range of months.
+ (Vaclav Petras)
+
+### Fixed
+
+- Avoid segmentation fault by using the weather coefficients only when
+ available. (Vaclav Petras)
+- Make the NetCDF time-series input option always present regardless
+ compilation settings which avoids use of uninitialized variable later
+ (and thus undefined behavior). Modules now produces error with
+ explanation when option is used but it was compiled without NetCDF
+ support. (Vaclav Petras)
+
+## 2018-06-04 - Mortality Addition
+
+### Added
+
- Mortality (Vaclav Petras)
- Enabled by a flag, mandatory mortality rate, optional start time
- Optional output series of accumulated dead tree counts per cell
@@ -15,7 +55,7 @@
and a provided value (Vaclav Petras)
- Multiply for image class is commutative (Vaclav Petras)
-## 2018-03-26 - March 2018 update
+## 2018-03-26 - March 2018 Update
### Changed
@@ -25,7 +65,7 @@
- Explicitly include necessary standard C++ headers
- Add formalities: basic documentation and a proper GRASS module name
-## 2017-09-05 - September 2017 update
+## 2017-09-05 - September 2017 Update
### Added
@@ -52,7 +92,7 @@
- Interface now checks if only one way of providing weather coefficients
was used. (Vaclav Petras)
-## 2017-01-28 - January 2017 status
+## 2017-01-28 - January 2017 Status
### Added
Deleted: grass-addons/grass7/raster/r.spread.sod/Img.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.cpp 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/Img.cpp 2018-06-21 20:09:01 UTC (rev 72867)
@@ -1,413 +0,0 @@
-/*
- * SOD model - raster manipulation
- *
- * Copyright (C) 2015-2017 by the authors.
- *
- * Authors: Zexi Chen (zchen22 ncsu edu)
- * Vaclav Petras (wenzeslaus gmail com)
- *
- * The code contained herein is licensed under the GNU General Public
- * License. You may obtain a copy of the GNU General Public License
- * Version 2 or later at the following locations:
- *
- * http://www.opensource.org/licenses/gpl-license.html
- * http://www.gnu.org/copyleft/gpl.html
- */
-
-
-// define to support NetCDF format directly
-// (requires linking to netcdf_c++)
-// #define SOD_GDAL_SUPPORT
-
-#include "Img.h"
-
-extern "C" {
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include <grass/raster.h>
-}
-
-#ifdef SOD_GDAL_SUPPORT
-#include <gdal/gdal.h>
-#include <gdal/gdal_priv.h>
-#endif
-
-#include <algorithm>
-#include <stdexcept>
-
-using std::string;
-using std::cerr;
-using std::endl;
-
-/* Iterate over two ranges and apply a binary function which modifies
- * the first parameter.
- */
-template<class InputIt1, class InputIt2, class BinaryOperation>
-BinaryOperation for_each_zip(InputIt1 first1, InputIt1 last1, InputIt2 first2, BinaryOperation f) {
- for (; first1 != last1; ++first1, ++first2) {
- f(*first1, *first2);
- }
- return f;
-}
-
-Img::Img()
-{
- width = 0;
- height = 0;
- w_e_res = 0;
- n_s_res = 0;
- data = NULL;
-}
-
-Img::Img(const Img& other)
-{
- 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];
- 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;
- height = other.height;
- w_e_res = other.w_e_res;
- n_s_res = other.n_s_res;
- data = other.data;
- other.data = nullptr;
-}
-
-Img::Img(int width, int height, int w_e_res, int n_s_res)
-{
- this->width = width;
- this->height = height;
- this->w_e_res = w_e_res;
- this->n_s_res = n_s_res;
- this->data = new int[width * height];
-}
-
-Img::Img(int width, int height, int w_e_res, int n_s_res, int value)
-{
- this->width = width;
- this->height = height;
- this->w_e_res = w_e_res;
- this->n_s_res = n_s_res;
- this->data = new int[width * height]{value};
-}
-
-Img::Img(const char *fileName)
-{
-#ifdef SOD_GDAL_SUPPORT
- GDALDataset *dataset;
- GDALRasterBand *dataBand;
-
- GDALAllRegister();
- dataset = (GDALDataset *) GDALOpen(fileName, GA_ReadOnly);
- double adfGeoTransform[6];
-
- if (!dataset) {
- cerr << "Can not open the image!" << endl;
- }
- else {
- width = dataset->GetRasterXSize();
- height = dataset->GetRasterYSize();
-
- if (dataset->GetGeoTransform(adfGeoTransform) == CE_None) {
- w_e_res = abs(adfGeoTransform[1]);
- n_s_res = abs(adfGeoTransform[5]);
- }
-
- //cout << width << "x" << height <<endl;
- //cout << w_e_res << "X" << n_s_res << endl;
-
- dataBand = dataset->GetRasterBand(1);
- data = new int[width * height];
-
- CPLErr error = dataBand->RasterIO(GF_Read, 0, 0, width, height,
- data, width, height,
- GDT_Int32, 0, 0);
- if (error == CE_Failure)
- throw std::runtime_error(string("Reading raster failed"
- " in GDAL RasterIO: ")
- + CPLGetLastErrorMsg());
- GDALClose((GDALDatasetH) dataset);
- }
-#else
- throw std::runtime_error("GDAL support not available");
-#endif
-}
-
-
-// TODO: add move constuctor
-Img Img::fromGrassRaster(const char *name)
-{
- int fd = Rast_open_old(name, "");
-
- Img img;
-
- img.width = Rast_window_cols();
- img.height = Rast_window_rows();
-
- Cell_head region;
- Rast_get_window(®ion);
- img.w_e_res = region.ew_res;
- img.n_s_res = region.ns_res;
-
- img.data = new int[img.height * img.width];
-
- for (int row = 0; row < img.height; row++) {
- Rast_get_c_row(fd, img.data + (row * img.width), row);
- }
-
- Rast_close(fd);
- return img;
-}
-
-/*
- const int ** Img::getData(){
- return this->data;
- }
- */
-
-Img& Img::operator=(const Img& other)
-{
- if (this != &other)
- {
- if (data)
- delete[] data;
- 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];
- std::copy(other.data, other.data + (width * height), data);
- }
- return *this;
-}
-
-Img& Img::operator=(Img&& other)
-{
- if (this != &other)
- {
- if (data)
- delete[] data;
- width = other.width;
- height = other.height;
- w_e_res = other.w_e_res;
- n_s_res = other.n_s_res;
- data = other.data;
- other.data = nullptr;
- }
- return *this;
-}
-
-Img Img::operator+(const Img& image) const
-{
- if (this->width != image.getWidth() || this->height != image.getHeight()) {
- cerr << "The height or width of one image do not match with that of the other one!" << endl;
- return Img();
- }
- else {
- 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);
-
- 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] + image.data[i * width + j];
- }
- }
- return out;
- }
-}
-
-Img Img::operator-(const Img& image) const
-{
- if (this->width != image.getWidth() || this->height != image.getHeight()) {
- cerr << "The height or width of one image do not match with that of the other one!" << endl;
- return Img();
- }
- else {
- 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);
-
- 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] - image.data[i * width + j];
- }
- }
- return out;
- }
-}
-
-Img Img::operator*(const Img& image) const
-{
- if (width != image.getWidth() || height != image.getHeight()) {
- throw std::runtime_error("The height or width of one image do"
- " not match with that of the other one.");
- }
- auto out = Img(width, height, w_e_res, n_s_res);
-
- std::transform(data, data + (width * height), image.data, out.data,
- [](const int& a, const int& b) { return a * b; });
- return out;
-}
-
-Img Img::operator/(const Img& image) const
-{
- if (width != image.getWidth() || height != image.getHeight()) {
- throw std::runtime_error("The height or width of one image do"
- " not match with that of the other one.");
- }
- auto out = Img(width, height, w_e_res, n_s_res);
-
- std::transform(data, data + (width * height), image.data, out.data,
- [](const int& a, const int& b) { return a / b; });
- return out;
-}
-
-Img Img::operator*(double value) const
-{
- auto out = Img(width, height, w_e_res, n_s_res);
-
- std::transform(data, data + (width * height), out.data,
- [&value](const int& a) { return a * value; });
- return out;
-}
-
-Img Img::operator/(double value) const
-{
- auto out = Img(width, height, w_e_res, n_s_res);
-
- std::transform(data, data + (width * height), out.data,
- [&value](const int& a) { return a / value; });
- return out;
-}
-
-Img& Img::operator+=(int value)
-{
- std::for_each(data, data + (width * height),
- [&value](int& a) { a += value; });
- return *this;
-}
-
-Img& Img::operator-=(int value)
-{
- std::for_each(data, data + (width * height),
- [&value](int& a) { a -= value; });
- return *this;
-}
-
-Img& Img::operator*=(double value)
-{
- std::for_each(data, data + (width * height),
- [&value](int& a) { a *= value; });
- return *this;
-}
-
-Img& Img::operator/=(double value)
-{
- std::for_each(data, data + (width * height),
- [&value](int& a) { a /= value; });
- return *this;
-}
-
-Img& Img::operator+=(const Img& image)
-{
- for_each_zip(data, data + (width * height), image.data,
- [](int& a, int& b) { a += b; });
- return *this;
-}
-
-Img& Img::operator-=(const Img& image)
-{
- for_each_zip(data, data + (width * height), image.data,
- [](int& a, int& b) { a -= b; });
- return *this;
-}
-
-Img& Img::operator*=(const Img& image)
-{
- for_each_zip(data, data + (width * height), image.data,
- [](int& a, int& b) { a *= b; });
- return *this;
-}
-
-Img& Img::operator/=(const Img& image)
-{
- for_each_zip(data, data + (width * height), image.data,
- [](int& a, int& b) { a /= b; });
- return *this;
-}
-
-Img operator*(double factor, const Img& image)
-{
- return image * factor;
-}
-
-Img::~Img()
-{
- if (data) {
- delete[] data;
- }
-}
-
-void Img::toGrassRaster(const char *name)
-{
- int fd = Rast_open_new(name, CELL_TYPE);
- for (int i = 0; i < height; i++)
- Rast_put_c_row(fd, data + (i * width));
- Rast_close(fd);
-}
-
-// ref_name file is used to retrieve transformation and projection
-// information from the known (input) file
-void Img::toGdal(const char *name, const char *ref_name)
-{
-#ifdef SOD_GDAL_SUPPORT
- const char *format = "GTiff";
-
- GDALAllRegister();
-
- // obtain information for output Geotiff images
- GDALDataset *inputDataset = (GDALDataset *) GDALOpen(ref_name,
- GA_ReadOnly);
- double inputAdfGeoTransform[6];
- inputDataset->GetGeoTransform(inputAdfGeoTransform);
-
- // setup driver
- GDALDriver *gdalDriver = GetGDALDriverManager()->GetDriverByName(format);
-
- // set output Dataset and create output geotiff
- char **papszOptions = NULL;
- GDALDataset *outDataset = gdalDriver->Create(name, width, height, 1,
- GDT_Byte, papszOptions);
- outDataset->SetGeoTransform(inputAdfGeoTransform);
- outDataset->SetProjection(inputDataset->GetProjectionRef());
- GDALRasterBand *outBand = outDataset->GetRasterBand(1);
- CPLErr error = outBand->RasterIO(GF_Write, 0, 0, width, height,
- data, width, height,
- GDT_Int32, 0, 0);
- if (error == CE_Failure)
- throw std::runtime_error(string("Writing raster failed"
- " in GDAL RasterIO: ")
- + CPLGetLastErrorMsg());
- GDALClose((GDALDatasetH) outDataset);
- GDALClose((GDALDatasetH) inputDataset);
- CSLDestroy(papszOptions);
-#else
- throw std::runtime_error("GDAL support not available");
-#endif
-}
Deleted: grass-addons/grass7/raster/r.spread.sod/Img.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.h 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/Img.h 2018-06-21 20:09:01 UTC (rev 72867)
@@ -1,127 +0,0 @@
-/*
- * SOD model - raster manipulation
- *
- * Copyright (C) 2015-2017 by the authors.
- *
- * Authors: Zexi Chen (zchen22 ncsu edu)
- * Vaclav Petras (wenzeslaus gmail com)
- *
- * The code contained herein is licensed under the GNU General Public
- * License. You may obtain a copy of the GNU General Public License
- * Version 2 or later at the following locations:
- *
- * http://www.opensource.org/licenses/gpl-license.html
- * http://www.gnu.org/copyleft/gpl.html
- */
-
-
-#ifndef IMG_H
-#define IMG_H
-
-#include <iostream>
-#include <string>
-#include <cmath>
-#include <algorithm>
-#include <stdlib.h>
-
-
-enum Direction
-{
- N = 0, NE = 45, E = 90, SE = 135, S = 180, SW = 225, W = 270, NW = 315, NONE // NO means that there is no wind
-};
-
-class Img
-{
-private:
- int width;
- int height;
- // the west-east resolution of the pixel
- int w_e_res;
- // the north-south resolution of the pixel
- int n_s_res;
- int *data;
-public:
- 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);
- Img(int width, int height, int w_e_res, int n_s_res, int value);
- Img& operator=(Img&& other);
- Img& operator=(const Img& other);
-
- int getWidth() const
- {
- return width;
- }
-
- int getHeight() const
- {
- return height;
- }
-
- int getWEResolution() const
- {
- return w_e_res;
- }
-
- int getNSResolution() const
- {
- return n_s_res;
- }
-
- void fill(int value)
- {
- std::fill(data, data + (width * height), value);
- }
-
- void zero()
- {
- std::fill(data, data + (width * height), 0);
- }
-
- template<class UnaryOperation>
- void for_each(UnaryOperation op)
- {
- std::for_each(data, data + (width * height), op);
- }
-
- const int& operator()(unsigned row, unsigned col) const
- {
- return data[row * width + col];
- }
-
- int& operator()(unsigned row, unsigned col)
- {
- return data[row * width + col];
- }
-
- Img operator+(const Img& image) const;
- Img operator-(const Img& image) const;
- Img operator*(const Img& image) const;
- Img operator/(const Img& image) const;
- Img operator*(double value) const;
- Img operator/(double value) const;
- Img& operator+=(int value);
- Img& operator-=(int value);
- Img& operator*=(double value);
- Img& operator/=(double value);
- Img& operator+=(const Img& image);
- 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);
- void toGdal(const char *name, const char *ref_name);
-
- static Img fromGrassRaster(const char *name);
-};
-
-#endif
Modified: grass-addons/grass7/raster/r.spread.sod/Spore.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.cpp 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.cpp 2018-06-21 20:09:01 UTC (rev 72867)
@@ -95,6 +95,19 @@
generator.seed(random_seed);
}
+void Sporulation::SporeRemove(Img& I, Img& S, const DImg temperature,
+ double critical_temperature)
+{
+ for (int i = 0; i < height; i++) {
+ for (int j = 0; j < width; j++) {
+ if (temperature(i, j) < critical_temperature) {
+ S(i, j) += I(i, j); // move back to suseptable pool
+ I(i, j) = 0; // remove all infection
+ }
+ }
+ }
+}
+
void Sporulation::SporeGen(const Img& I, const double *weather,
double weather_value, double rate)
{
Modified: grass-addons/grass7/raster/r.spread.sod/Spore.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.h 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.h 2018-06-21 20:09:01 UTC (rev 72867)
@@ -19,7 +19,7 @@
#ifndef SPORE_H
#define SPORE_H
-#include "Img.h"
+#include "raster.h"
#include <random>
@@ -29,6 +29,12 @@
CAUCHY, CAUCHY_MIX // NO means that there is no wind
};
+// NONE means that there is no wind
+enum Direction
+{
+ N = 0, NE = 45, E = 90, SE = 135, S = 180, SW = 225, W = 270, NW = 315, NONE
+};
+
class Sporulation
{
private:
@@ -42,6 +48,8 @@
std::default_random_engine generator;
public:
Sporulation(unsigned random_seed, const Img &size);
+ void SporeRemove(Img& I, Img &S, const DImg temperature,
+ double critical_temperature);
void SporeGen(const Img& I, const double *weather,
double weather_value, double rate);
void SporeSpreadDisp_singleSpecies(Img& S, Img& I, Img& I2,
Modified: grass-addons/grass7/raster/r.spread.sod/date.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/date.h 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/date.h 2018-06-21 20:09:01 UTC (rev 72867)
@@ -35,10 +35,12 @@
Date(const Date &d): year(d.year), month(d.month), day(d.day){}
Date(int y, int m, int d): year(y), month(m), day(d){}
Date(): year(2000), month(1), day(1){}
- void increasedByWeek();
- Date getYearEnd();
- Date getNextYearEnd();
- bool isYearEnd();
+ inline void increasedByWeek();
+ inline void increasedByMonth();
+ inline Date getYearEnd();
+ inline Date getNextYearEnd();
+ inline bool isYearEnd();
+ inline bool isLastMonthOfYear();
int getMonth() const {return month;}
int getYear() const { return year;}
int getDay() const {return day;}
@@ -45,12 +47,12 @@
void setMonth(int m){month = m;}
void setYear(int y){year = y;}
void setDay(int d){day = d;}
- int weeksFromDate(Date start);
- friend std::ostream& operator<<(std::ostream& os, const Date &d);
- friend bool operator> (const Date &d1, const Date &d2);
- friend bool operator>= (const Date &d1, const Date &d2);
- friend bool operator< (const Date &d1, const Date &d2);
- friend bool operator<= (const Date &d1, const Date &d2);
+ inline int weeksFromDate(Date start);
+ inline friend std::ostream& operator<<(std::ostream& os, const Date &d);
+ inline friend bool operator> (const Date &d1, const Date &d2);
+ inline friend bool operator>= (const Date &d1, const Date &d2);
+ inline friend bool operator< (const Date &d1, const Date &d2);
+ inline friend bool operator<= (const Date &d1, const Date &d2);
};
std::ostream& operator<<(std::ostream& os, const Date &d)
@@ -69,6 +71,12 @@
return false;
}
+bool Date::isLastMonthOfYear(){
+ if (month == 12)
+ return true;
+ return false;
+}
+
Date Date::getNextYearEnd(){
if (month == 1)
return Date(year, 12, 31);
@@ -151,6 +159,25 @@
}
}
+void Date::increasedByMonth()
+{
+ month += 1;
+ if (month > 12) {
+ year++;
+ month = 1;
+ }
+ if (year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)) {
+ if (day > day_in_month[1][month]) {
+ day = day_in_month[1][month];
+ }
+ }
+ else {
+ if (day > day_in_month[0][month]) {
+ day = day_in_month[0][month];
+ }
+ }
+}
+
int Date::weeksFromDate(Date start) {
int week = 0;
Modified: grass-addons/grass7/raster/r.spread.sod/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/main.cpp 2018-06-21 19:49:23 UTC (rev 72866)
+++ grass-addons/grass7/raster/r.spread.sod/main.cpp 2018-06-21 20:09:01 UTC (rev 72867)
@@ -21,7 +21,7 @@
// #define SOD_NETCDF_SUPPORT
#include "date.h"
-#include "Img.h"
+#include "raster.h"
#include "Spore.h"
extern "C" {
@@ -135,15 +135,11 @@
" value '" + text +"' provided");
}
-bool seasonality_from_string(const string& text)
+typedef std::pair<int, int> Season;
+
+inline Season seasonality_from_option(const Option* opt)
{
- if (text == "yes")
- return true;
- else if (text == "no")
- return false;
- else
- throw std::invalid_argument("seasonality_from_string: Invalid"
- " value '" + text +"' provided");
+ return {std::atoi(opt->answers[0]), std::atoi(opt->answers[1])};
}
void read_names(std::vector<string>& names, const char* filename)
@@ -242,7 +238,10 @@
{
struct Option *species, *lvtree, *infected, *outside_spores;
struct Option *nc_weather, *moisture_file, *temperature_file, *weather_value, *weather_file;
+ struct Option *lethal_temperature_value, *lethal_temperature_months;
+ struct Option *actual_temperature_file;
struct Option *start_time, *end_time, *seasonality;
+ struct Option *step;
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;
@@ -343,13 +342,11 @@
opt.wind->required = YES;
opt.wind->guisection = _("Weather");
-#ifdef SOD_NETCDF_SUPPORT
opt.nc_weather = G_define_standard_option(G_OPT_F_BIN_INPUT);
opt.nc_weather->key = "ncdf_weather";
opt.nc_weather->description = _("Weather data");
opt.nc_weather->required = NO;
opt.nc_weather->guisection = _("Weather");
-#endif
opt.moisture_file = G_define_standard_option(G_OPT_F_INPUT);
opt.moisture_file->key = "moisture_file";
@@ -387,6 +384,40 @@
opt.weather_value->required = NO;
opt.weather_value->guisection = _("Weather");
+ opt.lethal_temperature_value = G_define_option();
+ opt.lethal_temperature_value->type = TYPE_DOUBLE;
+ opt.lethal_temperature_value->key = "lethal_temperature";
+ opt.lethal_temperature_value->label =
+ _("Temperature when the pest or patogen dies");
+ opt.lethal_temperature_value->description =
+ _("The temerature unit must be the same as for the"
+ "temerature raster map (typically degrees of Celsius)");
+ opt.lethal_temperature_value->required = NO;
+ opt.lethal_temperature_value->multiple = NO;
+ opt.lethal_temperature_value->guisection = _("Weather");
+
+ opt.lethal_temperature_months = G_define_option();
+ opt.lethal_temperature_months->type = TYPE_INTEGER;
+ opt.lethal_temperature_months->key = "lethal_month";
+ opt.lethal_temperature_months->label =
+ _("Month when the pest or patogen dies due to low temperature");
+ opt.lethal_temperature_months->description =
+ _("The temerature unit must be the same as for the"
+ "temerature raster map (typically degrees of Celsius)");
+ // TODO: implement this as multiple
+ opt.lethal_temperature_months->required = NO;
+ opt.lethal_temperature_months->guisection = _("Weather");
+
+ // TODO: rename coefs in interface and improve their descs
+ opt.actual_temperature_file = G_define_standard_option(G_OPT_F_INPUT);
+ opt.actual_temperature_file->key = "actual_temperature_file";
+ opt.actual_temperature_file->label =
+ _("Input file with one temperature raster map name per line");
+ opt.actual_temperature_file->description =
+ _("The temperature should be in actual temperature units (typically degrees of Celsius)");
+ opt.actual_temperature_file->required = NO;
+ opt.actual_temperature_file->guisection = _("Weather");
+
opt.start_time = G_define_option();
opt.start_time->type = TYPE_INTEGER;
opt.start_time->key = "start_time";
@@ -406,12 +437,27 @@
opt.seasonality = G_define_option();
opt.seasonality->type = TYPE_STRING;
opt.seasonality->key = "seasonality";
- opt.seasonality->label = _("Seasonal spread");
- opt.seasonality->description = _("Spread limited to certain months (season)");
- opt.seasonality->options = "yes,no";
- opt.seasonality->answer = "yes";
+ opt.seasonality->label = _("Seasonal spread (from,to)");
+ opt.seasonality->description =
+ _("Spread limited to certain months (season), for example"
+ " 5,9 for spread starting at the beginning of May and"
+ " ending at the end of September");
+ opt.seasonality->key_desc = "from,to";
+ //opt.seasonality->options = "1-12";
+ opt.seasonality->answer = "1,12";
+ opt.seasonality->multiple = NO;
opt.seasonality->guisection = _("Time");
+ opt.step = G_define_option();
+ opt.step->type = TYPE_STRING;
+ opt.step->key = "step";
+ opt.step->label = _("Step of simulation");
+ opt.step->description = _("How often the simulation computes new step");
+ opt.step->options = "week,month";
+ opt.step->descriptions = _("week;Compute next simulation step each week;month;Compute next simulation step each month");
+ opt.step->required = YES;
+ opt.step->guisection = _("Time");
+
opt.spore_rate = G_define_option();
opt.spore_rate->type = TYPE_DOUBLE;
opt.spore_rate->key = "spore_rate";
@@ -558,6 +604,13 @@
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+#ifndef SOD_NETCDF_SUPPORT
+ if (opt.nc_weather->answer) {
+ G_fatal_error(_("Direct NetCDF support is not available in this"
+ " installation, import the data instead"));
+ }
+#endif
+
unsigned num_runs = 1;
if (opt.runs->answer)
num_runs = std::stoul(opt.runs->answer);
@@ -572,7 +625,7 @@
file_exists_or_fatal_error(opt.weather_file);
// Seasonality: Do you want the spread to be limited to certain months?
- bool ss = seasonality_from_string(opt.seasonality->answer);
+ Season season = seasonality_from_option(opt.seasonality);
Direction pwdir = direction_enum_from_string(opt.wind->answer);
@@ -604,11 +657,14 @@
cerr << "Start date must precede the end date!!!" << endl;
exit(EXIT_FAILURE);
}
+
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;
+ string step = opt.step->answer;
+
// mortality
bool mortality = false;
unsigned first_year_to_die = 1; // starts at 1 (same as the opt)
@@ -716,6 +772,24 @@
}
#endif
+ double use_lethal_temperature = false;
+ double lethal_temperature_value;
+ int lethal_temperature_month = 0; // invalid value for month
+ std::vector<string> actual_temperature_names;
+ std::vector<DImg> actual_temperatures;
+ if (opt.lethal_temperature_value->answer)
+ lethal_temperature_value = std::stod(opt.lethal_temperature_value->answer);
+ if (opt.lethal_temperature_months->answer)
+ lethal_temperature_month = std::stod(opt.lethal_temperature_months->answer);
+ if (opt.actual_temperature_file->answer) {
+ file_exists_or_fatal_error(opt.actual_temperature_file);
+ read_names(actual_temperature_names, opt.actual_temperature_file->answer);
+ for (string name : actual_temperature_names) {
+ actual_temperatures.push_back(DImg::fromGrassRaster(name.c_str()));
+ }
+ use_lethal_temperature = true;
+ }
+
const unsigned max_weeks_in_year = 53;
double *mcf = nullptr;
double *ccf = nullptr;
@@ -754,11 +828,30 @@
Date dd_current(dd_start);
// main simulation loop (weekly steps)
- for (int current_week = 0; ; current_week++, dd_current.increasedByWeek()) {
+ for (int current_week = 0; ; current_week++, step == "month" ? dd_current.increasedByMonth() : dd_current.increasedByWeek()) {
if (dd_current < dd_end)
- if (!ss || !(dd_current.getMonth() > 9))
+ if (season.first >= dd_current.getMonth() && dd_current.getMonth() <= season.second)
unresolved_weeks.push_back(current_week);
+ // removal is out of sync with the actual runs but it does
+ // not matter as long as removal happends out of season
+ if (use_lethal_temperature
+ && dd_current.getMonth() == lethal_temperature_month
+ && (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();
+ if (simulation_year >= actual_temperatures.size())
+ G_fatal_error(_("Not enough temperatures"));
+ #pragma omp parallel for num_threads(threads)
+ for (unsigned run = 0; run < num_runs; run++) {
+ sporulations[run].SporeRemove(inf_species_rasts[run],
+ sus_species_rasts[run],
+ actual_temperatures[simulation_year],
+ lethal_temperature_value);
+ }
+ }
+
// if all the oaks are infected, then exit
if (all_infected(S_species_rast)) {
cerr << "In the " << dd_current << " all suspectible oaks are infected!" << endl;
@@ -766,7 +859,7 @@
}
// check whether the spore occurs in the month
- if (dd_current.isYearEnd() || dd_current >= dd_end) {
+ if ((step == "month" ? dd_current.isLastMonthOfYear() : dd_current.isYearEnd()) || dd_current >= dd_end) {
if (!unresolved_weeks.empty()) {
unsigned week_in_chunk = 0;
@@ -797,7 +890,10 @@
unsigned week_in_chunk = 0;
// actual runs of the simulation per week
for (auto week : unresolved_weeks) {
- double *week_weather = weather + week_in_chunk * width * height;
+ double *week_weather = 0;
+ if (weather) {
+ week_weather = weather + week_in_chunk * width * height;
+ }
if (!weather_coeff && !weather_values.empty()) {
weather_value = weather_values[week];
}
Added: grass-addons/grass7/raster/r.spread.sod/raster.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/raster.h (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/raster.h 2018-06-21 20:09:01 UTC (rev 72867)
@@ -0,0 +1,433 @@
+#ifndef RASTER_H
+#define RASTER_H
+
+/*
+ * SOD model - raster manipulation
+ *
+ * Copyright (C) 2015-2018 by the authors.
+ *
+ * Authors: Vaclav Petras <wenzeslaus gmail com>
+ * Completely rewritten by Vaclav Petras based on
+ * version by Zexi Chen <zchen22 ncsu edu>.
+ *
+ * The code contained herein is licensed under the GNU General Public
+ * License. You may obtain a copy of the GNU General Public License
+ * Version 2 or later at the following locations:
+ *
+ * http://www.opensource.org/licenses/gpl-license.html
+ * http://www.gnu.org/copyleft/gpl.html
+ */
+
+#include <iostream>
+#include <string>
+#include <cmath>
+#include <algorithm>
+#include <stdlib.h>
+
+extern "C" {
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+}
+
+#include <algorithm>
+#include <stdexcept>
+
+using std::string;
+using std::cerr;
+using std::endl;
+
+/* Iterate over two ranges and apply a binary function which modifies
+ * the first parameter.
+ */
+template<class InputIt1, class InputIt2, class BinaryOperation>
+BinaryOperation for_each_zip(InputIt1 first1, InputIt1 last1, InputIt2 first2, BinaryOperation f) {
+ for (; first1 != last1; ++first1, ++first2) {
+ f(*first1, *first2);
+ }
+ return f;
+}
+
+template<typename Number>
+class Raster
+{
+private:
+ unsigned width;
+ unsigned height;
+ // the west-east resolution of the pixel
+ double w_e_res;
+ // the north-south resolution of the pixel
+ double n_s_res;
+ Number *data;
+public:
+ Raster()
+ {
+ width = 0;
+ height = 0;
+ w_e_res = 0;
+ n_s_res = 0;
+ data = NULL;
+ }
+
+ Raster(const Raster& other)
+ {
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = new Number[width * height];
+ std::copy(other.data, other.data + (width * height), data);
+ }
+
+ Raster(const Raster& other, Number value)
+ {
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = new Number[width * height]{value};
+ }
+
+ Raster(Raster&& other)
+ {
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = other.data;
+ other.data = nullptr;
+ }
+
+ Raster(int width, int height, int w_e_res, int n_s_res)
+ {
+ this->width = width;
+ this->height = height;
+ this->w_e_res = w_e_res;
+ this->n_s_res = n_s_res;
+ this->data = new Number[width * height];
+ }
+
+ // TODO: res are doubles
+ // TODO: size is unsigned?
+ Raster(int width, int height, int w_e_res, int n_s_res, int value)
+ {
+ this->width = width;
+ this->height = height;
+ this->w_e_res = w_e_res;
+ this->n_s_res = n_s_res;
+ this->data = new Number[width * height]{value};
+ }
+
+ int getWidth() const
+ {
+ return width;
+ }
+
+ int getHeight() const
+ {
+ return height;
+ }
+
+ int getWEResolution() const
+ {
+ return w_e_res;
+ }
+
+ int getNSResolution() const
+ {
+ return n_s_res;
+ }
+
+ void fill(Number value)
+ {
+ std::fill(data, data + (width * height), value);
+ }
+
+ void zero()
+ {
+ std::fill(data, data + (width * height), 0);
+ }
+
+ template<class UnaryOperation>
+ void for_each(UnaryOperation op)
+ {
+ std::for_each(data, data + (width * height), op);
+ }
+
+ const Number& operator()(unsigned row, unsigned col) const
+ {
+ return data[row * width + col];
+ }
+
+ Number& operator()(unsigned row, unsigned col)
+ {
+ return data[row * width + col];
+ }
+
+ Raster& operator=(const Raster& other)
+ {
+ if (this != &other)
+ {
+ if (data)
+ delete[] data;
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = new Number[width * height];
+ std::copy(other.data, other.data + (width * height), data);
+ }
+ return *this;
+ }
+
+ Raster& operator=(Raster&& other)
+ {
+ if (this != &other)
+ {
+ if (data)
+ delete[] data;
+ width = other.width;
+ height = other.height;
+ w_e_res = other.w_e_res;
+ n_s_res = other.n_s_res;
+ data = other.data;
+ other.data = nullptr;
+ }
+ return *this;
+ }
+
+ Raster operator+(const Raster& image) const
+ {
+ if (this->width != image.getWidth() || this->height != image.getHeight()) {
+ cerr << "The height or width of one image do not match with that of the other one!" << endl;
+ return Raster();
+ }
+ else {
+ auto re_width = this->width;
+ auto re_height = this->height;
+ auto out = Raster(re_width, re_height, this->w_e_res, this->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] + image.data[i * width + j];
+ }
+ }
+ return out;
+ }
+ }
+
+ Raster operator-(const Raster& image) const
+ {
+ if (this->width != image.getWidth() || this->height != image.getHeight()) {
+ cerr << "The height or width of one image do not match with that of the other one!" << endl;
+ return Raster();
+ }
+ else {
+ auto re_width = this->width;
+ auto re_height = this->height;
+ auto out = Raster(re_width, re_height, this->w_e_res, this->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] - image.data[i * width + j];
+ }
+ }
+ return out;
+ }
+ }
+
+ Raster operator*(const Raster& image) const
+ {
+ if (width != image.getWidth() || height != image.getHeight()) {
+ throw std::runtime_error("The height or width of one image do"
+ " not match with that of the other one.");
+ }
+ auto out = Raster(width, height, w_e_res, n_s_res);
+
+ std::transform(data, data + (width * height), image.data, out.data,
+ [](const Number& a, const Number& b) { return a * b; });
+ return out;
+ }
+
+ Raster operator/(const Raster& image) const
+ {
+ if (width != image.getWidth() || height != image.getHeight()) {
+ throw std::runtime_error("The height or width of one image do"
+ " not match with that of the other one.");
+ }
+ auto out = Raster(width, height, w_e_res, n_s_res);
+
+ std::transform(data, data + (width * height), image.data, out.data,
+ [](const Number& a, const Number& b) { return a / b; });
+ return out;
+ }
+
+ Raster operator*(double value) const
+ {
+ auto out = Raster(width, height, w_e_res, n_s_res);
+
+ std::transform(data, data + (width * height), out.data,
+ [&value](const Number& a) { return a * value; });
+ return out;
+ }
+
+ Raster operator/(double value) const
+ {
+ auto out = Raster(width, height, w_e_res, n_s_res);
+
+ std::transform(data, data + (width * height), out.data,
+ [&value](const Number& a) { return a / value; });
+ return out;
+ }
+
+ Raster& operator+=(Number value)
+ {
+ std::for_each(data, data + (width * height),
+ [&value](Number& a) { a += value; });
+ return *this;
+ }
+
+ Raster& operator-=(Number value)
+ {
+ std::for_each(data, data + (width * height),
+ [&value](Number& a) { a -= value; });
+ return *this;
+ }
+
+ Raster& operator*=(double value)
+ {
+ std::for_each(data, data + (width * height),
+ [&value](Number& a) { a *= value; });
+ return *this;
+ }
+
+ Raster& operator/=(double value)
+ {
+ std::for_each(data, data + (width * height),
+ [&value](Number& a) { a /= value; });
+ return *this;
+ }
+
+ Raster& operator+=(const Raster& image)
+ {
+ for_each_zip(data, data + (width * height), image.data,
+ [](Number& a, Number& b) { a += b; });
+ return *this;
+ }
+
+ Raster& operator-=(const Raster& image)
+ {
+ for_each_zip(data, data + (width * height), image.data,
+ [](Number& a, Number& b) { a -= b; });
+ return *this;
+ }
+
+ Raster& operator*=(const Raster& image)
+ {
+ for_each_zip(data, data + (width * height), image.data,
+ [](Number& a, Number& b) { a *= b; });
+ return *this;
+ }
+
+ Raster& operator/=(const Raster& image)
+ {
+ for_each_zip(data, data + (width * height), image.data,
+ [](Number& a, Number& b) { a /= b; });
+ return *this;
+ }
+
+ friend inline Raster operator*(double factor, const Raster& image)
+ {
+ return image * factor;
+ }
+
+ friend inline Raster pow(Raster image, double value) {
+ image.for_each([value](Number& a){a = std::pow(a, value);});
+ return image;
+ }
+ friend inline Raster sqrt(Raster image) {
+ image.for_each([](Number& a){a = std::sqrt(a);});
+ return image;
+ }
+
+
+ ~Raster()
+ {
+ if (data) {
+ delete[] data;
+ }
+ }
+
+ static inline Raster fromGrassRaster(const char *name)
+ {
+ int fd = Rast_open_old(name, "");
+
+ Raster img;
+
+ img.width = Rast_window_cols();
+ img.height = Rast_window_rows();
+
+ Cell_head region;
+ Rast_get_window(®ion);
+ img.w_e_res = region.ew_res;
+ img.n_s_res = region.ns_res;
+
+ img.data = new Number[img.height * img.width];
+
+ for (int row = 0; row < img.height; row++) {
+ Rast_get_d_row(fd, img.data + (row * img.width), row);
+ }
+
+ Rast_close(fd);
+ return img;
+ }
+
+ void inline toGrassRaster(const char *name)
+ {
+ int fd = Rast_open_new(name, DCELL_TYPE);
+ for (int i = 0; i < height; i++)
+ Rast_put_d_row(fd, data + (i * width));
+ Rast_close(fd);
+ }
+
+};
+
+template <>
+inline Raster<int> Raster<int>::fromGrassRaster(const char *name)
+{
+ int fd = Rast_open_old(name, "");
+
+ Raster img;
+
+ img.width = Rast_window_cols();
+ img.height = Rast_window_rows();
+
+ Cell_head region;
+ Rast_get_window(®ion);
+ img.w_e_res = region.ew_res;
+ img.n_s_res = region.ns_res;
+
+ img.data = new int[img.height * img.width];
+
+ for (int row = 0; row < img.height; row++) {
+ Rast_get_c_row(fd, img.data + (row * img.width), row);
+ }
+
+ Rast_close(fd);
+ return img;
+}
+
+template <>
+inline void Raster<int>::toGrassRaster(const char *name)
+{
+ int fd = Rast_open_new(name, CELL_TYPE);
+ for (int i = 0; i < height; i++)
+ Rast_put_c_row(fd, data + (i * width));
+ Rast_close(fd);
+}
+
+// convenient definitions, names for backwards compatibility
+typedef Raster<int> Img;
+typedef Raster<double> DImg;
+
+#endif
More information about the grass-commit
mailing list