Displaying properties of spatial objects

We begin by loading an example dataset of 100 forest stands distributed on points in the landscape:

data("examplepointslandscape")

We transform the data set to class sf because sp objects are deprecated from ver. 1.0.0 on.

pts_sf <- sp_to_sf(examplepointslandscape)
pts_sf
## Simple feature collection with 100 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 402000 ymin: 4640000 xmax: 429000 ymax: 4646000
## Projected CRS: ETRS89 / UTM zone 31N
## # A tibble: 100 × 9
##            geometry id    elevation slope aspect land_cover_type forest   soil  
##         <POINT [m]> <chr>     <dbl> <dbl>  <dbl> <chr>           <list>   <list>
##  1 (411000 4646000) 81065       516  9.94   92.7 wildland        <forest> <soil>
##  2 (412000 4646000) 81066       537 13.7   115.  wildland        <forest> <soil>
##  3 (413000 4646000) 81067       578 23.3   152.  wildland        <forest> <soil>
##  4 (415000 4646000) 81068       615 12.7   360.  wildland        <forest> <soil>
##  5 (416000 4646000) 81069       644 22.5   186.  wildland        <forest> <soil>
##  6 (423000 4646000) 81070       520  4.04  225.  wildland        <forest> <soil>
##  7 (426000 4646000) 81073       608 12.2   337.  wildland        <forest> <soil>
##  8 (428000 4646000) 81075       745  8.06  332.  wildland        <forest> <soil>
##  9 (402000 4645000) 81076       527 14.4   245.  wildland        <forest> <soil>
## 10 (404000 4645000) 81077       483 26.0    36.7 wildland        <forest> <soil>
## # … with 90 more rows, and 1 more variable: state <list>

Using plot_variable() functions for spatial landscape objects, we can draw maps of some variables using:

plot_variable(pts_sf, "basal_area")

The set of maps available can be known by inspecting the help of function extract_variables(). Alternatively, the package provides function shinyplot_land() to display maps interactively.

Simulations on a set of forest stands

Package medfate includes functions spwb(), growth() and fordyn() to simulate soil water balance, carbon balance and forest dynamics on a single forest stand, respectively. This section describe how to run simulations on a set of forest stands in one call. This is done using functions spwb_spatial(), growth_spatial() and fordyn_spatial(), respectively.

As an example, we will use function spwb_spatial(), which simulates soil plant water balance on forests distributed in particular locations. The function takes an object of class sf as input.

Weather inputs and local control parameters

As before, we need meteorological data, species parameters and control parameters for local simulations (which will apply to all forest stands):

data("examplemeteo")
data("SpParamsMED")
local_control <- defaultControl()

Calling the simulation function

The call to spwb_spatial() can be done as follows (here we restrict the dates for simplicity):

dates <- seq(as.Date("2001-01-01"), as.Date("2001-02-28"), by="day")
res <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo, 
                  dates = dates, local_control = local_control,  progress = FALSE)

Function spwb_spatial() first initializes model inputs by calling forest2spwbInput() for each forest stand described in the sf landscape object. Then it calls function spwb() for each forest stand and stores the result. The fact that we used examplemeteo as meteorological input involves that the same weather was applied to all forest stands, but different weather could have been specified for each one (see documentation of function spwb_spatial()).

The result of calling spwb_spatial() is an object of class sf with the following columns:

names(res)
## [1] "geometry" "id"       "state"    "result"

Column geometry contains the geometry given as input to simulations, column id contains the identification label of each stand, column state contains the spwbInput corresponding to each forest stand (which can be used in subsequent simulations) and column result contains the output of spwb() function for each forest stand (i.e. its elements are objects of the S3 class spwb).

Temporal summaries, plots and maps

The structure of the output of spwb_spatial() allows querying information for the simulation of any particular forest stand. For example, we may use function plot.spwb(), from package medfate, to display the simulation results on a particular plot:

plot(res$result[[1]], "Evapotranspiration")

Similarly, if we want a monthly summary of water balance for the first stand, we can use function summary.spwb() from package medfate:

