In MagicForrest/DGVMTools: DGVM Processing, Analysis and Plotting Tools

knitr::opts_chunk$set(echo=FALSE, fig.width=8, out.width = "100%", autodep=TRUE)
library(DGVMTools, quietly = TRUE, warn.conflicts = FALSE)

Introduction

DGVMTools is designed to provide core functionality for developing and utilising DGVMs (Dynamic Global Vegetation Models) and similar models (such as Land Surface Models). The basic motivation is that many people work with such models but there is no centralised effort to provide convenient, reliable and efficient tools to perform common tasks such as reading model data, aggregating it in space, time and across vegetation types, producing plots, comparing to standard observed datasets etc. Instead, people develop (and possibly share) scripts on an ad hoc basis, often resulting in duplicated effort and technical limitations hindering scientific analysis. DGVMTools aims to provide a common toolbox to enable such tasks.

Tutorial

Now we'll jump right in and start using the package to read, process and plot some DGVM output data and hopefully demonstrate the features, logic and power of the package. This tutorial uses small LPJ-GUESS output files included in the package. There is a more detailed description of the package (including implementation details) after the tutorial.

Step 1: Define a Source

The first job in pretty much any analysis is to define the source (or sources) of the data to be analysed. Here we use the defineSource() function to do exactly that, and it produces an object of type Source. To define a Source we need three key bits of information and these are the three essential arguments to defineSource() [^1]. These are:

1. the id argument, a unique identifier to label this source (should be alphanumeric with no spaces)
2. the dir argument, to specify the location of the data on disk (for the purposes of this tutorial we generate this with a function, but it is normally just a character string)
3. the format argument, which specifies the structure of the data on disk (ie what model did it come from or is it netCDF data)

The name argument is optional and can be used to provide a more human-friendly string to identify the run. If it is not provided, the id argument is used for the name. It is also possible just provide a name argument and not an id argument, in which case the id will be derived from the name.

# a little bit of magic to get the path to the example data included in the package on your system
example.run.directory <- system.file("extdata", "LPJ-GUESS_Runs", "CentralEurope", package = "DGVMTools")
print(example.run.directory)

GUESS.run <- defineSource(id = "LPJ-GUESS_Example",
                          dir = example.run.directory, # this would normally just be a character string containing a path
                          format = GUESS,
                          name = "LPJ-GUESS Example Run")

We can check the object's class take a look at it using the handy print function.

``` {r Source info, echo=TRUE} class(GUESS.run) print(GUESS.run)

That looks a little complicated, but it is mostly just meta-data designed to make things more convenient later on (or are needed for special cases such as grid cell coordinates that don't correspond to the centre of the grid cell), so we can ignore these for just now. Aside from the three obligatory things mentioned above,  the only thing that is information about the PFTs in the run, and even that is usually defined automatically in the `format` argument.  Note that at this stage we have only created a metadata object.  We actually retrieve data in the next step...

[^1]:  Experienced R users may have noticed that here we are using a specific function to make out `Source` object as opposed to `new()` function.  Whilst it is possible to construct a `Source` using `new()` it is __strongly recommended not to do so__ because using `defineSource()` is safer as it ensures the integrity and internal consistency of the `Source`.  It is also just _very_ much more convenient.

### Step 2: Reading data

To actually read data from disk with use the `getField()` function.  This can be as simple as the following:

``` {r Get data, echo=TRUE}
LAI.full <- getField(source = GUESS.run, 
                     quant = "lai")

This step read all the LAI data (across all available space and time) from the source and stores it as a Field object. We can have look at this new object using a print() command.

