[GRASS-SVN] r71505 - in grass-addons/grass7/raster: . r.spread.sod
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Sep 24 19:17:50 PDT 2017
Author: wenzeslaus
Date: 2017-09-24 19:17:49 -0700 (Sun, 24 Sep 2017)
New Revision: 71505
Added:
grass-addons/grass7/raster/r.spread.sod/
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/Makefile
grass-addons/grass7/raster/r.spread.sod/README.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
grass-addons/grass7/raster/r.spread.sod/r.spread.sod.html
Modified:
grass-addons/grass7/raster/Makefile
Log:
r.spread.sod: Sudden Oak Death spread model (for single host species)
Modified: grass-addons/grass7/raster/Makefile
===================================================================
--- grass-addons/grass7/raster/Makefile 2017-09-25 02:14:53 UTC (rev 71504)
+++ grass-addons/grass7/raster/Makefile 2017-09-25 02:17:49 UTC (rev 71505)
@@ -101,6 +101,7 @@
r.skyview \
r.smooth.seg \
r.soillossbare \
+ r.spread.sod \
r.stream.basins \
r.stream.channel \
r.stream.distance \
Property changes on: grass-addons/grass7/raster/r.spread.sod
___________________________________________________________________
Added: svn:ignore
+ OBJ.*
*.tmp.html
Added: grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,93 @@
+# Change Log
+
+All notable changes to this project should be documented in this file.
+
+The format is based on [Keep a Changelog](http://keepachangelog.com/).
+
+## 2017-09-05 - September 2017 update
+
+### Added
+
+- Long-range dispersal kernel (Anna Petrasova)
+ - Events are recorded.
+ - The affected points are exported as a vector map.
+- Weather coefficients as GRASS GIS raster maps (Vaclav Petras)
+ - Input weather coefficients are obtained from a file with list of map
+ names (ordering and naming is resolved separately, e.g. same raster
+ can be used multiple times, i.e. temporal oversampling is possible).
+ - Weather rasters are now automatically resampled on the fly to the
+ raster grid based the computational region.
+
+### Changed
+
+- Spread of SOD based on a single species (Anna Petrasova)
+ - Spread for UMCA and oak replaced by single species, assumed tanoak.
+
+### Fixed
+
+- Interface now checks if only one way of providing weather coefficients
+ was used. (Vaclav Petras)
+
+## 2017-01-28 - January 2017 status
+
+### Added
+
+- GRASS GIS interface added (Vaclav Petras)
+ - Main spatial inputs handled through GRASS GIS libraries so that
+ different extents and resolutions can be used. (Not implemented for
+ the weather time series.)
+ - Text and numerical inputs hardcoded in defines replaced by
+ description using GRASS GIS library. Command line parameters are now
+ parsed, checked including dependencies, and it is possible to open
+ GUI.
+ - All hardcoded parameters replaced by command line parameters.
+ - Basic descriptions for all parameters.
+ - GDAL input and output is in the code.
+- Multiple stochastic runs and parallelization (Vaclav Petras)
+ - Optionally outputs also stdandard deviation.
+ - Aggregates series output from runs.
+ - Creates OpenMP threads for chunks of weeks.
+ - User can set number of threads from command line (or GUI).
+- Simplified weather inputs (Vaclav Petras)
+ - Weather can be supplied as as a text file (spatially constant)
+ - Weather can be supplied as one variable (non-spatial and non-temporal)
+
+### Changed
+
+- Efficiency improvements (Vaclav Petras)
+ - Enable inlining of size getters of Img class which makes all_infected
+ function much faster.
+ - Move constructor and assignment operator added for cases when RVO
+ is not applied.
+ - Internal storage changed to one array (usually faster allocation).
+ - Use the interal data directly to write using GDAL output.
+- Code cleanup (Vaclav Petras)
+ - Time measurement code removed (use external tool such as time instead).
+ - Commented-out code for seeding by time removed.
+ - Use same API style for Von Mises as for std lib distributions.
+ - Creating data for Img outside of the object is avoided.
+ - Indexing the Img done using operator ().
+ - Using operators for all operations which fit semantically.
+ - Remove unused variables from the code.
+ - The 'using namespace' statement replaced by explicit use 'using' for
+ string and other classes or objects.
+- Date class API extended to provide readable comparison operators
+ replacing usage of method with unclear name (Anna Petrasova)
+- Random seed for generators is required or must be explicitly generated.
+- The cpl_string dependency was removed. (Vaclav Petras)
+- The sporulation object is now seeded only once in the beginning.
+ This changes the stochastic output of one run. (Anna Petrasova)
+
+### Fixed
+
+- Initialize memory for the sporulation object for the cases when it is
+ zero cases to fix conditional jump which depends on uninitialised
+ value. (Vaclav Petras)
+- Memory allocation and deallocation is done by the right pair of
+ malloc-free or new-delete. (Vaclav Petras)
+- Compilation output and other non-repository files removed from the
+ repository. (Vaclav Petras)
+- Von Mises distribution concentration parameter is float not integer.
+ (Vaclav Petras)
+- Copy of GDAL code was removed from the repository, using system GDAL
+ includes now. (Vaclav Petras)
Property changes on: grass-addons/grass7/raster/r.spread.sod/CHANGELOG.md
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/Img.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.cpp (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/Img.cpp 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,389 @@
+/*
+ * 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
+ */
+
+
+#include "Img.h"
+
+extern "C" {
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+}
+
+#include <gdal/gdal.h>
+#include <gdal/gdal_priv.h>
+
+#include <algorithm>
+
+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(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)
+{
+ 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);
+ }
+}
+
+
+// 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 factor) 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);
+
+ 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;
+ }
+ }
+ 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::~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)
+{
+ 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);
+}
Property changes on: grass-addons/grass7/raster/r.spread.sod/Img.cpp
___________________________________________________________________
Added: svn:mime-type
+ text/x-c++src
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/Img.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Img.h (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/Img.h 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,122 @@
+/*
+ * 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);
+ //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 factor) 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);
+ ~Img();
+
+ void toGrassRaster(const char *name);
+ void toGdal(const char *name, const char *ref_name);
+
+ static Img fromGrassRaster(const char *name);
+};
+
+#endif
Property changes on: grass-addons/grass7/raster/r.spread.sod/Img.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/Makefile
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/Makefile 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,17 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.spread.sod
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(VECTORLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP) $(VECTORDEP)
+EXTRA_LIBS = $(GDALLIBS) $(OMPLIB)
+EXTRA_CFLAGS = $(GDALCFLAGS) $(OMPCFLAGS) $(VECT_CFLAGS) -std=c++11
+EXTRA_INC = $(VECT_INC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+default: cmd
+endif
Property changes on: grass-addons/grass7/raster/r.spread.sod/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/README.md
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/README.md (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/README.md 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,44 @@
+# SOD-modeling-cpp
+recoding the model to create a c++ version of the SOD-model base on https://github.com/f-tonini/SOD-modeling.
+This repository contains the c++ version scripts used to develop a stochastic landscape spread model of forest pathogen *P. ramorum*.
+
+The reference paper: Ross K. Meentemeyer, Nik J. Cunniffe, Alex R. Cook, Joao A. N. Filipe, Richard D. Hunter, David M. Rizzo, and Christopher A. Gilligan 2011. Epidemiological modeling of invasion in heterogeneous landscapes: spread of sudden oak death in California (1990–2030). *Ecosphere* 2:art17. [http://dx.doi.org/10.1890/ES10-00192.1] (http://www.esajournals.org/doi/abs/10.1890/ES10-00192.1)
+
+## The files
+The main.cpp contains the main program to run.
+
+## To run the model
+
+You can use Linux to run the model in the following way.
+
+Open an terminal and install dependencies:
+
+ sudo apt-get install libgdal-dev libnetcdf-dev
+
+Download the model code as ZIP or using Git:
+
+ git clone ...
+
+Change the current directory to the model directory:
+
+ cd ...
+
+Compile:
+
+ make
+
+Run:
+
+ ./a.out > result.txt
+
+## Authors
+
+* Francesco Tonini (original R version)
+* Zexi Chen (C++ version)
+* Vaclav Petras (parallelization, GRASS interface)
+* Anna Petrasova (single species simulation)
+
+## License
+
+This program is free software under the GNU General Public License
+(>=v2). Read the file LICENSE for details.
Property changes on: grass-addons/grass7/raster/r.spread.sod/README.md
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/Spore.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.cpp (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.cpp 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,319 @@
+/*
+ * SOD model - spore simulation
+ *
+ * Copyright (C) 2015-2017 by the authors.
+ *
+ * Authors: Zexi Chen (zchen22 ncsu edu)
+ * Vaclav Petras (wenzeslaus gmail com)
+ * Anna Petrasova (kratochanna 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
+ */
+
+
+#include "Spore.h"
+
+#include <cmath>
+#include <tuple>
+
+// PI is used in the code and M_PI is not guaranteed
+// fix it, but prefer the system definition
+#ifndef M_PI
+ #define M_PI 3.14159265358979323846
+#endif
+#ifndef PI
+ #define PI M_PI
+#endif
+
+using std::cerr;
+using std::endl;
+
+/*
+Von Mises Distribution(Circular data distribution)
+
+mu is the mean angle, expressed in radians between 0 and 2*pi,
+and kappa is the concentration parameter, which must be greater
+than or equal to zero. If kappa is equal to zero, this distribution
+reduces to a uniform random angle over the range 0 to 2*pi
+*/
+class von_mises_distribution
+{
+public:
+ von_mises_distribution(double mu, double kappa)
+ : mu(mu), kappa(kappa), distribution(0.0, 1.0)
+ {}
+ template<class Generator>
+ double operator ()(Generator& generator)
+ {
+ double a, b, c, f, r, theta, u1, u2, u3, z;
+
+ if (kappa <= 1.e-06)
+ return 2 * PI * distribution(generator);
+
+ a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa);
+ b = (a - sqrt(2.0 * a)) / (2.0 * kappa);
+ r = (1.0 + b * b) / (2.0 * b);
+
+ while (true) {
+ u1 = distribution(generator);
+ z = cos(PI * u1);
+ f = (1.0 + r * z) / (r + z);
+ c = kappa * (r - f);
+ u2 = distribution(generator);
+ if (u2 <= c * (2.0 - c) || u2 < c * exp(1.0 - c))
+ break;
+ }
+
+ u3 = distribution(generator);
+ if (u3 > 0.5) {
+ theta = fmod(mu + acos(f), 2 * PI);
+ }
+ else {
+ theta = fmod(mu - acos(f), 2 * PI);
+ }
+ return theta;
+ }
+private:
+ double mu;
+ double kappa;
+ std::uniform_real_distribution<double> distribution;
+};
+
+Sporulation::Sporulation(unsigned random_seed, const Img& size)
+ :
+ width(size.getWidth()),
+ height(size.getHeight()),
+ w_e_res(size.getWEResolution()),
+ n_s_res(size.getNSResolution()),
+ sp(width, height, w_e_res, n_s_res)
+{
+ generator.seed(random_seed);
+}
+
+void Sporulation::SporeGen(const Img& I, const double *weather,
+ double weather_value, double rate)
+{
+ double lambda = 0;
+ for (int i = 0; i < height; i++) {
+ for (int j = 0; j < width; j++) {
+ if (I(i, j) > 0) {
+ if (weather)
+ lambda = rate * weather[i * width + j];
+ else
+ lambda = rate * weather_value;
+ int sum = 0;
+ std::poisson_distribution<int> distribution(lambda);
+
+ for (int k = 0; k < I(i, j); k++) {
+ sum += distribution(generator);
+ }
+ sp(i, j) = sum;
+ }
+ else {
+ sp(i, j) = 0;
+ }
+ }
+ }
+}
+
+void Sporulation::SporeSpreadDisp_singleSpecies(Img& S, Img& I,
+ const Img& lvtree_rast,
+ std::vector<std::tuple<int, int> >& outside_spores,
+ Rtype rtype, const double *weather,
+ double weather_value, double scale1,
+ double kappa, Direction wdir, double scale2,
+ double gamma)
+{
+ std::cauchy_distribution < double >distribution_cauchy_one(0.0, scale1);
+ std::cauchy_distribution < double >distribution_cauchy_two(0.0, scale2);
+
+ std::bernoulli_distribution distribution_bern(gamma);
+ std::uniform_real_distribution < double >distribution_uniform(0.0, 1.0);
+
+ if (wdir == NONE)
+ kappa = 0;
+ von_mises_distribution vonmisesvariate(wdir * PI / 180, kappa);
+
+ double dist = 0;
+ double theta = 0;
+ bool scale2_used = false;
+
+ for (int i = 0; i < height; i++) {
+ for (int j = 0; j < width; j++) {
+ if (sp(i, j) > 0) {
+ for (int k = 0; k < sp(i, j); k++) {
+
+ // generate the distance from cauchy distribution or cauchy mixture distribution
+ if (rtype == CAUCHY) {
+ dist = abs(distribution_cauchy_one(generator));
+ }
+ else if (rtype == CAUCHY_MIX) {
+ // use bernoulli distribution to act as the sampling with prob(gamma,1-gamma)
+ if (distribution_bern(generator)) {
+ dist = abs(distribution_cauchy_one(generator));
+ scale2_used = false;
+ }
+ else {
+ dist = abs(distribution_cauchy_two(generator));
+ scale2_used = true;
+ }
+ }
+ else {
+ cerr <<
+ "The paramter Rtype muse be set as either CAUCHY OR CAUCHY_MIX"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ theta = vonmisesvariate(generator);
+
+ int row = i - round(dist * cos(theta) / n_s_res);
+ int col = j + round(dist * sin(theta) / w_e_res);
+
+ if (row < 0 || row >= height || col < 0 || col >= width) {
+ // export only spores coming from long-range dispersal kernel outside of modeled area
+ if (scale2_used)
+ outside_spores.emplace_back(std::make_tuple(row, col));
+ continue;
+ }
+ if (S(row, col) > 0) {
+ double prob_S =
+ (double)(S(row, col)) /
+ lvtree_rast(row, col);
+ double U = distribution_uniform(generator);
+
+ if (weather)
+ prob_S *= weather[row * width + col];
+ else
+ prob_S *= weather_value;
+ if (U < prob_S) {
+ I(row, col) += 1;
+ S(row, col) -= 1;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+
+void Sporulation::SporeSpreadDisp(Img& S_umca, Img& S_oaks, Img& I_umca,
+ Img& I_oaks, const Img& lvtree_rast,
+ Rtype rtype, const double *weather,
+ double weather_value, double scale1,
+ double kappa, Direction wdir, double scale2,
+ double gamma)
+{
+ std::cauchy_distribution < double >distribution_cauchy_one(0.0, scale1);
+ std::cauchy_distribution < double >distribution_cauchy_two(0.0, scale2);
+
+ std::bernoulli_distribution distribution_bern(gamma);
+ std::uniform_real_distribution < double >distribution_uniform(0.0, 1.0);
+
+ if (wdir == NONE)
+ kappa = 0;
+ von_mises_distribution vonmisesvariate(wdir * PI / 180, kappa);
+
+ double dist = 0;
+ double theta = 0;
+
+ for (int i = 0; i < height; i++) {
+ for (int j = 0; j < width; j++) {
+ if (sp(i, j) > 0) {
+ for (int k = 0; k < sp(i, j); k++) {
+
+ // generate the distance from cauchy distribution or cauchy mixture distribution
+ if (rtype == CAUCHY) {
+ dist = abs(distribution_cauchy_one(generator));
+ }
+ else if (rtype == CAUCHY_MIX) {
+ // use bernoulli distribution to act as the sampling with prob(gamma,1-gamma)
+ if (distribution_bern(generator))
+ dist = abs(distribution_cauchy_one(generator));
+ else
+ dist = abs(distribution_cauchy_two(generator));
+ }
+ else {
+ cerr <<
+ "The paramter Rtype muse be set as either CAUCHY OR CAUCHY_MIX"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ theta = vonmisesvariate(generator);
+
+ int row = i - round(dist * cos(theta) / n_s_res);
+ int col = j + round(dist * sin(theta) / w_e_res);
+
+ if (row < 0 || row >= height)
+ continue;
+ if (col < 0 || col >= width)
+ continue;
+
+ if (row == i && col == j) {
+ if (S_umca(row, col) > 0 ||
+ S_oaks(row, col) > 0) {
+ double prob =
+ (double)(S_umca(row, col) +
+ S_oaks(row, col)) /
+ lvtree_rast(row, col);
+
+ double U = distribution_uniform(generator);
+
+ if (weather)
+ prob = prob * weather[row * width + col];
+ else
+ prob = prob * weather_value;
+
+ // if U < prob, then one host will become infected
+ if (U < prob) {
+ double prob_S_umca =
+ (double)(S_umca(row, col)) /
+ (S_umca(row, col) +
+ S_oaks(row, col));
+ double prob_S_oaks =
+ (double)(S_oaks(row, col)) /
+ (S_umca(row, col) +
+ S_oaks(row, col));
+
+ std::bernoulli_distribution
+ distribution_bern_prob(prob_S_umca);
+ if (distribution_bern_prob(generator)) {
+ I_umca(row, col) += 1;
+ S_umca(row, col) -= 1;
+ }
+ else {
+ I_oaks(row, col) += 1;
+ S_oaks(row, col) -= 1;
+ }
+ }
+ }
+ }
+ else {
+ if (S_umca(row, col) > 0) {
+ double prob_S_umca =
+ (double)(S_umca(row, col)) /
+ lvtree_rast(row, col);
+ double U = distribution_uniform(generator);
+
+ if (weather)
+ prob_S_umca *= weather[row * width + col];
+ else
+ prob_S_umca *= weather_value;
+ if (U < prob_S_umca) {
+ I_umca(row, col) += 1;
+ S_umca(row, col) -= 1;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+}
Property changes on: grass-addons/grass7/raster/r.spread.sod/Spore.cpp
___________________________________________________________________
Added: svn:mime-type
+ text/x-c++src
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/Spore.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/Spore.h (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/Spore.h 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,61 @@
+/*
+ * SOD model - spore simulation
+ *
+ * Copyright (C) 2015-2017 by the authors.
+ *
+ * Authors: Zexi Chen (zchen22 ncsu edu)
+ * Vaclav Petras (wenzeslaus gmail com)
+ * Anna Petrasova (kratochanna 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 SPORE_H
+#define SPORE_H
+
+#include "Img.h"
+
+#include <random>
+
+
+enum Rtype
+{
+ CAUCHY, CAUCHY_MIX // NO means that there is no wind
+};
+
+class Sporulation
+{
+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;
+ Img sp;
+ std::default_random_engine generator;
+public:
+ 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,
+ 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,
+ Direction wdir = NONE, double scale2 = 0.0,
+ double gamma = 0.0);
+ void SporeSpreadDisp(Img& S_umca, Img& S_oaks, Img& I_umca,
+ Img& I_oaks, const Img& lvtree_rast, Rtype rtype,
+ const double *weather, double weather_value,
+ double scale1, double kappa = 2,
+ Direction wdir = NONE, double scale2 = 0.0,
+ double gamma = 0.0);
+};
+
+#endif
Property changes on: grass-addons/grass7/raster/r.spread.sod/Spore.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/date.h
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/date.h (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/date.h 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,164 @@
+/*
+ * SOD model - date manipulation
+ *
+ * Copyright (C) 2015-2017 by the authors.
+ *
+ * Authors: Zexi Chen (zchen22 ncsu edu)
+ * Anna Petrasova
+ *
+ * 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 DATE
+#define DATE
+
+#include <iostream>
+
+class Date{
+
+private:
+ int year;
+ int month;
+ int day;
+ int day_in_month[2][13] = {
+ {0,31,28,31,30,31,30,31,31,30,31,30,31},
+ {0,31,29,31,30,31,30,31,31,30,31,30,31}
+ };
+
+public:
+ 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();
+ int getMonth() const {return month;}
+ int getYear() const { return year;}
+ int getDay() const {return day;}
+ 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);
+};
+
+std::ostream& operator<<(std::ostream& os, const Date &d)
+{
+ os << d.year << '-' << d.month << '-' << d.day;
+ return os;
+}
+
+Date Date::getYearEnd() {
+ return Date(year, 12, 31);
+}
+
+bool Date::isYearEnd(){
+ if (month == 12 && (day + 7) > 31)
+ return true;
+ return false;
+}
+
+Date Date::getNextYearEnd(){
+ if (month == 1)
+ return Date(year, 12, 31);
+ else
+ return Date(year + 1, 12, 31);
+}
+
+bool operator> (const Date &d1, const Date &d2)
+{
+ if(d1.year < d2.year)
+ return false;
+ else if (d1.year > d2.year)
+ return true;
+ else {
+ if (d1.month < d2.month)
+ return false;
+ else if (d1.month > d2.month)
+ return true;
+ else {
+ if (d1.day <= d2.day)
+ return false;
+ else
+ return true;
+ }
+ }
+}
+
+bool operator<= (const Date &d1, const Date &d2)
+{
+ return !(d1 > d2);
+}
+
+bool operator< (const Date &d1, const Date &d2)
+{
+ if(d1.year > d2.year)
+ return false;
+ else if (d1.year < d2.year)
+ return true;
+ else {
+ if (d1.month > d2.month)
+ return false;
+ else if (d1.month < d2.month)
+ return true;
+ else {
+ if (d1.day >= d2.day)
+ return false;
+ else
+ return true;
+ }
+ }
+}
+
+bool operator>= (const Date &d1, const Date &d2)
+{
+ return !(d1 < d2);
+}
+
+void Date::increasedByWeek()
+{
+ day += 7;
+ if (year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)) {
+ if (day > day_in_month[1][month]) {
+ day = day - day_in_month[1][month];
+ month++;
+ if (month > 12) {
+ year++;
+ month = 1;
+ }
+ }
+ }
+ else {
+ if (day > day_in_month[0][month]) {
+ day = day - day_in_month[0][month];
+ month++;
+ if (month > 12) {
+ year++;
+ month = 1;
+ }
+ }
+ }
+}
+
+int Date::weeksFromDate(Date start) {
+
+ int week = 0;
+ while (start <= *this) {
+ week++;
+ start.increasedByWeek();
+ }
+ return week - 1;
+}
+
+#endif // DATE
Property changes on: grass-addons/grass7/raster/r.spread.sod/date.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/main.cpp (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/main.cpp 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,777 @@
+/*
+ * SOD model
+ *
+ * Copyright (C) 2015-2017 by the authors.
+ *
+ * Authors: Zexi Chen (zchen22 ncsu edu)
+ * Vaclav Petras (wenzeslaus gmail com)
+ * Anna Petrasova (kratochanna 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_NETCDF_SUPPORT
+
+#include "date.h"
+#include "Img.h"
+#include "Spore.h"
+
+extern "C" {
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/vector.h>
+#include <grass/raster.h>
+}
+
+#ifdef SOD_NETCDF_SUPPORT
+#include <netcdfcpp.h>
+#endif
+
+#include <map>
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <fstream>
+#include <sstream>
+#include <string>
+
+using std::string;
+using std::cout;
+using std::cerr;
+using std::endl;
+
+
+#define DIM 1
+
+// Initialize infected trees for each species
+// needed unless empirical info is available
+static Img initialize(Img& img1,Img& img2) {
+ if (img1.getWidth() != img2.getWidth() ||
+ img2.getHeight() != img2.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 = img1.getWidth();
+ auto re_height = img1.getHeight();
+ auto out = Img(re_width, re_height, img1.getWEResolution(), img1.getNSResolution());
+
+ for (int i = 0; i < re_height; i++) {
+ for (int j = 0; j < re_width; j++) {
+ if (img2(i, j) > 0) {
+ if (img1(i, j) > img2(i, j))
+ out(i, j) =
+ img1(i, j) <
+ (img2(i, j) *
+ 2) ? img1(i, j) : (img2(i, j) * 2);
+ else
+ out(i, j) = img1(i, j);
+ }
+ else {
+ out(i, j) = 0;
+ }
+ }
+ }
+ return out;
+ }
+}
+
+string generate_name(const string& basename, const Date& date)
+{
+ // counting on year being 4 digits
+ auto year = G_double_to_basename_format(date.getYear(), 4, 0);
+ auto month = G_double_to_basename_format(date.getMonth(), 2, 0);
+ auto day = G_double_to_basename_format(date.getDay(), 2, 0);
+ auto sep = G_get_basename_separator();
+ string name = basename + sep + year + "_" + month + "_" + day;
+ return name;
+}
+
+Direction direction_enum_from_string(const string& text)
+{
+ std::map<string, Direction> mapping{
+ {"N", N}, {"NE", NE}, {"E", E}, {"SE", SE}, {"S", S},
+ {"SW", SW}, {"W", W}, {"NW", NW}, {"NONE", NONE}
+ };
+ try {
+ return mapping.at(text);
+ }
+ catch (const std::out_of_range&) {
+ throw std::invalid_argument("direction_enum_from_string: Invalid"
+ " value '" + text +"' provided");
+ }
+}
+
+Rtype radial_type_from_string(const string& text)
+{
+ if (text == "cauchy")
+ return CAUCHY;
+ else if (text == "cauchy_mix")
+ return CAUCHY_MIX;
+ else
+ throw std::invalid_argument("radial_type_from_string: Invalid"
+ " value '" + text +"' provided");
+}
+
+bool seasonality_from_string(const string& text)
+{
+ if (text == "yes")
+ return true;
+ else if (text == "no")
+ return false;
+ else
+ throw std::invalid_argument("seasonality_from_string: Invalid"
+ " value '" + text +"' provided");
+}
+
+void read_names(std::vector<string>& names, const char* filename)
+{
+ std::ifstream file(filename);
+ string line;
+ while (std::getline(file, line)) {
+ names.push_back(line);
+ }
+}
+
+std::vector<double> weather_file_to_list(const string& filename)
+{
+ std::ifstream input(filename);
+ std::vector<double> output;
+ string line;
+ while (std::getline(input, line))
+ {
+ double m, c;
+ std::istringstream stream(line);
+ stream >> m >> c;
+ output.push_back(m * c);
+ }
+ return output;
+}
+
+bool all_infected(Img& S_rast)
+{
+ bool allInfected = true;
+ for (int j = 0; j < S_rast.getHeight(); j++) {
+ for (int k = 0; k < S_rast.getWidth(); k++) {
+ if (S_rast(j, k) > 0)
+ allInfected = false;
+ }
+ }
+ return allInfected;
+}
+
+#ifdef SOD_NETCDF_SUPPORT
+void get_spatial_weather(NcVar *mcf_nc, NcVar *ccf_nc, double* mcf, double* ccf, double* weather, int width, int height, int step)
+{
+ // read the weather information
+ if (!mcf_nc->set_cur(step, 0, 0)) {
+ cerr << "Can not read the coefficients from the mcf_nc pointer to mcf array " << step << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (!ccf_nc->set_cur(step, 0, 0)) {
+ cerr << "Can not read the coefficients from the ccf_nc pointer to ccf array "<< step << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (!mcf_nc->get(mcf, 1, height, width)) {
+ cerr << "Can not get the record from mcf_nc " << step << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (!ccf_nc->get(ccf, 1, height, width)) {
+ cerr << "Can not get the record from ccf_nc " << step << endl;
+ exit(EXIT_FAILURE);
+ }
+ for (int j = 0; j < height; j++) {
+ for (int k = 0; k < width; k++) {
+ weather[j * width + k] = mcf[j * width + k] * ccf[j * width + k];
+ }
+ }
+}
+#endif
+
+// TODO: create image which is float/double/template
+void array_from_grass_raster(const char* name, double* array,
+ unsigned width, unsigned height)
+{
+ int fd = Rast_open_old(name, "");
+
+ for (int row = 0; row < height; row++) {
+ Rast_get_d_row(fd, array + (row * width), row);
+ }
+
+ Rast_close(fd);
+}
+
+void weather_rasters_to_array(const string& moisture_name,
+ const string& temperature_name,
+ double* mcf, double* ccf, double* weather,
+ unsigned width, unsigned height)
+{
+ // read the weather information
+ array_from_grass_raster(moisture_name.c_str(), mcf, width, height);
+ array_from_grass_raster(temperature_name.c_str(), ccf, width, height);
+ for (unsigned j = 0; j < height; j++) {
+ for (unsigned k = 0; k < width; k++) {
+ weather[j * width + k] = mcf[j * width + k] * ccf[j * width + k];
+ }
+ }
+}
+
+struct SodOptions
+{
+ struct Option *species, *lvtree, *infected, *outside_spores;
+ struct Option *nc_weather, *moisture_file, *temperature_file, *weather_value, *weather_file;
+ struct Option *start_time, *end_time, *seasonality;
+ struct Option *spore_rate, *wind;
+ struct Option *radial_type, *scale_1, *scale_2, *kappa, *gamma;
+ struct Option *seed, *runs, *threads;
+ struct Option *output, *output_series;
+ struct Option *stddev, *stddev_series;
+};
+
+struct SodFlags
+{
+ struct Flag *generate_seed;
+};
+
+
+int main(int argc, char *argv[])
+{
+ SodOptions opt;
+ SodFlags flg;
+
+ G_gisinit(argv[0]);
+
+ struct GModule *module = G_define_module();
+
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("spread"));
+ G_add_keyword(_("model"));
+ G_add_keyword(_("disease"));
+ module->description = _("Stochastic landscape spread model of forest pathogen - Sudden Oak Death (SOD)");
+
+ opt.species = G_define_standard_option(G_OPT_R_INPUT);
+ opt.species->key = "species";
+ opt.species->description = _("Input infected species raster map");
+ opt.species->guisection = _("Input");
+
+ opt.lvtree = G_define_standard_option(G_OPT_R_INPUT);
+ opt.lvtree->key = "lvtree";
+ opt.lvtree->description = _("Input live tree (all) raster map");
+ opt.lvtree->guisection = _("Input");
+
+ // TODO: is this oaks?
+ opt.infected = G_define_standard_option(G_OPT_R_INPUT);
+ opt.infected->key = "infected";
+ opt.infected->description = _("Initial sources of infection raster map");
+ opt.infected->guisection = _("Input");
+
+ opt.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.output->guisection = _("Output");
+
+ opt.output_series = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
+ opt.output_series->key = "output_series";
+ opt.output_series->description = _("Basename for output series");
+ opt.output_series->required = NO;
+ opt.output_series->guisection = _("Output");
+
+ opt.stddev = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.stddev->key = "stddev";
+ opt.stddev->description = _("Standard deviations");
+ opt.stddev->required = NO;
+ opt.stddev->guisection = _("Output");
+
+ opt.stddev_series = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
+ opt.stddev_series->key = "stddev_series";
+ opt.stddev_series->description
+ = _("Basename for output series of standard deviations");
+ opt.stddev_series->required = NO;
+ opt.stddev_series->guisection = _("Output");
+
+ opt.outside_spores = G_define_standard_option(G_OPT_V_OUTPUT);
+ opt.outside_spores->key = "outside_spores";
+ opt.outside_spores->required = NO;
+ opt.outside_spores->guisection = _("Output");
+
+ opt.wind = G_define_option();
+ opt.wind->type = TYPE_STRING;
+ opt.wind->key = "wind";
+ opt.wind->label = _("Prevailing wind direction");
+ opt.wind->description = _("NONE means that there is no wind");
+ opt.wind->options = "N,NE,E,SE,S,SW,W,NW,NONE";
+ 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";
+ opt.moisture_file->label =
+ _("Input file with one moisture map name per line");
+ opt.moisture_file->description =
+ _("Moisture coefficient");
+ opt.moisture_file->required = NO;
+ opt.moisture_file->guisection = _("Weather");
+
+ opt.temperature_file = G_define_standard_option(G_OPT_F_INPUT);
+ opt.temperature_file->key = "temperature_file";
+ opt.temperature_file->label =
+ _("Input file with one temperature map name per line");
+ opt.temperature_file->description =
+ _("Temperature coefficient");
+ opt.temperature_file->required = NO;
+ opt.temperature_file->guisection = _("Weather");
+
+ opt.weather_file = G_define_standard_option(G_OPT_F_INPUT);
+ opt.weather_file->key = "weather_file";
+ opt.weather_file->label = _("Text file with weather");
+ opt.weather_file->description =
+ _("Moisture and temperature");
+ opt.weather_file->required = NO;
+ opt.weather_file->guisection = _("Weather");
+
+ opt.weather_value = G_define_option();
+ opt.weather_value->type = TYPE_INTEGER;
+ opt.weather_value->key = "weather_value";
+ opt.weather_value->label = _("Value to be used as weather coeficient");
+ opt.weather_value->description =
+ _("Spatially and temporally constant weather coeficient"
+ " (usually moisture times temperture)");
+ opt.weather_value->required = NO;
+ opt.weather_value->guisection = _("Weather");
+
+ opt.start_time = G_define_option();
+ opt.start_time->type = TYPE_INTEGER;
+ opt.start_time->key = "start_time";
+ opt.start_time->label = _("Start year for the simulation");
+ opt.start_time->description = _("The first day of the year will be used");
+ opt.start_time->required = YES;
+ opt.start_time->guisection = _("Time");
+
+ opt.end_time = G_define_option();
+ opt.end_time->type = TYPE_INTEGER;
+ opt.end_time->key = "end_time";
+ opt.end_time->label = _("End year for the simulation");
+ opt.end_time->description = _("The last day of the year will be used");
+ opt.end_time->required = YES;
+ opt.end_time->guisection = _("Time");
+
+ 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->guisection = _("Time");
+
+ opt.spore_rate = G_define_option();
+ opt.spore_rate->type = TYPE_DOUBLE;
+ opt.spore_rate->key = "spore_rate";
+ opt.spore_rate->label = _("Spore production rate per week for each infected tree");
+ opt.spore_rate->answer = "4.4";
+ opt.spore_rate->guisection = _("Spores");
+
+ opt.radial_type = G_define_option();
+ opt.radial_type->type = TYPE_STRING;
+ opt.radial_type->key = "radial_type";
+ opt.radial_type->label = _("Radial distribution type");
+ opt.radial_type->answer = "cauchy";
+ opt.radial_type->options = "cauchy,cauchy_mix";
+ opt.radial_type->guisection = _("Spores");
+
+ opt.scale_1 = G_define_option();
+ opt.scale_1->type = TYPE_DOUBLE;
+ opt.scale_1->key = "scale_1";
+ opt.scale_1->label = _("Scale parameter for the first Cauchy distribution");
+ opt.scale_1->answer = "20.57";
+ opt.scale_1->guisection = _("Spores");
+
+ opt.scale_2 = G_define_option();
+ opt.scale_2->type = TYPE_DOUBLE;
+ opt.scale_2->key = "scale_2";
+ opt.scale_2->label = _("Scale parameter for the second Cauchy distribution");
+ opt.scale_2->guisection = _("Spores");
+
+ opt.kappa = G_define_option();
+ opt.kappa->type = TYPE_DOUBLE;
+ opt.kappa->key = "kappa";
+ opt.kappa->label = _("Concentration parameter for the von Mises distribution");
+ opt.kappa->answer = "2";
+ opt.kappa->guisection = _("Spores");
+
+ opt.gamma = G_define_option();
+ opt.gamma->type = TYPE_DOUBLE;
+ opt.gamma->key = "gamma";
+ opt.gamma->label = _("Gamma parameter for Bernoulli distribution");
+ opt.gamma->description = _("Probability of using the first Cauchy distribution");
+ opt.gamma->options = "0-1";
+ opt.gamma->guisection = _("Spores");
+
+ opt.seed = G_define_option();
+ opt.seed->key = "random_seed";
+ opt.seed->type = TYPE_INTEGER;
+ opt.seed->required = NO;
+ opt.seed->label = _("Seed for random number generator");
+ opt.seed->description =
+ _("The same seed can be used to obtain same results"
+ " or random seed can be generated by other means.");
+ opt.seed->guisection = _("Randomness");
+
+ flg.generate_seed = G_define_flag();
+ flg.generate_seed->key = 's';
+ flg.generate_seed->label =
+ _("Generate random seed (result is non-deterministic)");
+ flg.generate_seed->description =
+ _("Automatically generates random seed for random number"
+ " generator (use when you don't want to provide the seed option)");
+ flg.generate_seed->guisection = _("Randomness");
+
+ opt.runs = G_define_option();
+ opt.runs->key = "runs";
+ opt.runs->type = TYPE_INTEGER;
+ opt.runs->required = NO;
+ opt.runs->label = _("Number of simulation runs");
+ opt.runs->description =
+ _("The individual runs will obtain different seeds"
+ " and will be avaraged for the output");
+ opt.runs->guisection = _("Randomness");
+
+ opt.threads = G_define_option();
+ opt.threads->key = "nprocs";
+ opt.threads->type = TYPE_INTEGER;
+ opt.threads->required = NO;
+ opt.threads->description =
+ _("Number of threads for parallel computing");
+ opt.threads->options = "1-";
+ opt.threads->guisection = _("Randomness");
+
+ G_option_exclusive(opt.seed, flg.generate_seed, NULL);
+ G_option_required(opt.seed, flg.generate_seed, NULL);
+
+ // weather
+ G_option_exclusive(opt.moisture_file, opt.nc_weather,
+ opt.weather_file, opt.weather_value, NULL);
+ G_option_exclusive(opt.temperature_file, opt.nc_weather,
+ opt.weather_file, opt.weather_value, NULL);
+ G_option_collective(opt.moisture_file, opt.temperature_file, NULL);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ unsigned num_runs = 1;
+ if (opt.runs->answer)
+ num_runs = std::stoul(opt.runs->answer);
+
+ unsigned threads = 1;
+ if (opt.threads->answer)
+ threads = std::stoul(opt.threads->answer);
+
+ // Seasonality: Do you want the spread to be limited to certain months?
+ bool ss = seasonality_from_string(opt.seasonality->answer);
+
+ Direction pwdir = direction_enum_from_string(opt.wind->answer);
+
+ // set the spore rate
+ double spore_rate = std::stod(opt.spore_rate->answer);
+ Rtype rtype = radial_type_from_string(opt.radial_type->answer);
+ double scale1 = std::stod(opt.scale_1->answer);
+ double scale2 = 0;
+ if (rtype == CAUCHY_MIX && !opt.scale_2->answer)
+ G_fatal_error(_("The option %s is required for %s=%s"),
+ opt.scale_2->key, opt.radial_type->key,
+ opt.radial_type->answer);
+ else if (opt.scale_2->answer)
+ scale2 = std::stod(opt.scale_2->answer);
+ double kappa = std::stod(opt.kappa->answer);
+ double gamma = 0.0;
+ if (rtype == CAUCHY_MIX && !opt.gamma->answer)
+ G_fatal_error(_("The option %s is required for %s=%s"),
+ opt.gamma->key, opt.radial_type->key,
+ opt.radial_type->answer);
+ else if (opt.gamma->answer)
+ gamma = std::stod(opt.gamma->answer);
+
+ // initialize the start Date and end Date object
+ // options for times are required ints
+ int start_time = std::stoi(opt.start_time->answer);
+ int end_time = std::stoi(opt.end_time->answer);
+ if (start_time > end_time) {
+ 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);
+
+ unsigned seed_value;
+ if (opt.seed->answer) {
+ seed_value = std::stoul(opt.seed->answer);
+ G_verbose_message(_("Read random seed from %s option: %ud"),
+ opt.seed->key, seed_value);
+ } else {
+ // flag of option is required
+ std::random_device rd;
+ seed_value = rd();
+ G_verbose_message(_("Generated random seed (-%c): %ud"),
+ flg.generate_seed->key, seed_value);
+ }
+
+ // read the suspectible UMCA raster image
+ Img species_rast = Img::fromGrassRaster(opt.species->answer);
+
+ // read the living trees raster image
+ Img lvtree_rast = Img::fromGrassRaster(opt.lvtree->answer);
+
+ // read the initial infected oaks image
+ Img I_species_rast = Img::fromGrassRaster(opt.infected->answer);
+
+ // create the initial suspectible oaks image
+ Img S_species_rast = species_rast - I_species_rast;
+
+ // SOD-immune trees image
+ //Img SOD_rast = umca_rast + oaks_rast;
+ //Img IMM_rast = lvtree_rast - SOD_rast;
+
+ // retrieve the width and height of the images
+ int width = species_rast.getWidth();
+ int height = species_rast.getHeight();
+
+#ifdef SOD_NETCDF_SUPPORT
+ std::shared_ptr<NcFile> weather_coeff = nullptr;
+#else
+ bool weather_coeff = false;
+#endif
+ std::vector<string> moisture_names;
+ std::vector<string> temperature_names;
+ double weather_from_rasters = false;
+ std::vector<double> weather_values;
+ double weather_value = 0;
+
+#ifdef SOD_NETCDF_SUPPORT
+ if (opt.nc_weather->answer)
+ weather_coeff = std::make_shared<NcFile>(opt.nc_weather->answer, NcFile::ReadOnly);
+#else
+ if (false)
+ ;
+ // TODO: switch order of ifs to avoid this
+#endif
+ else if (opt.moisture_file->answer && opt.temperature_file->answer) {
+ read_names(moisture_names, opt.moisture_file->answer);
+ read_names(temperature_names, opt.temperature_file->answer);
+ weather_from_rasters = true;
+ }
+ else if (opt.weather_file->answer)
+ weather_values = weather_file_to_list(opt.weather_file->answer);
+ else if (opt.weather_value->answer)
+ weather_value = std::stoi(opt.weather_value->answer);
+ else
+ weather_value = 1; // no change (used in multiplication)
+
+#ifdef SOD_NETCDF_SUPPORT
+ if (weather_coeff && !weather_coeff->is_valid()) {
+ cerr << "Can not open the weather coefficients file(.cn)!" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ // TODO: do we need to free this? docs is 404
+ NcVar *mcf_nc = nullptr;
+ NcVar *ccf_nc = nullptr;
+
+ if (weather_coeff && !(mcf_nc = weather_coeff->get_var("Mcoef"))) {
+ cerr << "Can not read the moisture coefficients from the file!" <<
+ endl;
+ exit(EXIT_FAILURE);
+ }
+
+ if (weather_coeff && !(ccf_nc = weather_coeff->get_var("Ccoef"))) {
+ cerr << "Can not read the temperature coefficients from the file!" <<
+ endl;
+ exit(EXIT_FAILURE);
+ }
+#endif
+
+ const unsigned max_weeks_in_year = 53;
+ double *mcf = nullptr;
+ double *ccf = nullptr;
+ double *weather = nullptr;
+ if (weather_coeff || weather_from_rasters) {
+ mcf = new double[height * width];
+ ccf = new double[height * width];
+ weather = new double[max_weeks_in_year * height * width];
+ }
+
+ // build the Sporulation object
+ 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);
+ sporulations.reserve(num_runs);
+ for (unsigned i = 0; i < num_runs; ++i)
+ sporulations.emplace_back(seed_value++, I_species_rast);
+ std::vector<std::vector<std::tuple<int, int> > > outside_spores(num_runs);
+
+ std::vector<unsigned> unresolved_weeks;
+ unresolved_weeks.reserve(max_weeks_in_year);
+
+ // 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))
+ 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;
+ break;
+ }
+
+ // check whether the spore occurs in the month
+ if (dd_start.isYearEnd() || dd_start >= dd_end) {
+ if (!unresolved_weeks.empty()) {
+
+ unsigned week_in_chunk = 0;
+ // get weather for all the weeks
+ for (auto week : unresolved_weeks) {
+ double *week_weather = weather + week_in_chunk * width * height;
+#ifdef SOD_NETCDF_SUPPORT
+ if (weather_coeff) {
+ get_spatial_weather(mcf_nc, ccf_nc, mcf, ccf, week_weather, width, height, week);
+ }
+#else
+ if (false)
+ ;
+ // TODO: switch order of ifs to avoid this
+#endif
+ else if (weather_from_rasters) {
+ weather_rasters_to_array(moisture_names[week],
+ temperature_names[week],
+ mcf, ccf, week_weather,
+ width, height);
+ }
+ ++week_in_chunk;
+ }
+
+ // stochastic simulation runs
+ #pragma omp parallel for num_threads(threads)
+ for (unsigned run = 0; run < num_runs; run++) {
+ 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;
+ if (!weather_coeff && !weather_values.empty()) {
+ weather_value = weather_values[week];
+ }
+ sporulations[run].SporeGen(inf_species_rasts[run], week_weather, weather_value, spore_rate);
+
+ sporulations[run].SporeSpreadDisp_singleSpecies(sus_species_rasts[run], inf_species_rasts[run],
+ lvtree_rast, outside_spores[run],
+ rtype, week_weather,
+ weather_value, scale1, kappa, pwdir, scale2,
+ gamma);
+ ++week_in_chunk;
+ }
+ }
+ unresolved_weeks.clear();
+ }
+ if (opt.output_series->answer || opt.stddev_series->answer) {
+ // aggregate
+ I_species_rast.zero();
+ for (unsigned i = 0; i < num_runs; i++)
+ I_species_rast += inf_species_rasts[i];
+ I_species_rast /= num_runs;
+ // write result
+ // date is always end of the year, even for seasonal spread
+ if (opt.output_series->answer) {
+ string name = generate_name(opt.output_series->answer, dd_start);
+ I_species_rast.toGrassRaster(name.c_str());
+ }
+ }
+ if (opt.stddev_series->answer) {
+ Img stddev(I_species_rast.getWidth(), I_species_rast.getHeight(),
+ I_species_rast.getWEResolution(), I_species_rast.getNSResolution(), 0);
+ for (unsigned i = 0; i < num_runs; i++) {
+ Img tmp = inf_species_rasts[i] - I_species_rast;
+ stddev += tmp * tmp;
+ }
+ stddev /= num_runs;
+ stddev.for_each([](int& a){a = std::sqrt(a);});
+ string name = generate_name(opt.stddev_series->answer, dd_start);
+ stddev.toGrassRaster(name.c_str());
+ }
+ }
+
+ if (dd_start >= dd_end)
+ break;
+ }
+
+ // aggregate
+ I_species_rast.zero();
+ for (unsigned i = 0; i < num_runs; i++)
+ I_species_rast += inf_species_rasts[i];
+ I_species_rast /= num_runs;
+ // write final result
+ I_species_rast.toGrassRaster(opt.output->answer);
+
+ if (opt.stddev->answer) {
+ Img stddev(I_species_rast.getWidth(), I_species_rast.getHeight(),
+ I_species_rast.getWEResolution(), I_species_rast.getNSResolution(), 0);
+ for (unsigned i = 0; i < num_runs; i++) {
+ Img tmp = inf_species_rasts[i] - I_species_rast;
+ stddev += tmp * tmp;
+ }
+ stddev /= num_runs;
+ stddev.for_each([](int& a){a = std::sqrt(a);});
+ stddev.toGrassRaster(opt.stddev->answer);
+ }
+ if (opt.outside_spores->answer) {
+ Cell_head region;
+ Rast_get_window(®ion);
+ struct Map_info Map;
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ if (Vect_open_new(&Map, opt.outside_spores->answer, WITHOUT_Z) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"), opt.outside_spores->answer);
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ for (unsigned i = 0; i < num_runs; i++) {
+ for (unsigned j = 0; j < outside_spores[i].size(); j++) {
+ int row = std::get<0>(outside_spores[i][j]);
+ int col = std::get<1>(outside_spores[i][j]);
+ double n = Rast_row_to_northing(row, ®ion);
+ double e = Rast_col_to_easting(col, ®ion);
+ Vect_reset_line(Points);
+ Vect_reset_cats(Cats);
+ Vect_append_point(Points, e, n, 0);
+ Vect_cat_set(Cats, 1, i + 1);
+ Vect_write_line(&Map, GV_POINT, Points, Cats);
+ }
+ }
+ Vect_build(&Map);
+ Vect_close(&Map);
+ Vect_destroy_line_struct(Points);
+ Vect_destroy_cats_struct(Cats);
+ }
+
+ // clean the allocated memory
+ if (weather) {
+ delete[] mcf;
+ delete[] ccf;
+ delete[] weather;
+ }
+
+ return 0;
+}
Property changes on: grass-addons/grass7/raster/r.spread.sod/main.cpp
___________________________________________________________________
Added: svn:mime-type
+ text/x-c++src
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.spread.sod/r.spread.sod.html
===================================================================
--- grass-addons/grass7/raster/r.spread.sod/r.spread.sod.html (rev 0)
+++ grass-addons/grass7/raster/r.spread.sod/r.spread.sod.html 2017-09-25 02:17:49 UTC (rev 71505)
@@ -0,0 +1,77 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.spread.sod</em> simulates spread of a plant pathogen, specifically
+it simulates spread of Sudden Oak Death disease in California and Oregon
+(USA).
+
+
+<h2>NOTES</h2>
+
+The directions of wind consider north (N) to be grid north, if your
+true north is different direction, you need to make an adjustment.
+
+
+<h2>EXAMPLES</h2>
+
+Example of creating file with list of input maps (unix-like command
+line):
+
+<div class="code"><pre>
+g.list type=raster pattern="moisture_*" mapset=PRISM -m > moistures.txt
+g.list type=raster pattern="temperature_*" mapset=PRISM -m > temperatures.txt
+</pre></div>
+
+Note that the above assumes that the names will be ordered by time.
+You may need to pipe the result through <em>sort</em> with <tt>-n</tt>
+(<tt>| sort -n</tt>) to achieve this if filenames contain ordinal
+numbers but are not zero-padded (<tt>temperature_001</tt> versus
+<tt>temperature_1</tt>).
+
+Example of the run of the model (unix-like command line):
+
+<div class="code"><pre>
+r.spread.sod species=lide lvtree=all infected=infected_2005 \
+ moisture_file=moistures.txt temperature_file=temperatures.txt \
+ output=spread_sod start_time=2005 end_time=2010 \
+ rtype=cauchy wind=NE rseed=4
+</pre></div>
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>
+ Ross K. Meentemeyer, Nik J. Cunniffe, Alex R. Cook,
+ Joao A. N. Filipe, Richard D. Hunter, David M. Rizzo,
+ and Christopher A. Gilligan 2011.
+ Epidemiological modeling of invasion in heterogeneous landscapes:
+ spread of sudden oak death in California (1990-2030).
+ <em>Ecosphere</em> 2:art17.
+ <a href="http://dx.doi.org/10.1890/ES10-00192.1">DOI: 10.1890/ES10-00192.1</a>
+<li>
+ Tonini, Francesco, Douglas Shoemaker, Anna Petrasova, Brendan Harmon,
+ Vaclav Petras, Richard C. Cobb, Helena Mitasova,
+ and Ross K. Meentemeyer.
+ Tangible geospatial modeling for collaborative solutions
+ to invasive species management.
+ <em>Environmental Modelling & Software</em> 92 (2017): 176-188.
+ <a href="https://doi.org/10.1016/j.envsoft.2017.02.020">DOI: 10.1016/j.envsoft.2017.02.020</a>
+</ul>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+ <a href="r.spread.html">r.spread</a>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Francesco Tonini* (original R version)<br>
+Zexi Chen* (C++ version)<br>
+Vaclav Petras* (parallelization, GRASS interface)<br>
+Anna Petrasova* (single species simulation)<br>
+
+<br>
+* <a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a>
Property changes on: grass-addons/grass7/raster/r.spread.sod/r.spread.sod.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list