summary(res$result[[1]], freq="months",FUN=sum, output="WaterBalance")
##                 PET Precipitation      Rain      Snow    NetRain  Snowmelt
## 2001-01-01 45.07945      74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521       4.99943  2.457859  2.541571  0.9699644  5.552842
##            Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01    52.730084      0   40.6162567           34.03849    18.461760
## 2001-02-01     6.522806      0    0.1569381           21.18700     1.487895
##            SoilEvaporation PlantExtraction Transpiration
## 2001-01-01       5.5785170        9.996655      9.998214
## 2001-02-01       0.9530808       18.740241     18.746027
##            HydraulicRedistribution
## 2001-01-01              0.03169547
## 2001-02-01              0.00000000

However, a more convenient way of generating summaries is by calculating them on all forest stands in one step, using function simulation_summary() on objects issued from simulations:

res_sum <- simulation_summary(res, summary_function = summary.spwb, 
                              freq="months", output="WaterBalance")

where the argument summary_function points to the function to be used to generate local summaries and the remaining arguments are those of the local summary function. The result of using simulation_summary() is again an object of class sf that contains the spatial geometry and the list of summaries for all stands:

names(res_sum)
## [1] "geometry" "id"       "summary"

The summary for the first stand can now be accessed through the first element of column summary:

res_sum$summary[[1]]
##                 PET Precipitation      Rain      Snow    NetRain  Snowmelt
## 2001-01-01 45.07945      74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521       4.99943  2.457859  2.541571  0.9699644  5.552842
##            Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01    52.730084      0   40.6162567           34.03849    18.461760
## 2001-02-01     6.522806      0    0.1569381           21.18700     1.487895
##            SoilEvaporation PlantExtraction Transpiration
## 2001-01-01       5.5785170        9.996655      9.998214
## 2001-02-01       0.9530808       18.740241     18.746027
##            HydraulicRedistribution
## 2001-01-01              0.03169547
## 2001-02-01              0.00000000

Summary objects are handy because their plot_summary() function allows us to display maps of summaries for specific dates:

plot_summary(res_sum, "Transpiration", "2001-01-01", limits=c(0,45))

plot_summary(res_sum, "Transpiration", "2001-02-01", limits=c(0,45))

To avoid displaying maps one by one, the package includes function shinyplot_land() that allows displaying maps of temporal summaries interactively.

Advanced simulation features

Simulation with integrated temporal summaries

If one needs to save memory, it is possible with spwb_spatial() to generate temporal summaries automatically after the simulation of soil water balance of each stand, and storing those summaries instead of all the output of function spwb().

For example the following code will keep temporal summaries of water balance components instead of simulation results:

res_2 <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo, 
                  dates = dates, local_control = local_control,                  
                  keep_results = FALSE, progress = FALSE,
                  summary_function = summary.spwb, summary_arguments = list(freq="months"))

Parameter keep_results = FALSE tells spwb_spatial() not to keep the simulation results of forest stands, whereas summary_function = summary.spwb tells spwb_spatial() to perform and store summaries before discarding the results of any stand. The output has slightly different column names:

names(res_2)
## [1] "geometry" "id"       "state"    "result"   "summary"

In particular, result is not included. Now the temporal summaries can be directly accessed through the column summary:

res_2$summary[[1]]
##                 PET Precipitation      Rain      Snow    NetRain  Snowmelt
## 2001-01-01 45.07945      74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521       4.99943  2.457859  2.541571  0.9699644  5.552842
##            Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01    52.730084      0   40.6162567           34.03849    18.461760
## 2001-02-01     6.522806      0    0.1569381           21.18700     1.487895
##            SoilEvaporation PlantExtraction Transpiration
## 2001-01-01       5.5785170        9.996655      9.998214
## 2001-02-01       0.9530808       18.740241     18.746027
##            HydraulicRedistribution
## 2001-01-01              0.03169547
## 2001-02-01              0.00000000

And one can produce maps with summary results directly from the output of the simulation function:

plot_summary(res_2, "Transpiration", "2001-02-01", limits=c(0,45))

The possibility of running a summary function after the simulation of each stand is not limited to summary.spwb(). Users can define their own summary functions, provided the first argument is object, which will contain the result of the simulation (i.e., the result of calling spwb(), growth() or fordyn()). For example, the following function calls droughtStress():

f_stress <- function(object, ...) {
  return(droughtStress(object, ..., draw = FALSE))
}

Now we can call again spwb_spatial:

res_3 <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo, 
                  dates = dates, local_control = local_control,                  
                  keep_results = FALSE, progress = FALSE,
                  summary_function = f_stress, 
                  summary_arguments = list(index = "ADS", freq = "months", bySpecies=TRUE))

The drought stress summary of stand #1 is:

res_3$summary[[1]]
##            Dorycnium spp. Genista spp. Hedera helix Juniperus communis
## 2001-01-01            NaN          NaN  0.006165402        0.003861844
## 2001-02-01            NaN          NaN  0.007618106        0.004739411
##            Lavandula spp. Phillyrea angustifolia Pinus nigra Pistacia lentiscus
## 2001-01-01     0.01720418            0.003735971 0.009425589        0.003528457
## 2001-02-01     0.02141699            0.004658094 0.011569840        0.004327903
##            Quercus coccifera Quercus faginea Quercus ilex Rubus spp.
## 2001-01-01       0.003229448             NaN   0.01092188        NaN
## 2001-02-01       0.004017091             NaN   0.01272611        NaN
##            Salvia rosmarinus Viburnum spp.
## 2001-01-01               NaN           NaN
## 2001-02-01               NaN           NaN

And we can draw a map with the monthly drought stress of Quercus ilex across space using:

plot_summary(res_3, "Quercus ilex", "2001-02-01")

Simulation in several steps

The result of a simulation includes an element state, which stores the state of soil and stand variables at the end of the simulation. This information can be used to perform a new simulation from the point where the first one ended. In order to do so, we need to update the state variables in spatial object with their values at the end of the simulation, using function update_landscape():

pts_sf_mod <- update_landscape(pts_sf, res)

The resulting object w is the same as y except for the state variables. For example we can compare the water potential in the first soil layer:

plot_variable(pts_sf, "psi1")

plot_variable(pts_sf_mod, "psi1")

By using w as input we can now simulate water balance in the set of stands for an extra month:

dates <- seq(as.Date("2001-03-01"), as.Date("2001-03-31"), by="day")
res_3 <- spwb_spatial(pts_sf_mod, SpParamsMED, examplemeteo, 
                      dates = dates, local_control = local_control, 
                      summary_function = summary.spwb, summary_arguments = list(freq = "months"),
                      progress = FALSE)

And display a map with the resulting month transpiration:

plot_summary(res_3, "Transpiration", "2001-03-01", limits=c(0,45))

Simulation of watershed hydrology

Simulation of watershed hydrology involves describing a watershed in a gridded mode and simulating soil water balance water while including additional compartments and hydrological processes, such as routing runoff from one cell to the other.

Example data set

To illustrate this kind of simulation, which is carried out using function spwb_land(), we first load a small example watershed included with the package, which we transform to sf:

data("examplewatershed")
ws_sf <- sp_to_sf(examplewatershed)
ws_sf
## Simple feature collection with 66 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 401380 ymin: 4671820 xmax: 402880 ymax: 4672620
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 19
##            geometry    id elevation slope aspect land_cover_type forest   soil  
##         <POINT [m]> <int>     <dbl> <dbl>  <dbl> <chr>           <list>   <list>
##  1 (402630 4672570)     1      1162 11.3   79.2  wildland        <forest> <soil>
##  2 (402330 4672470)     2      1214 12.4   98.7  agriculture     <forest> <soil>
##  3 (402430 4672470)     3      1197 10.4  102.   wildland        <forest> <soil>
##  4 (402530 4672470)     4      1180  8.12  83.3  wildland        <forest> <soil>
##  5 (402630 4672470)     5      1164 13.9   96.8  wildland        <forest> <soil>
##  6 (402730 4672470)     6      1146 11.2    8.47 agriculture     <forest> <soil>
##  7 (402830 4672470)     7      1153  9.26 356.   agriculture     <forest> <soil>
##  8 (402230 4672370)     8      1237 14.5   75.1  wildland        <forest> <soil>
##  9 (402330 4672370)     9      1213 13.2   78.7  wildland        <forest> <soil>
## 10 (402430 4672370)    10      1198  8.56  75.6  agriculture     <forest> <soil>
## # … with 56 more rows, and 11 more variables: state <list>, waterOrder <int>,
## #   waterQ <list>, queenNeigh <list>, channel <lgl>, depth_to_bedrock <dbl>,
## #   bedrock_conductivity <dbl>, bedrock_porosity <dbl>, snowpack <dbl>,
## #   aquifer <dbl>, represented_area <dbl>