``` {r Field info, echo=TRUE} print(LAI.full)

This is a lot output but it basically describes everything that you might ever need to know about this `Field`.  The most most important part is the data, which is stored in what looks like an R `data.frame` (actually it is a `data.table` which is like a `data.frame` in many ways but often much faster) and each 'layer' (in this case each PFT) is stored as a column.  There is a bunch of meta-data included too, mostly about the spatial, temporal and annual dimensions of the data.  An important piece of metadata in the information about the type of data represented, this is stored in a `Quantity` object.  The `Quantity` was looked up automatically from the `quant = "lai"` argument.  

###  Step 3: Aggregating Data

Whilst reading in the entire data set may be useful, it is often to large and unwieldy to work with or to plot.  When analysing DGVM output, one commonly aggregates the data across space or time. To get the data in a more manageable form we are going to call `getField()` two more times, this time making use of its built in  __aggregate__ options.  First we open the file and simultaneously aggregate across all years using the `year.aggregate.method` (in this case we take the mean, but there are other options such as the sum or the variance)

``` {r Get data with year aggregation, echo=TRUE}
LAI.year.mean <- getField(source = GUESS.run, 
                          quant = "lai", 
                          year.aggregate.method = "mean")

To see how this is different to the previous case, we use the getDimInfo.

``` {r Check dimensions, echo = TRUE} getDimInfo(LAI.full) getDimInfo(LAI.year.mean)

We can see that the full data (`LAI.full`) has Lon, Lat and Year dimensions, but the yearly aggregated data (`LAI.year.mean`) has only Lon and Lat, since Year has been averaged away. 

We can now do a similar thing and aggregate across _space_.


``` {r Get data with spatial aggregation, echo = TRUE}
LAI.spatial.mean <- getField(source = GUESS.run, 
                             quant = "lai", 
                             spatial.aggregate.method = "mean")
getDimInfo(LAI.spatial.mean)

Note that this time we have aggregated over space, so the Lon and Lat dimensions have gone, but Year remains.

Step 4: Plotting data

Having got some sensibly aggregated data, we can now plot it using plotTemporal and plotSpatial.

Spatial plots (maps)

First we us plotSpatial to produce a map from the annually-averaged data. Note the underlying plot system is ggplot and plotSpatial returns a ggplot2 object which we then need to print in order to actually see it.

``` {r Spatial plot, fig.asp = 1, echo = TRUE} print(plotSpatial(LAI.year.mean))

save this plot with bigger text as an example plot for the reference publication, also omittingthe tropical layers which are all zero in europe

spatial_example_plot <- plotSpatial(LAI.year.mean, layers = c("BNE", "BINE", "BNS", "IBS", "TeBS", "TeBE", "C3G", "C4G", "Total"), text.multiplier = 3.5)

By default all layers are plotted.  We can choose one (or more) using the `layers` argument.
``` {r Spatial plot layers, fig.asp = 0.5, echo = TRUE}
print(plotSpatial(LAI.year.mean, layers = c("TeBS", "Total")))

The plotSpatial function comes with a lot of arguments which provide a great deal of flexibility. Also, since the the function returns a ggplot object, they can be further modified using standard ggplot commands, which also adds piles of flexibility. There are more examples of this later on.

Temporal plots (time series)

The spatially aggregated Field can be plotted as a time series using plotTemporal.

``` {r Temporal plot, fig.asp = 1, echo = TRUE} print(plotTemporal(LAI.spatial.mean))

use the long

save a temporal example plot for the reference publication

temporal_example_plot <- plotTemporal(LAI.spatial.mean, col.by = NULL, text.multiplier = 3.5, layers = c("BNE", "BINE", "BNS", "IBS", "TeBS", "TeBE", "C3G", "Total"))

And so we see a pretty straightforward times series plot of LAI averaged over the simulated area.

Note that trying to make plots of data with the incorrect dimensions should fail gracefully:

``` {r Fail plot, fig.asp = 1, echo = TRUE}
attempted.spatial.plot <- plotSpatial(LAI.spatial.mean)
attempted.temporal.plot <- plotTemporal(LAI.year.mean)

Analysis

Having read, aggregated and visualised data, the next steps would likely be further analysis such as aggregating different layers (eg combining all tree PFTs), deriving biomes and comparing to data. This section details such operations.

Combining layers using layerOp()

