Chapter 6 Advanced usage
Chapters 3, 4 & 5 illustrated a simple workflow, but IsoriX allows for much more complex workflows. The goal of this chapter is to illustrate such possibilities.
6.1 Working with plots
6.1.1 Fiddling with the plotting functions
If you want to change anything on the existing plot, have a look at ?plots.
The help file should provide you with all the details you need to customize your plots.
Here is are a couples of example:
plot(EuropeIsoscape,
sources = list(pch = 3, col = "orange"),
borders = list(col = "white"),
mask = list(fill = "darkgrey"),
palette = list(range = c(-130, 10), step = 1, n_labels = 10, fn = "rainbow"))
plot(EuropeIsoscape,
title = "H Isoscape",
sources = list(draw = FALSE),
borders = list(borders = NULL),
mask = list(fill = "black"),
palette = list(range = c(-130, 20), step = 30, fn = NULL)) 
You can see that it is possible to provide a function defining the colors using the argument fn from the list palette.
Note that when it is set to NULL, the famous palette viridis is used instead of our default palette.
If you need to change things on the plots that are not covered by the plotting functions from IsoriX, you will have to learn how to use lattice. The package lattice is very powerful but not always very intuitive to use, so we are going to list here a few tips.
The easiest way to modify a plot is to start saving the plot into an object:
Then, the (generic) function update() allows you to change many aspects of the plot whenever applied to your plot object (here isoscape_plot).
For example, you can increase the size of the text in the previous plot as follows:
As an entry point to see what can be changed in a lattice plot, you should consult the help page that shows up when typing ?lattice::xyplot in your R console.
Multipanel plots created after assignments can be modified using the same principle. For example, you can rename the title of each panel and modify the organisation of the panels as follows:
assignment_plot <- plot(AssignedBats,
who = 1:14,
sources = list(draw = FALSE),
calibs = list(draw = FALSE),
assigns = list(draw = FALSE)
)
update(assignment_plot,
par.settings = list(fontsize = list(text = 8)),
strip = lattice::strip.custom(factor.levels = paste("Bat ", 1:14)),
layout = c(col = 2, row = 7)
)
If you want to add things on top of the plots, you can do it as you would do it for any plot created using lattice. For example, consider you want to add the point that is the most compatible with the unknown origin of bats (see chapter 5). The first step is to recover the coordinates for such a point:
library(terra)
coord <- crds(AssignedBats2$group$pv)
MaxLocation <- coord[which.max(values(AssignedBats2$group$pv)), ]
Maximum <- data.frame(long = MaxLocation[1], lat = MaxLocation[2])
MaximumWe can then plot this information on top of the assignment plot by simply typing:
plot(AssignedBats2, who = "group", plot = FALSE) +
xyplot(Maximum$lat ~ Maximum$long,
pch = 13, col = "orange", cex = 5, lwd = 2, panel = panel.points)
6.1.2 Why you should save plots & how to export nice looking plots?
Displaying the plot directly using R is very inefficient when it comes to high resolution isoscapes. It may take a very long time and it may not fully work.
It is better to directly plot the figures into files.
Geeky note: Since the lattice system is not part of ggplot2, the function ggsave() should not be used.
6.1.2.1 Using base R graphics devices functions
A first option is to save your plot using one of the base R function such as png() or tiff().
For example, you could do:
png(filename = "output/Myisoscape.png",
width = 1920,
height = 1080,
res = 200)
plot(EuropeIsoscape)
dev.off()or
tiff(filename = "output/Myisoscape.tiff",
width = 1920,
height = 1080,
res = 200)
plot(EuropeIsoscape)
dev.off()As you can see, you first initialize the creation of the plot with the function png() or tiff(), then you call your plot, then you tell your computer that you are done with dev.off().
As arguments, you probably want to specify the dimensions and the resolution of the file.
The height and width are here considered to be in pixels (default setting, you can choose other units if you want using the argument units). Here we chose the so-called Full-HD or 1080p resolution (1080x1920) because we wanted to observe the isoscape carefully on a monitor of that resolution.
If your screen is Full-HD, try it and display the plot in full screen to get better results (if the plot does not match the definition of your screen, things can get ugly), if your screen is not, try another resolution.
If you don’t know what resolution your screen has, you can visit https://screenresolutiontest.com.
The parameter res is very useful as it rescales the line and fonts in the plot.
So if everything is too small just increase the value and if everything looks too bold and ugly, decreases it.
You also should specify the filename.
If like us, you don’t indicate the path in the filename, the file will be stored in your working directory which you can easily see by typing getwd().
So after running such code, simply go in your working directory and open the right file.
We prefer the PNGs files since they are much lighter than their TIFF counterpart, but one format may lead to a better output for you. The nice thing about these functions is that there do not require you to install anything extra. They also allow you to set the resolution, which is sometimes asked by journals. So if that works for you do stick to that.
Unfortunately, depending on your system, usual plotting functions provided with R may not always work well for high resolution rasters.
For example, using png() or tiff() may lead to display artifacts showing some random white lines.
Changing the value of the argument res is often enough to get rid of these lines.
Another more effective way is to use the package Cairo.
6.1.2.2 Using Cairo graphics devices functions
One set of plotting functions that seem to never suffer from white line artifacts are those provided by the package Cairo.
To use such functions, you first must make sure Cairo is installed and can be loaded:
If it complains, try to see why. If you have not installed it, just install the package. One noticeable source of trouble comes with MacOS. On this operating system, the package Cairo does not always install/load successfully. Problems happen when the program xquartz is not present on the system. This program used to be shipped by defaults on old version of MacOS but it is no longer the case. So if you are in trouble, please start by installing xquartz (outside R).
If Cairo is installed and loaded successfully, typing ?Cairo will show you all possible file formats you can save your plots into (PNG, JPG, TIFF, PDF…).
Here we will show how to save our isoscape both as a PNG and as a PDF.
Let us start by creating a PNG file with the main isoscape:
CairoPNG(filename = "output/Myisoscape.png",
height = 1080,
width = 1920,
res = 200)
plot(EuropeIsoscape)
dev.off()As you can see, once you have managed to install Cairo it becomes as simple as using the base R functions.
PDFs derived from raster objects tend to be heavy and to be badly rendered by most viewing software, so we recommend you to stick to PNGs. Yet, if you must, here is how to create PDFs using Cairo:
The arguments slightly differ with the function creating PNGs.
The argument setting the name of the file is now called file and not filename, the resolution has now to be provided in inches and the argument res does not exist for PDFs, but default results are usually fine.
Geeky note: You can also use the Cairo graphics device to improve the rendering of documents produced using knitr, which is precisely what we are doing to render the plots of this bookdown.
The trick is to simply add dev='CairoPNG' in the option of the chunks producing plots.
6.2 Exporting spatial objects to GIS
It is straightforward to export all spatial objects created by IsoriX into formats compatible with main software for Geographic Information System (GIS). This can be done using multiple R packages. Here is an example of how to export a GTiff raster using the package terra:
6.3 Weighted isoscapes
For several applications it is recommend to produce weighted isoscapes. We will consider here the case of precipitation-weighted annual average isoscapes, which is often recommended when isoscapes are derived from the isotopic composition measured in precipitation water. Yet, the instructions that follow could also be adapted to other type of isoscapes that must be weighted for other reasons.
The rational for using such precipitation-weighted isoscape is that the amount of precipitation influences how much specific isotopologues of water molecules end up in the soil, with repercussion on the entire food web. Here the weights we will use are the monthly cumulative amount of precipitation water averaged over a year, but you could also adapt the procedure to use a different temporal resolution.
You can follow two different approaches leading to weighted isoscapes. The two approaches differ according to when the weighting is being implemented in comparison to when the statistical models behind the isoscapes are being fitted:
the first strategy is to weigh the source data used to fit the isoscape. In this case, only the step 3.1 & 3.2 must be performed differently, but the rest of the entire workflow presented in chapters 3, 4 & 5 remain the same.
the second strategy is to build several non-weighted isoscapes (e.g. one for each month of the year) and then to combine and merge those isoscapes while applying the weights.
Note: It is not clear to us which of these two approaches is better (please let us know if you do!). To our knowledge the first approach is the one used by most, but this could be because it is easier to implement. We suspect the second approach to be better, but at the time of the writing, IsoriX does not yet allow to use that method in workflows involving assignments. So for now, your choice is simple: if you must perform an assignment stick to the first method; otherwise, you may try both and compare the obtained results.
6.3.1 Method 1: weighting before fitting
6.3.1.1 Preparing the source data
To create precipitation-weighted annual average isoscapes, the preparation of source data – the data containing the isotopic composition found in precipitation – must differ a little from what is described in section 3.2 since we need to retain the information about precipitation amounts.
rows_missing_or_unreliable_info <- is.na(rawGNIP$H2) |
is.na(rawGNIP$day.span) |
rawGNIP$day.span > -25 |
rawGNIP$day.span < -35 |
is.na(rawGNIP$Precipitation)
columns_to_keep <- c("Name.of.Site", "Latitude", "Longitude", "Altitude",
"Year", "Month", "H2", "Precipitation")
GNIPData_with_precip <- rawGNIP[!rows_missing_or_unreliable_info, columns_to_keep]
colnames(GNIPData_with_precip) <- c("source_ID", "lat", "long", "elev", "year", "month", "source_value", "precip")Which yields the following data:
We can check how many rows where discarded due to missing precipitation amount by comparing the object GNIPData created in section 3.2 to the object GNIPData_with_precip just created:
## [1] 62279## [1] 55803The row number is still high although 6476 were lost in the process.
6.3.1.2 Processing the source data with weighting
We can now process the source data we just prepared.
For now the function prepsources() does not handle the weighting.
We will add this feature soon, but we will prepare the source data by hand for the time being.
We will use the package dplyr for that since it makes aggregation much more straightforward than otherwise.
We start by aggregating and weighting the data over the different months per year:
library(dplyr) ## install this package beforehand if you don't have it
GNIPData_with_precip |>
filter(long > -30 & long < 60 & lat > 30 & lat < 70) |>
group_by(source_ID, year) |>
filter(n() > 6) |> ## only select years with more than 6 months of data
mutate(w = precip/sum(precip)) |>
summarize(total_precip = sum(precip),
mean_source_value = sum(w*source_value),
var_source_value = (n()/(n() - 1)) * sum(w*(source_value - mean(source_value))^2),
n_source_value = n(),
lat = unique(lat),
long = unique(long),
elev = unique(elev)) |>
as.data.frame() -> GNIPData_with_precipEUagg_yearlySo the only difference with what we did in section 3.2 is that we did not aggregate across the years and that we computed a weighted mean and variance of the isotopic measurements.
We can check the data we just created:
We then aggregate the data further, for each location, by computing averages over the different years. That way, we attempt to represent the situation for an average year. We can again use dplyr to do this, as follows:
GNIPData_with_precipEUagg_yearly |>
group_by(source_ID) |>
summarize(mean_source_value = mean(mean_source_value),
var_source_value = mean(var_source_value),
n_source_value = sum(n_source_value),
long = unique(long),
lat = unique(lat),
elev = unique(elev)) |>
as.data.frame() -> GNIPData_with_precipEUaggWe can check the data we just created:
6.3.1.3 Building the precipitation-weighted annual average isoscapes
We can now use these data and follow the rest of the workflow as described in more details in chapter 3, 4 & 5.
For example, you can create the precipitation-weighted annual average isoscapes exactly as described in chapter 3 with the only difference being the data provided to isofit():
EuropeFit_weighted <- isofit(data = GNIPData_with_precipEUagg,
mean_model_fix = list(elev = TRUE, lat_abs = TRUE))
EuropeIsoscape_weighted <- isoscape(raster = ElevEurope, isofit = EuropeFit_weighted)Let us now compare the mean isoscapes produced with and without weighting by precipitation amount:
The two isoscapes look very similar, but we can clarify what the differences are by creating a map showing the differences between the two mean isoscapes:
levelplot(EuropeIsoscape_weighted$isoscapes$mean - EuropeIsoscape$isoscapes$mean,
margin = FALSE,
main = "mean EuropeIsoscape_weighted - mean EuropeIsoscape") +
layer(lpolygon(CountryBorders, border = "white")) +
layer(lpolygon(OceanMask, col = "black"))
Let us similarly compare the isoscapes showing the prediction variance and the residual variance around the mean isoscape values:
levelplot(EuropeIsoscape_weighted$isoscapes$mean_predVar - EuropeIsoscape$isoscapes$mean_predVar,
margin = FALSE,
main = "predVar EuropeIsoscape_weighted - predVar EuropeIsoscape") +
layer(lpolygon(CountryBorders, border = "white")) +
layer(lpolygon(OceanMask, col = "black"))
levelplot(EuropeIsoscape_weighted$isoscapes$mean_residVar - EuropeIsoscape$isoscapes$mean_residVar,
margin = FALSE,
main = "residVar EuropeIsoscape_weighted - residVar EuropeIsoscape") +
layer(lpolygon(CountryBorders, border = "white")) +
layer(lpolygon(OceanMask, col = "black"))
As you can see the residual variance is quite different because it now represents the (weighted) variation in isotopic values between months for an average year and no longer the total variation combining both the variation across months and years.
You could then continue with the calibration and assignment steps exactly as described in the chapters 4 & 5.
For this you simply need to replace the object EuropeIsoscape used in these chapters by the object EuropeIsoscape_weighted just created.
6.3.2 Method 2: weighting after fitting
This second method takes a very different approach as the one just illustrated to produce weighted isoscapes.
It is still highly experiemental and should probably not been used yet since statistical details need to be figured out (computation of variances may be wrong).
As mentioned earlier, this method also does not (for now) allows for workflow involving assignments, but we will attempt to do this in the future.
The principle of this second weighting method remains simple: instead of feeding weighted isotopic values to isofit(), the idea here is to actually build different isoscapes (one per month in the context of precipitation-weighted annual average isoscapes) and to merge them in a fashion accounting for weights.
6.3.2.1 Preparing & processing the source data
Here the preparation and processing of the source data can be done at once with the use of the function prepsources() using the argument split_by = "month":
GNIPDataEU12agg <- prepsources(data = GNIPData,
long_min = -30, long_max = 60,
lat_min = 30, lat_max = 70,
split_by = "month")Let us display the data:
As you can see, we now have observation for each month.
The column mean_source_value gives the average for each month across the years.
The column var_source_value gives the between-year variance in monthly measurements.
Contrary to what we did in the first weighting method, we do not need to manipulate precipitation amounts at this stage.
6.3.3 Downloading the precipitation data
For this second weighting method we need more precipitation data than those included in the GNIP raw data.
Indeed, we actually need predictions for the precipitation amount all over the landscape covered by the isoscapes.
One source of such data is https://worldclim.org/data/worldclim21.html which provides datasets at different spatial resolution.
You can download and unzip manually the file you want, but IsoriX provides the function getprecip() which directly downloads and unzips the file presenting the highest resolution (30s or ~ 1 km\(^2\)).
You can use getprecip() as follows:
old_opt <- options(timeout = 500) ## modify R options
getprecip(path = "input/") ## for now do not forget the final '/' (we will fix that)
options(old_opt) ## restore original R optionsGeeky note: You have to modify the R options because, by default, R aborts any download lasting more than 1 min. In the future, we will try to do this behind the scene so you don’t have to worry about it.
6.3.4 Preparing the precipitation data
Once the precipitation data downloaded, you can now import them into R and reshape them to the right format and resolution.
This is what the function prepcipitate() is there for:
## NULLWarning: As you can see, you need to have prepared a structural raster (as shown in section 3.6.1) to run prepcipitate().
This implies that there is at the moment some circularity in the workflow we proposed since in section 3.6.1 we used a set of fitted models stored in the object EuropeFit which you may not have at this stage.
To circumvent this circularity, either perform first a simple workflow (i.e. produce a non-weighted isoscape) before producing a weighted isoscape, or use the argument manual_crop instead of isofit when calling prepraster().
Indeed, the object EuropeFit used in section 3.6.1 when calling prepraster() with its argument isofit = EuropeFit is only used to define the coordinates of the boundaries of the structural raster and that can be done manually instead using e.g. manual_crop = c(-27.34, 57.1, 30.08, 68.96).
We will try to smooth this out in a future version of IsoriX.
Let us examine the object PrecipitationBrickEU just prepared:
## class : SpatRaster
## dimensions : 623, 1352, 12 (nrow, ncol, nlyr)
## resolution : 0.06873321, 0.06873649 (x, y)
## extent : -31.55346, 61.37383, 28.08577, 70.9086 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : month_1, month_2, month_3, month_4, month_5, month_6, ...
## min values : 0.0000, NaN, NaN, NaN, NaN, NaN, ...
## max values : 396.2898, NaN, NaN, NaN, NaN, NaN, ...This set of rasters (combined in an object of class RasterBrick) contains average precipitation amounts for each month and for each location defined by the structural raster.
We can for instance visualise the distribution of such precipitation amounts for January and June using the function levelplot() (plot() also works but the output is different):
## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf## Error in colorkey$space: $ operator is invalid for atomic vectors6.3.5 Fitting the models for a multiple set of isoscapes
We now use the function isomultifit() for fitting one set of isoscapes per month.
Remember that GNIPDataEU12agg indeed contains data at the monthly level (averged across years).
As for isofit(), several parameters can be adjusted to fit different models in isomultifit(), and we will use here the settings seen previously:
EuropeFit12 <- isomultifit(data = GNIPDataEU12agg,
mean_model_fix = list(elev = TRUE, lat_abs = TRUE))Geeky note: For the time being, the function isomultifit() fits sequentially the sets of isoscapes.
So for the example used here, it means that 12 sets of isoscapes are fitted sequentially, which thus takes more computing time than building a single set of isoscape.
In the future our plan is to revise the internal structure of this function for it to be able to use several CPUs simultaneously and thus benefit from the IsoriX option used to define the number of CPUs to work with (see section 2.5).
6.3.6 Building the precipitation-weighted annual average isoscapes
We can now use the set of isoscape models to build the isoscapes using the function isomultiscape() to which we pass the raster brick storing the precipitation amounts as weights:
EuropeIsoscape_weighted2 <- isomultiscape(raster = ElevEurope,
isofit = EuropeFit12,
weighting = PrecipitationBrickEU)## ### Multilayered raster containing the isoscapes
## class : SpatRaster
## dimensions : 623, 1352, 8 (nrow, ncol, nlyr)
## resolution : 0.06873321, 0.06873649 (x, y)
## extent : -31.55346, 61.37383, 28.08577, 70.9086 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source(s) : memory
## names : mean, mean_predVar, mean_~idVar, mean_respVar, disp, disp_predVar, ...
## min values : NaN, NaN, NaN, NaN, NaN, NaN, ...
## max values : NaN, NaN, NaN, NaN, NaN, NaN, ...
##
## ### Points of the sampled sources
## class : SpatVector
## geometry : points
## dimensions : 326, 0 (geometries, attributes)
## extent : -27.34, 57.1, 30.08, 68.96 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## NULL
## NULLYou can now plot this isoscape:
## Warning in max(var): no non-missing arguments to max; returning -Inf## Warning in min(var): no non-missing arguments to min; returning Inf## Error in seq.default(min(range, na.rm = TRUE), (max(range, na.rm = TRUE)), : 'from' must be a finite numberWarning: In this second weighting method, the meaning of the residual variance is difference from that of non-weighted isoscapes and from that of isoscapes weighted using the first method. We will revisit this in a future version of IsoriX.