Here is a map of elevation values:

plot_variable(ws_sf, "elevation")

Since the landscape contains agricultural lands, we need to define crop factors, which will determine transpiration flow as a proportion of potential evapotranspiration:

ws_sf$crop_factor = NA
ws_sf$crop_factor[ws_sf$land_cover_type=="agriculture"] = 0.75

Simulation function

As for the call to function spwb_spatial() we will use the same weather (i.e. examplemeteo) across the watershed. To speed up calculations we call function spwb_land() for a single month.

dates <- seq(as.Date("2001-01-01"), as.Date("2001-01-31"), by="day")
res_ws <- spwb_land(ws_sf, SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month")
## 
## ------------ spwb_land ------------
## Grid cells: 66, patchsize: 100 m2, area: 0.66 ha
## Cell land use wildland: 48 agriculture: 17 artificial: 0 rock: 1 water: 0
## Cells with soil: 65
## Meteorological input class: data.frame
## Number of days to simulate: 31
## Number of summaries: 1
## Number of outlet cells: 1
## 
## Preparing spwb input...
## done.
## 
## Initial average soil water content (mm): 530.8
## Initial average snowpack water content (mm): 0
## Initial average aquifer water content (mm): 0
## Initial watershed water content (mm): 522.76
## 
## Performing daily simulations:
## .+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++done
## 
## Final average soil water content (mm): 530.8
## Final average snowpack water content (mm): 0.05
## Final average aquifer water content (mm): 53.17
## Final watershed water content (mm): 522.76
## 
## Change in snowpack water content (mm): 0.05
## Snowpack water balance result (mm): 3.19
## Snowpack water balance components:
##   Snow fall (mm) 16.65 Snow melt (mm) 13.46
## 
## Change in soil water content (mm): 0
## Soil water balance result (mm): -10.58
## Soil water balance components:
##   Net rainfall (mm) 37.92 Snow melt (mm) 13.47
##   Runon (mm) 1.1 Runoff (mm) 0
##   Subsurface input (mm) 0  Subsurface output (mm) 0
##   Deep drainage (mm) 53.99 Aquifer discharge (mm) 0
##   Soil evaporation (mm) 5.79 Plant transpiration (mm) 3.29
## 
## Change in aquifer water content (mm): 53.17
## Aquifer water balance result (mm): 53.17
## Aquifer water balance components:
##   Deep drainage (mm) 53.17 Aquifer discharge (mm) 0
## 
## Change in watershed water content (mm): 0
## Watershed water balance result (mm): 45.94
## Watershed water balance components:
##   Precipitation (mm) 74.75
##   Interception (mm) 19.87 Soil evaporation (mm) 5.7 Plant Transpiration (mm) 3.24
##   Export (mm) 0
## Watershed lateral flows:
##   Subsurface flow (mm) 0
##   Groundwater flow (mm) 7.77
## 
## ------------ spwb_land ------------

Although simulations are performed using daily temporal steps, parameter summary_frequency allows storing results at coarser temporal scales, to reduce the amount of memory in spatial results.

Summaries and plots

Unlike spwb_spatial() where summaries could be arbitrarily generated a posteriori from simulation results, with spwb_land() the summaries are always fixed and embedded with the simulation result. For example, we can inspect the summaries for a given landscape cell using:

res_ws$sf$summary[[1]]
##            Runon Runoff Infiltration     Rain NetRain     Snow Snowmelt
## 2001-01-01     0      0     66.83459 58.09884 53.1979 16.65065  13.6367
##            Interception DeepDrainage AquiferDischarge SaturationExcess
## 2001-01-01     4.900943     57.47486                0                0
##            SoilEvaporation Transpiration SWE  SoilVol        Psi1  WTD      DTA
## 2001-01-01        6.995272       3.28633   0 532.3416 -0.03300001 4000 14.78647
##            InterflowInput InterflowOutput BaseflowInput BaseflowOutput
## 2001-01-01              0               0      5.935362       4.507225

Several plots can be drawn from the result of function spwb_land() in a similar way as done for spwb_spatial(). As an example we display a map of soil water volume for the simulated month:

plot_summary(res_ws$sf, "SoilVol", "2001-01-01")