The layerOp() function enables operations on layers, including arithmetic operations (+,-,*,/) and statistical calculations (mean, max, var). Furthermore the . (dot) notation allows layers to be selected based on the properties defined Layer objects. For example here we select all layers when have the string "Tree" defined in any slot of their Layer object:

``` {r layerOp, echo=TRUE}

calculate tree total the long way around

LAI.year.mean <- layerOp(x = LAI.year.mean, operator = "+", layers = c("BINE", "BNE", "BNS", "IBS", "TeBS", "TeBE", "TeNE", "TrBE", "TrIBE", "TrBR"), new.layer = "TreeTotal1")

calculate tree total using dot notation

LAI.year.mean <- layerOp(x = LAI.year.mean, operator = "+", layers = ".Tree", new.layer = "TreeTotal2")

check the difference (subtract, print and plot)

LAI.year.mean <- layerOp(x = LAI.year.mean, operator = "-", layers = c("TreeTotal1", "TreeTotal2"), new.layer = "TreeTotalDiff") print(LAI.year.mean) print(plotSpatial(LAI.year.mean, c("TreeTotal1", "TreeTotal2", "TreeTotalDiff")))

Note that to use this shorthand, the layers must be defined and present in the `defined.layers` slot of the `Source` object.  Selections can also be made on the phenology, climate zone and leaf form of PFT layers.  The shorthand `".PFT"` can also be used to select all PFT layers.  To calculate the maximum PFT in each gridcell:

``` {r layerOp max, echo=TRUE}

# calculate maximum PFT 
LAI.year.mean <- layerOp(LAI.year.mean, "max.layer", ".PFT", "MaxPFT")

# calculate maximum  Tree
LAI.year.mean <- layerOp(LAI.year.mean, "max.layer", ".Tree", "MaxTree")

# plot
print(plotSpatial(LAI.year.mean, c("MaxPFT", "MaxTree")))

# calculate total tree and grass totals.  Note that further to `"+"` as used above, addition of layers can be defined with the with character string `"sum"` or and arbitary function i.e. `sum`
LAI.year.mean <- layerOp(LAI.year.mean, "sum", ".Grass", "Grass")
LAI.year.mean <- layerOp(LAI.year.mean, sum, ".Tree", "Tree")


# Calculate if trees or grasses are dominant 
LAI.year.mean <- layerOp(LAI.year.mean, "max.layer", c("Tree", "Grass"), "MaxLifeform")

# plot
print(plotSpatial(LAI.year.mean, c("MaxLifeform")))

Fractions can also be useful:

``` {r layerOp fraction, echo=TRUE}

calculate maximum PFT

LAI.year.mean <- layerOp(LAI.year.mean, "/", c("Tree", "Total"), "TreeFraction") LAI.year.mean <- layerOp(LAI.year.mean, "/", c("Grass", "Total"), "GrassFraction")

plot

print(plotSpatial(LAI.year.mean, c("TreeFraction", "GrassFraction")))

The `layerOp()` functions is very flexible (it even works with suitable arbitrary functions), see the function documentation for further details.  Note that the "-" and "/" operations only work with exactly two input layers given in the `layers` argument, and like the rest of the package they follow a "one-two" convention for ordering ie. "layer one minus layer two" or "layer one over layer two".

#### Reading data and comparing layers

The package can also read netCDF files using the `NetCDF` format.  As netCDF files can be structured in myriad different ways, it is not absolutely guaranteed that all files will read correctly.  However, most files that conform closely to the CF convention should work.  There are a couple of small example netCDF files included in the package, one of which we will use here and compare to another run of LPJ-GUESS, this time cover a part of central-western Africa.  Please contact the author if you have a request to read netCDF with particular features which aren't currently supported.

``` {r open data, fig.height=8, echo=TRUE}

# new Source (different model run, this time over Africa)
GUESS.Africa.run <- defineSource(id = "LPJ-GUESS_Example_Africa",
                                 dir = system.file("extdata", "LPJ-GUESS_Runs", "CentralAfrica", package = "DGVMTools"), 
                                 format = GUESS,
                                 name = "LPJ-GUESS Africa Example Run")

