SOD Spread tutorial: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(Old PoPS repo is PoPS Core now)
 
(12 intermediate revisions by one other user not shown)
Line 1: Line 1:
[https://grass.osgeo.org/grass74/manuals/addons/r.spread.sod.html r.spread.sod] is a model for stochastic landscape spread of the forest pathogen Sudden Oak Death.
{{AddonCmd|r.pops.spread}} is a model for stochastic landscape spread of the pest and pathogens. It uses [https://github.com/ncsu-landscape-dynamics/pops-core PoPS] (Pest or Pathogen Spread) Core library. In this tutorial we will use it specifically to model the spread of Sudden Oak Death tree disease.
This tutorial shows how to prepare data and run the model. You can download the [http://fatra.cnr.ncsu.edu/pops/data/Oregon.zip sample dataset] and simulate sudden oak death spread in the [https://goo.gl/maps/b1ZPn3abtW82 Rouge River-Siskiyou National Forest region] of western Oregon.
This tutorial shows how to run the model. You can download the [http://fatra.cnr.ncsu.edu/pops/data/PoPS_SOD_tutorial.zip sample dataset] and simulate sudden oak death spread in the [https://goo.gl/maps/b1ZPn3abtW82 Rouge River-Siskiyou National Forest region] of western Oregon.


= Software =
= Software =
Required software includes:
Required software includes:
* GRASS GIS >=7.2
* GRASS GIS >=7.8
* Addon module:
* Addon module:
** {{AddonCmd|r.spread.sod}}
** {{AddonCmd|r.pops.spread}}


= Input data used in this tutorial =
= Input data used in this tutorial =
Download the [http://fatra.cnr.ncsu.edu/pops/data/Oregon.zip sample dataset] containing:
Download the [http://fatra.cnr.ncsu.edu/pops/data/PoPS_SOD_tutorial.zip sample dataset] containing:
* digital elevation model
* digital elevation model
* forest density
* orthophoto of study area
* orthophoto of study area
* study boundaries
* host layer
* mapset containing moisture and temperature values
* layer of all trees
* roads and rivers
* mapset containing weather coefficients
The sample dataset is a GRASS GIS Location, so it goes into your GRASS GIS Database which is usually a directory called grassdata in your home directory or your Documents directory.
The sample dataset is a GRASS GIS Location, so it goes into your GRASS GIS Database which is usually a directory called grassdata in your home directory or your Documents directory.


= Workflow =
= Workflow =
Download sample data and unzip it. Launch GRASS GIS and select the unzipped <tt>Oregon</tt> directory, Location <tt>Oregon</tt> and create a new Mapset <tt>tutorial</tt>.
Download sample data and unzip it into your grassdata. Launch GRASS GIS and select Location <tt>PoPS_SOD_tutorial</tt> and Mapset <tt>tutorial</tt>.


Change working directory:
Change working directory by
 
typing <tt>cd</tt> (stands for change directory) into the GUI Console and hit Enter. You can create a directory for this exercise (outside of grassdata) anywhere where you can write and read. Navigate to that directory.
Type <tt>cd</tt> (stands for change directory) into the GUI Console and hit Enter. Navigate to where ever you have saved the example data and change directory to the <tt>tutorial</tt> mapset.


Install addon:
Install addon:
<source lang="bash">
<source lang="bash">
g.extension r.spread.sod
g.extension r.pops.spread
</source>
</source>


First, we will set computational region of our analyses to an extent covering our study area and so that the cells are aligned with one of our rasters:
First, we will set computational region of our analyses to an extent covering our study area:
<source lang="bash">
<source lang="bash">
g.region vector=EU1_12x8 -p
g.region region=small_study_area -p
</source>
</source>


== Create files with lists of input maps ==
== Create files with lists of input maps ==
List all maps in the weather mapset in one text file using g.list. Search for ppt for precipitation and tmean or temperature mean for temperature:
Here we used already prepared weather coefficients, these can be created using this [https://github.com/ncsu-landscape-dynamics/weather-coefficient workflow].
Mapset weather includes raster layers per each week of a simulation. With that, we need to prepare a text file. List and write the maps in a file using g.list:
<source lang="bash">
g.list -m type=raster pattern="average_weather_*" mapset=weather output=weather.txt
</source>
The file should be written to your current working directory.
 
== Display infected data ==
Download [https://raw.githubusercontent.com/ncsu-landscape-dynamics/pops-intro-grass-notebook/master/color_infected.txt color ramp] into your current working directory to be used for better visualization of infected cells. Display the infected data. When displaying the infection, do not display zeros.
 
<source lang="bash">
<source lang="bash">
g.list -m type=raster pattern=*ppt* mapset=weather output=precipitation.txt
r.colors map=eu_infection_2019 rules=color_infected.txt
g.list -m type=raster pattern=*tmean* mapset=weather output=temperature.txt
d.rast map=ortho
d.rast map=eu_infection_2019 values=0 -i
d.vect map=NHDFlowline where="FCODE >= 46006" color=30:144:255
d.vect map=roads where="FULLNAME is not NULL" color=165:159:159 width=2
</source>
</source>


== Compute the spread of SOD for default values ==
== Compute the spread of SOD for default values ==
[[File:Spread sod.png|150px|thumb|Output of SOD spread for default r.spread.sod parameters for the years 2016-2021. Overlaid onto the orthophoto image.]]
Run the model using the text files created and setting only the required parameters. For this analysis we used a wind direction of NE and are looking at the first 5 years of spread from the initial 2019 infection discovery.
Run the model using the text files created and setting only the required parameters. For this analysis we used a wind direction of NE and are looking at the first 5 years of spread from the initial 2016 infection discovery.
 
We export a result from a single stochastic simulation (with specified random seed).
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_sod output_frequency=yearly runs=1 random_seed=1
</source>
 
We list newly created output layers representing infected trees in each year of the simulation and we set a custom color ramp.
 
<source lang="bash">
g.list type=raster pattern="spread_sod*" output=series.txt
r.colors rules=color_infected.txt file=series.txt
</source>
 
You can visualize the spread series in Animation Tool, do not display 0 values for infected maps.
 
[[File:SOD_spread_example.png|upright|center|Example output of SOD spread. Overlaid onto the orthophoto image.]]
 
We can run multiple stochastic runs and aggregate the results into a probability layer (0.1 if cell was infected once in 10 runs), average layer (average number of infected trees per cell) and standard deviation layer.
 
Here we run the process 10x and we use 2 cores for parallel processing (providing module was compiled with OpenMP support and cores are available).
 
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 average_series=average probability_series=probability stddev_series=stddev output_frequency=yearly runs=10 nprocs=2 random_seed=1
</source>
 
Set color ramp for probability:
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sod wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 random_seed=4
g.list type=raster pattern="probability*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
</source>


Add the initial infected area and ortho photo layers to the display to view the spread. Rearrange layers so ortho is on the bottom.
Explore results in Map Display or using Animation Tool, omit displaying 0 values of probability.
 
== Effect of dispersal kernel ==
The choice and parametrization of dispersal kernel significantly influences the spread and should be informed by calibration. The natural dispersal kernel (required) typically represents wind dispersal, additionally, we can optionally add the anthropogenic kernel which represents more human affected spread possibly over longer distances.
 
Each kernel is defined by type (cauchy, exponential), direction (none, N, NE, E, ...), direction strength (concentration around the direction using [https://en.wikipedia.org/wiki/Von_Mises_distribution von mises distribution]) and scale (distance).
 
<source lang="bash">
<source lang="bash">
d.rast inf_2016_null
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=400 natural_direction_strength=3 single_series=spread_distance output_frequency=yearly runs=1 random_seed=1
d.rast ortho
 
g.list type=raster pattern="spread_distance*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt
</source>
</source>


Set all spread values of 0 to null to view orthophoto layer underneath. Then add the layer.
Compare our initial run with run with increased kernel scale.
 
We can then also change the direction to E and increase the strength of direction:
 
<source lang="bash">
<source lang="bash">
r.null map=spread_sod setnull=0
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=E natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=10 single_series=spread_direction output_frequency=yearly runs=1 random_seed=1
d.rast spread_sod
 
g.list type=raster pattern="spread_direction*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt
</source>
</source>


[[File:Spread sodk1.png|150px|thumb|Output of SOD spread for kappa=0 for the years 2016-2021. Overlaid onto the orthophoto image.]]
Finally, we can select a different kernel type.
[[File:Spread sodk2.png|150px|thumb|Output of SOD spread for kappa=4 for the years 2016-2021. Overlaid onto the orthophoto image.]]
== Working with the optional parameters ==
===Kappa(k)===
A measure of concentration, a reciprocal measure of dispersion. For this model the Von Mises circular normal distribution equation was used. The default value of k is 2 so the effect wind dispersal has on the spores shows an average distribution.  


If k is set to 0 the wind dispersal will have no effect on distribution.
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sodk1 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 kappa=0 random_seed=4
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6  weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=cauchy natural_distance=242 natural_direction_strength=3 single_series=spread_type output_frequency=yearly runs=1 random_seed=1
 
g.list type=raster pattern="spread_type*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt
</source>
</source>


If k is larger than 2 the distribution becomes very concentrated around the prevailing wind direction.
== Effect of reproductive rate ==
Similarly to kernel, reproductive rate should be informed by calibration.
This model uses the [https://en.wikipedia.org/wiki/Poisson_distribution Poisson distribution] to generate spores.
In this example we double it:
 
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sodk2 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 kappa=4 random_seed=4
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=3 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_rate output_frequency=yearly runs=1 random_seed=1
 
g.list type=raster pattern="spread_rate*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt
</source>
</source>


===Gamma(y)===
== Treatments ==
[[File:Spread sodsr1.png|150px|thumb|Output of SOD spread for spore_rate=0 for the years 2016-2021. Overlaid onto the orthophoto image.]]
[[File:Spread sodsr2.png|150px|thumb|Output of SOD spread for spore_rate=10 for the years 2016-2021. Overlaid onto the orthophoto image.]]
A parameter between 0 and 1 used in the equation for Bernoulli distribution. This parameter is then factored into the Cauchy distribution model for the probability of using the first model. In this instance, a change in the gamma value doesn’t visually affect the output. An example console input for this factor would be:
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sodg1 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 gamma=0.5 random_seed=4
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=3 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_rate output_frequency=yearly runs=1 random_seed=1
 
g.list type=raster pattern="probability*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
</source>


===Spore_rate===
We treat the initial infection and a buffer around it. The treatments are applied at the end of the year.
The spore production rate per week for each infected tree. This model uses the Poisson distribution model to determine the spore rate value. At the time of the initial infection the spores were calculated to spread at a rate of 4.4 spores per week.
<source lang="bash">
r.buffer -z input=eu_infection_2019 output=buffer_A distances=200
r.mapcalc "treatment_A = if (isnull(buffer_A), 0, 1)"
</source>
 
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentA output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_A treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all
 
g.list type=raster pattern="probtreatmentA*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
 
Here we increase the buffer size:
<source lang="bash">
r.buffer -z input=eu_infection_2019 output=buffer_B distances=500
r.mapcalc "treatment_B = if (isnull(buffer_B), 0, 1)"
</source>
 
 
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentB output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_B treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all
 
g.list type=raster pattern="probtreatmentB*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
 
Now we create a large 1km wide barrier in an attempt to stop the spread. For this scenario, we assume the treatment is not 100% effective, but rather only 75% of host is removed.
 
<source lang="bash">
r.mapcalc "treatment_C = if (y() > 4687000 && y() < 4688000, 0.75, 0 )"
</source>
 
We will see that for the final year of our simulation, the disease spread in several stochastic runs over the barrier:
 
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentC output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_C treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all
 
g.list type=raster pattern="probtreatmentC*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
 
Finally, we manage with the 200m buffer treatments A in 2019 and with the barrier in 2021:
 
<source lang="bash">
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentAC output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_A,treatment_C treatment_date=2019-12-31,2021-12-31 treatment_length=0,0 treatment_application=ratio_to_all
 
g.list type=raster pattern="probtreatmentAC*" output=series.txt --o
r.colors color=magma file=series.txt
</source>
 
 
<!--
== Wind dispersal ==
=== Direction ===
Kappa influences the wind dispersal direction.  
Wind dispersal uses the [https://en.wikipedia.org/wiki/Von_Mises_distribution Von Mises] circular normal distribution.
Large kappa makes the distribution very concentrated about the angle, specifically the wind direction specified.
If kappa is set to 0 the wind dispersal direction will have no effect on distribution.
 
<source lang="bash">
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_k0 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week kappa=0 random_seed=4
</source>
 
If kappa is larger than 2 the distribution becomes very concentrated around the specified wind direction.
<source lang="bash">
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_k4 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week kappa=4 random_seed=4
</source>
 
=== Distance ===
Parameter short_distance_scale decides how far the spores can travel and it is the input to [https://en.wikipedia.org/wiki/Cauchy_distribution Cauchy distribution].
<source lang="bash">
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_sc5 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week short_distance_scale=5 random_seed=4
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_sc30 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week short_distance_scale=30 random_seed=4
</source>
 
 
== Spore rate ==
The spore production rate per week for each infected tree. This model uses the [https://en.wikipedia.org/wiki/Poisson_distribution Poisson distribution] to generate spores.


Decreasing the rate down to 1 scales down the spread a great deal. However there are still new areas of spread.
Decreasing the rate down to 1 scales down the spread a great deal. However there are still new areas of spread.
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sodsr1 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 spore_rate=1 random_seed=4
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_sr1 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week reproductive_rate=1 random_seed=4
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_sr6 wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week reproductive_rate=6 random_seed=4
</source>
</source>


Increasing the rate to 10 drastically increases the spread, as to be expected.
 
== Treatments ==
We will treat the area by removing the host. We will develop several scenarios.
First, no treatment is applied:
 
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sodsr2 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 spore_rate=10 random_seed=4
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_notreatment wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week  random_seed=4
</source>
</source>


Updating these calculations on a regular basis can be done quite easily through the model to better monitor the spread. These values can also be adjusted for distinct locations or different forms of sudden oak death.
We treat the initial infection and a buffer around it. The treatments are applied at the end of the year.
 
[[File:Spread sods1 1.png|150px|thumb|Output of SOD spread for scale_1=1 for the years 2016-2021. Overlaid onto the orthophoto image.]]
[[File:Spread sods1 2.png|150px|thumb|Output of SOD spread for scale_1=50 for the years 2016-2021. Overlaid onto the orthophoto image.]]
===Scale_1 and Scale_2===
Parameters dealing with the Cauchy distribution, which is the equation that calculates the continuous random distribution of the spores based on inputs. These scale factors effect the outcome by increasing or decreasing the distribution of the spores. Scale_1 is default set to 20.57 which was found to be the initial spread scale in the western Oregon SOD case.  


Setting scale_1 factor to 1 so it has no effect in the equation.
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sods1_1 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 scale_1=1 random_seed=4
r.buffer -z input=inf_2016 output=buffer_A distances=200
r.mapcalc "treatment_A = if (isnull(buffer_A), 0, 1)"
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_treatA wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week  treatments=treatment_A treatment_year=2016 random_seed=4
</source>
</source>


Setting scale_1 factor to 50 so it has a greater effect on the spread.
Here we increase the buffer size:
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sods1_2 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 scale_1=50 random_seed=4
r.buffer -z input=inf_2016 output=buffer_B distances=500
r.mapcalc "treatment_B = if (isnull(buffer_B), 0, 1)"
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_treatB wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week  treatments=treatment_B treatment_year=2016 random_seed=4
</source>
</source>


Scale_2 is undefined in the equation when the model is run on default. In this instance scale_2 has no visual effect on the output. An example console input for this factor would be:
Finally, we use the buffer from last example and create a large 1km wide barrier in an attempt to stop the spread. Inspect the annual results to see how effective the solution is:
<source lang="bash">
<source lang="bash">
r.spread.sod species=lide_den_int lvtree=all_den_int infected=inf_2016 output=spread_sods2 wind=NE moisture_file=precipitation.txt temperature_file=temperature.txt start_time=2016 end_time=2021 scale_2=20 random_seed=4
r.mapcalc "barrier = if (y() > 4684000 && y() <4685000, 1, 0 )"
r.pops.spread host=host total_plants=all_trees infected=inf_2016 output=spread_sod_treatB_barrier output_series=spread_sod_treatB_barrier wind=N moisture_coefficient_file=precipitation.txt temperature_coefficient_file=temperature.txt start_time=2016 end_time=2021 step=week  treatments=treatment_B,barrier treatment_year=2016,2018 random_seed=4
</source>
</source>
-->

Latest revision as of 20:35, 7 September 2020

r.pops.spread is a model for stochastic landscape spread of the pest and pathogens. It uses PoPS (Pest or Pathogen Spread) Core library. In this tutorial we will use it specifically to model the spread of Sudden Oak Death tree disease. This tutorial shows how to run the model. You can download the sample dataset and simulate sudden oak death spread in the Rouge River-Siskiyou National Forest region of western Oregon.

Software

Required software includes:

Input data used in this tutorial

Download the sample dataset containing:

  • digital elevation model
  • orthophoto of study area
  • host layer
  • layer of all trees
  • roads and rivers
  • mapset containing weather coefficients

The sample dataset is a GRASS GIS Location, so it goes into your GRASS GIS Database which is usually a directory called grassdata in your home directory or your Documents directory.

Workflow

Download sample data and unzip it into your grassdata. Launch GRASS GIS and select Location PoPS_SOD_tutorial and Mapset tutorial.

Change working directory by typing cd (stands for change directory) into the GUI Console and hit Enter. You can create a directory for this exercise (outside of grassdata) anywhere where you can write and read. Navigate to that directory.

Install addon:

g.extension r.pops.spread

First, we will set computational region of our analyses to an extent covering our study area:

g.region region=small_study_area -p

Create files with lists of input maps

Here we used already prepared weather coefficients, these can be created using this workflow. Mapset weather includes raster layers per each week of a simulation. With that, we need to prepare a text file. List and write the maps in a file using g.list:

g.list -m type=raster pattern="average_weather_*" mapset=weather output=weather.txt

The file should be written to your current working directory.

Display infected data

Download color ramp into your current working directory to be used for better visualization of infected cells. Display the infected data. When displaying the infection, do not display zeros.

r.colors map=eu_infection_2019 rules=color_infected.txt
d.rast map=ortho
d.rast map=eu_infection_2019 values=0 -i
d.vect map=NHDFlowline where="FCODE >= 46006" color=30:144:255
d.vect map=roads where="FULLNAME is not NULL" color=165:159:159 width=2

Compute the spread of SOD for default values

Run the model using the text files created and setting only the required parameters. For this analysis we used a wind direction of NE and are looking at the first 5 years of spread from the initial 2019 infection discovery.

We export a result from a single stochastic simulation (with specified random seed).

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_sod output_frequency=yearly runs=1 random_seed=1

We list newly created output layers representing infected trees in each year of the simulation and we set a custom color ramp.

g.list type=raster pattern="spread_sod*" output=series.txt
r.colors rules=color_infected.txt file=series.txt

You can visualize the spread series in Animation Tool, do not display 0 values for infected maps.

Example output of SOD spread. Overlaid onto the orthophoto image.
Example output of SOD spread. Overlaid onto the orthophoto image.

We can run multiple stochastic runs and aggregate the results into a probability layer (0.1 if cell was infected once in 10 runs), average layer (average number of infected trees per cell) and standard deviation layer.

Here we run the process 10x and we use 2 cores for parallel processing (providing module was compiled with OpenMP support and cores are available).

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 average_series=average probability_series=probability stddev_series=stddev output_frequency=yearly runs=10 nprocs=2 random_seed=1

Set color ramp for probability:

g.list type=raster pattern="probability*" output=series.txt --o
r.colors color=magma file=series.txt

Explore results in Map Display or using Animation Tool, omit displaying 0 values of probability.

Effect of dispersal kernel

The choice and parametrization of dispersal kernel significantly influences the spread and should be informed by calibration. The natural dispersal kernel (required) typically represents wind dispersal, additionally, we can optionally add the anthropogenic kernel which represents more human affected spread possibly over longer distances.

Each kernel is defined by type (cauchy, exponential), direction (none, N, NE, E, ...), direction strength (concentration around the direction using von mises distribution) and scale (distance).

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=400 natural_direction_strength=3 single_series=spread_distance output_frequency=yearly runs=1 random_seed=1

g.list type=raster pattern="spread_distance*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt

Compare our initial run with run with increased kernel scale.

We can then also change the direction to E and increase the strength of direction:

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=E natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=10 single_series=spread_direction output_frequency=yearly runs=1 random_seed=1

g.list type=raster pattern="spread_direction*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt

Finally, we can select a different kernel type.

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6  weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=cauchy natural_distance=242 natural_direction_strength=3 single_series=spread_type output_frequency=yearly runs=1 random_seed=1

g.list type=raster pattern="spread_type*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt

Effect of reproductive rate

Similarly to kernel, reproductive rate should be informed by calibration. This model uses the Poisson distribution to generate spores. In this example we double it:

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=3 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_rate output_frequency=yearly runs=1 random_seed=1

g.list type=raster pattern="spread_rate*" output=series.txt --o
r.colors rules=color_infected.txt file=series.txt

Treatments

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=3 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 single_series=spread_rate output_frequency=yearly runs=1 random_seed=1

g.list type=raster pattern="probability*" output=series.txt --o
r.colors color=magma file=series.txt

We treat the initial infection and a buffer around it. The treatments are applied at the end of the year.

r.buffer -z input=eu_infection_2019 output=buffer_A distances=200
r.mapcalc "treatment_A = if (isnull(buffer_A), 0, 1)"
r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentA output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_A treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all

g.list type=raster pattern="probtreatmentA*" output=series.txt --o
r.colors color=magma file=series.txt

Here we increase the buffer size:

r.buffer -z input=eu_infection_2019 output=buffer_B distances=500
r.mapcalc "treatment_B = if (isnull(buffer_B), 0, 1)"


r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentB output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_B treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all

g.list type=raster pattern="probtreatmentB*" output=series.txt --o
r.colors color=magma file=series.txt

Now we create a large 1km wide barrier in an attempt to stop the spread. For this scenario, we assume the treatment is not 100% effective, but rather only 75% of host is removed.

r.mapcalc "treatment_C = if (y() > 4687000 && y() < 4688000, 0.75, 0 )"

We will see that for the final year of our simulation, the disease spread in several stochastic runs over the barrier:

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentC output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_C treatment_date=2019-12-31 treatment_length=0 treatment_application=ratio_to_all

g.list type=raster pattern="probtreatmentC*" output=series.txt --o
r.colors color=magma file=series.txt

Finally, we manage with the 200m buffer treatments A in 2019 and with the barrier in 2021:

r.pops.spread host=host total_plants=max_host infected=eu_infection_2019 start_date=2019-01-01 end_date=2023-12-31 step_unit=week reproductive_rate=1.6 weather_coefficient_file=weather.txt natural_direction=NE natural_dispersal_kernel=exponential natural_distance=242 natural_direction_strength=3 probability_series=probtreatmentAC output_frequency=yearly runs=10 nprocs=2 random_seed=1 treatments=treatment_A,treatment_C treatment_date=2019-12-31,2021-12-31 treatment_length=0,0 treatment_application=ratio_to_all

g.list type=raster pattern="probtreatmentAC*" output=series.txt --o
r.colors color=magma file=series.txt