# get standard variable vegC_std between the years 2000 and 2010
GUESS.vegC <- getField(GUESS.Africa.run,
                       "vegC_std",
                       first.year = 2000,
                       last.year = 2010,
                       year.aggregate.method = "mean")

# calculate Tree total and have a plot
GUESS.vegC <- layerOp(GUESS.vegC, "+", ".Tree", "Tree")
print(plotSpatial(GUESS.vegC, "Tree"))


# define the data Source (Saatchi biomass data @ HD)
Saatchi.dataset <- defineSource(id = "Saatch2011",
                                name = "Saatchi et al. 2011 biomass",
                                format = NetCDF,
                                dir = system.file("extdata", "NetCDF", "Saatchi2011", "HD", package = "DGVMTools"))

# get the data field
Saatchi.vegC <- getField(source = Saatchi.dataset,
                         quant = "vegC_std")

# NOTE: the metadata from the netCDF file always over-writes the specified `Quantity` data with a warning (as you can see above).

# have a look and plot
print(Saatchi.vegC)
print(plotSpatial(Saatchi.vegC))

Now we formally compare the layers using the compareLayers() function, which also gives us some comparison metrics. The resulting Comparison object can be plotting with plotSpatialComparison().

``` {r comparison, fig.asp = 1, echo=TRUE}

compare layers to produce a Comparison object

vegC.comparison <- compareLayers(field1 = GUESS.vegC, field2 = Saatchi.vegC, layers1 = "Tree", layers2 = "vegC_std", override.quantity = TRUE)

have a look at this comparison object (lots of metadata tracked)

print(vegC.comparison )

plot difference map

print(plotSpatialComparison(vegC.comparison))

make a prettier version with country lines, an cleaner title and a legend title

print(plotSpatialComparison(vegC.comparison, text.multiplier = 1.5, title = "Tree C Mass Difference: LPJ-GUESS - Saatchi et al. 2011", map.overlay = "world"))

version with the reference publication with slightly larger text size

comparison_example_plot <- plotSpatialComparison(vegC.comparison, text.multiplier = 3.5, title = "Tree C Mass Difference: LPJ-GUESS - Obs.", subtitle = "Western Central Africa", map.overlay = "world")

also percentage difference map

print(plotSpatialComparison(vegC.comparison, type = "percentage.difference", limits = c(-300, 300), text.multiplier = 1.5, title = "Tree C Mass Difference: LPJ-GUESS - Saatchi et al. 2011", map.overlay = "world"))

plot side-by-side

print(plotSpatialComparison(vegC.comparison, type = "values"))

#### Selecting gridcells and plotting the subannual cycle

It is also possible to select gridcells using `selectGridcells()` or the `spatial.extent` argument to getField().  The `spatial.extent` can be a single gridcell (as a numeric two-vector of Lon and Lat), a data.frame or a data.table with columns Lon and Lat, or shapefile (as a SpatialPolygonDataFrame) and anything which can be coerced to a `raster::Extent` or a `terra::Ext`.  To select rectangular spatial extents, the `raster::crop()` method has been extended. Here we first read a Field with only one gridcell and then use the `plotSubannual()` function to plot the seasonal cycle.  Note that to get the monthly of LAI we are using the `"mlai"` quantity instead of the `"lai"` quantity.

``` {r select, fig.asp = 1, echo=TRUE}

# define a gridcell and the get a Field of the monthly LAI for the gridcell
gridcell <- data.frame(Lon = c(0.25), Lat = c(51.25))
London.mlai <- getField(source = GUESS.run, 
                        quant = "mlai",
                        spatial.extent = gridcell,
                        spatial.extent.id = "Approximately London")


# plot LAI not that by default the data are colour by year
print(plotSubannual(field = London.mlai))

# add the representation mean year 
print(plotSubannual(field = London.mlai, 
                    summary.function = mean,
                    summary.function.label = "Mean Year"))

# also avoid colouring by year by instead specifying to colour by the Source
print(plotSubannual(field = London.mlai, 
                    summary.function = mean,
                    summary.function.label = "Mean Year",
                    col.by = "Source"))


# Store example plot for reference publication
subannual_example_plot <- plotSubannual(field = London.mlai, 
                                        summary.function = mean,
                                        summary.function.label = "Mean Year",
                                        text.multiplier = 3.5)

Calculating biomes and categorical comparisons

The getScheme() function can be used to classify model output into some pre-defined classification categories. Here we classify vegetation into global biomes.

``` {r model biomes, echo=TRUE}

read and plot global biomes

biomes.model <- getScheme(GUESS.run, Smith2014BiomeScheme, year.aggregate.method = "mean")

print(plotSpatial(biomes.model))

There is also a global map on potential natural vegetation (PNV) biomes available in the package that can be compared to the model output.  First we define the data as a `Source` and read it.

``` {r data biomes, echo=TRUE}

# data biomes source
biomes.source <- defineSource(id = "DataBiomes",
                              dir = system.file("extdata", "NetCDF", "HandP_PNV", "HD", package = "DGVMTools"), 
                              format = NetCDF,
                              name = "Haxeltine and Prentice Biomes")
# note that because the dataset is already in the form of a categorical variable, applying classification is not necessary,
# therefore we can just call with getField() instead of getScheme() 
biomes.data <- getField(source = biomes.source, quant = "Smith2014")
print(plotSpatial(biomes.data))

Then we can compare. Note that compareLayers() function automatically detects categorical quantities (such as the biomes in this case) and treats them accordingly.

``` {r biomes comparison, echo=TRUE}

biome.comparison <- compareLayers(field1 = biomes.model, field2 = biomes.data, layers1 = "Smith2014", layers2 = "Smith2014")

have a look at this comparison object

print(biome.comparison )

plot side-by-side

print(plotSpatialComparison(biome.comparison, type = "values"))

plot difference map

print(plotSpatialComparison(biome.comparison, type = "difference"))

### Customising your plots (ggplot2)

Since the plotting functionality in DGVMTools is based entirely on ggplot2, the functions all return ggplot2 objects. This is handy because they can then be modified using the full power and flexibility of ggplot2.  Here are some examples, but remember the full array of ggplot2 functionality is at your disposal. If in doubt, seach the .

#### Changing facet layouts

You can provide the ggplot2-style arguments "nrow" and "ncol" directly to all the plot functions. Additionally,  `plotSpatial` adds a column called `Facet` by which to facet the data so you can also adjust the facetting after the fact with a ggplot2::facet_wrap() command. Examples follow:


``` {r facetting, echo=TRUE}
# make a plot with silly facets
lifeform.plot <- plotSpatial(LAI.year.mean, layers = c("Tree", "Grass", "Total"), ylim = c(35,45))
print(lifeform.plot)


# arrange facets with three rows using nrow in the plotSpatial() call
lifeform.plot <- plotSpatial(LAI.year.mean, layers = c("Tree", "Grass", "Total"), ylim = c(35,45), nrow =3)
print(lifeform.plot)

# arrange facets into two columns with a subsequent call to facet_wrap()
lifeform.plot <- lifeform.plot + facet_wrap(~Facet, ncol = 2)
print(lifeform.plot)

Formatting legends and text

The theme() and guides() functions are very handy. See ggplot2 documentation for more info, and here is an example:

``` {r legends, fig.height=6, echo=TRUE}

a make and view the base lot

biome.plot <- plotSpatialComparison(biome.comparison, type = "values") print(biome.plot)

move and format legend, and make the text slightly smaller

biome.plot <- biome.plot + theme(legend.position = "bottom", text = element_text(size = theme_get()$text$size * 0.75)) biome.plot <- biome.plot+ guides(fill = guide_legend(ncol = 4)) print(biome.plot)

## Further Details
Whilst the tutorial above will get you started with the package, the following will give you further information and a deeper understanding of the package and is therefore recommended reading.

### Objects/classes

Concepts in a typical DGVM workflow are encapsulated into objects (technically these are S4 classes).  For example, data always orginate from a particular source, typically a model simulation or a dataset.  This is represented as a `Source` object, which includes the location of the data on disk, an identifying string and other pertinent pieces of information.  Once a `Source` has been defined, data can be extracted from it which will be stored as a `Field` object.  Conceptually, a `Field` is the data of one particular variable (say biomass) over a period of time and an extent in space (say monthly over the whole globe for the years 1961 to 1990).  `Fields` are the most important objects in DGVMTools because they are manipulated by the users and they contain data (most other classes are metadata). One can think of a `Field` as being similar to raster Stacks or Bricks, but with more metadata and implicit temporal information.  Like Stacks/Bricks, they also have layers, for example one layer per PFT.  These layers can be operated upon, for example to find the total of all of the layers, or the fraction of one layer divided another, or many other operations.  They can also be plotted, and two `Field` objects can be compared to produce a `Comparison` object which has data on the difference between the two fields and statistical summary metrics.

Apart from the all-important `Field`and `Source` (and `Comparison`) objects, there are a few other objects defined in DGVMTools.  Each `Source` object includes `Format` object which describes the format of that data on disk - i.e. how to read data from it, which PFTs one might expect to find in a simulation etc.  Different physical quantities are described by `Quantity` objects which include units, a long name and a default colour scale for plotting. A `Layer` object contains metadata about a component of a `Field`, for example a particular plant functional type or C pool.  It can include metadata such as its growth form, phenology, leaf type or climate zone, but also a default colour for consistency across plots.  The `Period` class contains information about sub-annual time periods such as months or seasons.  Users can also define custom periods to define special growing seasons, or 8-day MODIS intervals, or whatever.  The `BiomeScheme` class contains data about different discrete biome classifications (including rules for how to derive biomes from model output).  Also, spatial-temporal-annual information is encapsulated in an `STAInfo` meta-data object but this mostly is used behind the scenes.

Many objects have an `id` which should be a unique alphanumeric character string, typically used for making filenames, and a `name` character string which can be formatted more prettily for use on graphics (although it should also be unique within one analysis).

### Functions/methods

Obviously there are quite a lot of the fuctions in the package as it aims to go from raw model output data to publication-quality plots.  Typically a script will start with calls to `defineSource()` to initialise one or more `Source` objects.  The next step might be to call `getField()` which actually extracts data from the `Source` and builds a `Field` objects.  The `getField()` function has options for selecting data in space and time, and also for aggregating it (for example to calculate the average of all years, or the total across a region).  However you can further refine the spatial-temporal dimensions of `Field` by applying the `aggregate` family of functions (`aggregateSpatial()`, `aggregateYears()` and `aggregateSubannual()`) or the `select` family of function (`selectGridcells()`, `selectDays()`,  `selectMonths()`, `selectYears()`, `selectSeasons()`) and the `crop()` method (which crops the spatial extent exactly like the `raster::crop` function, indeed, it even work with raster and terra objects).

Having sorted out the spatial and temporal dimemsions, you can then manipulation the `Field` and plot it.  Rember that data is stored as 'layers' similar to GIS and raster?  To combine layers and perform other operations one can use the `layerOp()` function which supports arithmetic operators (+,-,/,*), finding maximum and minimums of a group of layers, and the mean, standard deviation and variance across layers.  One can plot layers in spatial layers (i.e. maps) with `plotSpatial()` and temporal layers (time series) with `plotTemporal()`.  Two layers can be quickly compared with a scatter plot using `plotScatter()`.  For a more rigorous comparison, one can use `compareLayers()` to a produce a `Comparison` objects which can be plotted in a variety of ways with dedicated functions.

### `whichLayers()` and the `.Layer` (dot-layer) notation for the `layerOp()` function

One special feature of DGVMTools that deserves it own mention at this point is the `.Layer` (dot-layer) short-hand notation that can be used with `layerOp()`.  The allows the selecting of layers based on PFT metadata.  For example `.Tree` will select all layers corresponding to a tree PFT.  `.Evergreen` will select all layers with evergreen phenology, `.Boreal` will choose all PFTs associated with the boreal zone, etc.  This is done by a hidden call to the useful function `whichLayers()` which allows you to list all PFTs in a `Field`, or to list all PFTs based on a criteria (say 'Tree', 'Evergreen' or 'Boreal' as above.)



### Applicability

The package treats model output and observed data in an equivalent fashion and so is equally useful for model-to-data, model-to-model and data-to-data analyses.  A companion R package, DGVMData, processes a selection of data to different resolution and outputs the data as netCDF files which can be read easily with the `NetCDF` Format included in DGVMTools.

Much of the package could be utilised for any data which follows comes distributed over space (longitude, latitude), time (multiple years of data with subannual time resolution of daily, monthly or annual) and and different conceptual layers (for example different vegetation types or soil carbon pools).  However it has several additional features which are focussed on vegetation modelling, namely: the ability to read data from selected DGVMs in their native formats (currently LPJ-GUESS(-SPITFIRE), aDGVM(1) and aDGVM2 are supported); facilities for handling metadata about Plant Functional Types (PFTs); and the means to classify vegetation model output into biomes.  


### Miscellaneous details

* DGVMTools objects can be really easily converted into rasters or data.frames, so you can mix and match the functionality you want and integrate DGVMTools into existing scripts.
* For speed, the data is stored internally as `data.table` objects which are similar to `data.frames` but faster.
* `Field` objects can be stored on disk and automatically re-read by the `getField()` function, potentionally saving a lot of time (since often the bottle neck in analysis scripts is reading the raw data from disk).
* Plotting is done using the package ggplot2 which is especially useful as plots can be modified after they have been created (much more flexibly than base or lattice graphics).  This maintains relative simplicity in the DGVMTools plot function, but high levels of flexibility through subsequent modifications with ggplot2 commands.
* Unlike raster objects, `Field`s do not need to be on grids of equally spaced longitude and latitude, or even on grids at all!  This is handy when working with sites, or with grid from climate models which often have unevenly spaced latitudes.
* Regridding is not supported by DGVMTools.


### Example plots for publications

The following code lays out the example plots for the publication describing DGVMTools.

``` {r Publication plots, fig.height=8, echo = TRUE}


# add a small margins to the plots before combining

spatial_example_plot <- spatial_example_plot + theme(plot.margin = margin(t = 0.2,  r = 0.2, b = 1.2, l = 0.2, unit = "cm"))
temporal_example_plot <- temporal_example_plot + theme(plot.margin = margin(t = 0.2,  r = 0.2, b = 1.2, l = 0.2, unit = "cm"))
subannual_example_plot <- subannual_example_plot + theme(plot.margin = margin(t = 1.2,  r = 0.2, b = 0.5, l = 0.2, unit = "cm"))
comparison_example_plot <- comparison_example_plot + theme(plot.margin = margin(t = 1.2,  r = 0.2, b = 0.5, l = 0.2, unit = "cm"))


if(require(ggpubr, quietly = TRUE)) {


  combined_plot <- ggarrange(spatial_example_plot, 
                             temporal_example_plot, 
                             subannual_example_plot, 
                             comparison_example_plot, 
                             ncol = 2, nrow = 2, 
                             labels = list("(a)", "(b)", "(c)", "(d)"), 
                             font.label = list(size = 35),
                             hjust = -0.4, vjust = c(1.5, 1.5, 0.3, 0.3)) 


  if(require(Cairo, quietly = TRUE)){

    Cairo::Cairo(file = "~/DGVMTools_Example_Figure.png", type = "png", width = 2000, height = 2000)
    suppressWarnings(print(combined_plot))
    dev.off()
  }

}

MagicForrest/DGVMTools documentation built on Aug. 23, 2024, 8:05 a.m.