Spatial output analysis
< Regional output analysis | Overview | Creating a new realization >
Level: ⚫⚫⚫⚪⚪ Advanced
Requirements
Local copy of the MAgPIE model (https://github.com/magpiemodel/magpie)
Have R installed (https://www.r-project.org/)
Have R packages
madrat,magclass,luplotandterrainstalledContent
Know where to find spatially explicit output data (NetCDF & magclass-format)
Data visualisation with the Panoply NetCDF Data Viewer
Learn how to make animated time series with output data
Use the
luplotpackage for spatial visualisation of magclass objects in RTransform
magclass-object to aSpatRaster-object of theterrapackageOverview
Introduction
In the previous tutorial, we have learned how to analyse MAgPIE results at the regional scale. The next step, covered by this tutorial, is the analysis of spatially explicit output data at the 0.5 degree level. We therefore build on knowledge covered by the previous tutorial on how to access the output folder of your local MAgPIE version. If you are unsure, how to access the output folder, or if you have not yet conducted a MAgPIE simulation, please, have a look at the previous tutorial(s).
MAgPIE itself does not operate at the 0.5 degree level, yet many procesess like land competition, agricultural production or irrigation are based on input data at the 0.5 degree grid cell level, which is aggregated to spatial simulation units (clusters) in order to address computational constraints.
Disaggregation scripts
After optimisation, some of the model outputs can be disaggregated back to the 0.5 degree level via dedicated disaggregation scripts that can be found under "<your MAgPIE folder>/scripts/output/extra". In the default MAgPIE configuration only the script "extra/disaggregation" is selected via cfg$output:
cfg$output <- c("output_check", "extra/disaggregation", "rds_report")
Other disaggregation scripts can be added by changing cfg$output before conducting a model run, e.g. by adding "extra/disaggregation_LUH2".
cfg$output <- c("output_check", "extra/disaggregation", "rds_report", "extra/disaggregation_LUH2")
Visualisation of spatially explicit MAgPIE output with the Panoply NetCDF Data Viewer
Spatially explicit MAgPIE output is usually stored in two different data formats in the output folder, including the MAgPIE-own magclass-format with the ending ‘.mz’ and the NetCDF format with the ending ‘.nc’. This section covers the handling and visualisation of spatial ouputs in the NetCDF format by using the NetCDF data visualisation tool Panoply. The magclass-format can be read-in and processed in R by using the R-packages magclass or luplot for visualisation. We will come back to this, however, at a later stage.
Installling Panoply
In order to download Panoply go to the website https://www.giss.nasa.gov/tools/panoply/. Scroll down until you find the header Get Panoply and click on Download Panoply. This will lead you to the following page:
Please, choose the right version for your operation system, select a destination folder and proceed with the download. After the download has finished, follow the instructions for installing Panoply that can be found on the website for your operating system.
Make sure Java 11 (or later) is installed
If you encounter an error while trying to execute Panoply after installation, make sure you have Java 11 (or later) available on your computer. In case you need to install/update Java, go to the website https://www.oracle.com/uk/java/technologies/downloads/#jdk21-windows and download the appropriate version for your operating system (e.g. JDK 21).
After the installing Java, please make sure that the JAVA_HOME is set correctly. To set the JAVA_HOME variable, locate your Java installation directory. If you chose JDK 21, it will be something like C:\Program Files\Java\jdk-21. On Windows 10/11 you can search for Edit the system environment variables. Select Environment Variables… and under System variables see whether the JAVA_HOME variable is set correctly. If it’s not yet set, press New… and enter JAVA_HOME in Variable name: and under Variable value: enter the path to your Java installation, here C:\Program Files\Java\jdk-21. Press OK to set the JAVA_HOME variable.
Using Panoply to visualise MAgPIE output
Navigate to the folder, in which Panoply is stored and click on Panoply.exe to execute the software. Each time a new session is started in Panoply, a new window opens that allows you to select your data. Navigate to the output folder of your local MAgPIE version and select a subfolder containing the results of a succesful model run. Panoply will display all NetCDF files that are available in this folder.
Select the file ‘cell.land_0.5.nc’ and press open. ‘cell.land_0.5.nc’ includes a spatially explicit time-series of the land area (Mha) covered by the different land-types (cropland, pasture, etc.) per 0.5 degree grid cell. We now see a window that lists all information that is stored in this file, including the names and the type of data that are available. We can also look at the panel on the right-hand side to see some meta data, e.g. the units.
We can now select the item ‘crop’ by left-clicking it. Then click Create Plot, accept all defaults and press Create. This will create a plot similar to this example:
The plot shows the cropland area in each grid cell given in Mha. By default, Panoply will plot the data in a smoothed form. Uncheck interpolate in the Plot Controls Window to view the actual cell values.
Panoply also allows us to perform basic calculations and plot them. Now, let us compute the difference in cropland area per grid cell between the initial time step and 2050.
We go back to our data window and click on ‘crop’ once again.
But instead of pressing Create plot, this time we click on Combine plot. We now find two array tabs on the top left-hand side. Also, Panoply will prompt the Arrays-window for us. Under Plot select ‘Array 2 - Array 1’ and select the year 2050 in Array 2: crop, as shown below, and hit ‘Enter’.
We can also change the appearance of the plot by reprojection the map to a different coordinate reference system (CRS), or by changing the color scale.
In order to change the CRS, we can select Map Projection in the plot controls
and then, for example, change to the ‘Equal Earth’ projection by clicking on the Projection drop-down menu in the Map projection window.
After changing the CRS, the map looks like this:
Now let us improve the looks of it a bit and work on the color scale. Therefore, select the option Scale under Window.
This will open the following window, where we can tweak a range of map properties.
We can set the range of the legend to [-0.25; 0.25] and select the color bar ‘CB_PiYG_09.cpt’ or any color bar to your liking. However, consider that this is a change map and using a neutral color like grey or white in the center of your color range is recommended to facilitate the interpretation of the data. Also set the label format to ‘%.2f’. This will change the legend label format and provides us with two decimal places behind the comma. We can also add a custom caption by selecting Custom and enter a suitable text. Under Window -> Label we can also change/remove global plot labels. In the end out plot looks like this:
Creating an animated time series based on MAgPIE projections using Panoply
Panoply can also be used to create an animated time series (e.g. in MP4-format) based on spatially explicit MAgPIE projections. In the following, we will create an animated time series of changes in cropland area per grid cell, based on the map that we have created in the previous steps.
To create an animation go to Export Animation.
This will open up the following window:
In order to display relative changes, set Array 1: crop to ‘None (leave fixed)’ and select a frame rate of e.g. 1 fps, then click on Okay. Choose a file destination, choose a file name and press Save. Panoply will now create an MP4-animation of the complete time series.
Analysis of spatial outputs in R using the magclass and luplot libraries
For further analysis and processing of spatial MAgPIE projections beyond visual interpretation, spatial outputs can be imported and analysed in R using the packages magclass (https://github.com/pik-piam/magclass) and luplot (https://github.com/pik-piam/luplot). magclass provides a range of tools for handling MAgPIE’s own magclass-format (ending ‘.mz’) for storing spatio-temporal data. The magclass package also includes a vignette that describes the functionality of the package, including examples on how to handle spatio-temporal MAgPIE data.
In addition, the magclass package also provides the functionality for transforming the maglcass-format to the SpatRaster-format of the terra package, which offers a wide range of functions for further processing spatial data.
We can make ourselves familiar with the magclass package by opening an R (RStudio) session. Find the output folder of your local MAgPIE version and set the working directory to the subfolder containing the results of a finished model run.
setwd("/path/to/your/magpie/output/subfolder")
We can call the library and look at the help pages:
library(magclass)
?magclass
In order to get an overview of the functionality, we can click on the index and look for helpful functions, e.g.
read.magpie and navigate to the help/documentation page, by clicking on it.
We can now try to read in ‘cell.land_0.5.mz’ by typing
x <- read.magpie('cell.land_0.5.mz')
Then we can use str(x) to find out how the data is structured in x. We can analyse x and perform a range of calculations. Data dimensions in magclass objects are indexed in the following three-dimensional form
x['<cell>', '<time>', '<data>']
For example, we can compute the difference in the primary forest (‘primforest’) area in all grid cells between 2015 and 2050:
primforestDiff <- x[, 2050, 'primforest'] - x[, 2015, 'primforest']
We can also compute the difference for all land classes simultaneously by typing
allDiff <- x[, 2050,] - x[, 2015,]
Now, let us load the library luplot and look at the help pages:
library(luplot)
?luplot
The library contains helpful functions for plotting magclass objects in different ways, e.g. plotmap2 that transfers data stored in a magclass object to a simple map. However, we cannot plot all data at the same time. Therefore, we need to select a time step and a variable that we want to look at. Here, we choose the year y2050 and the variable ‘cropland’ again. In plotmap2 can also define the upper and lower boundary of the legend via legend_range.
luplot::plotmap2(x[, 2050, 'cropland'], legend_range = c(0, 0.3))
The data can also be transformed to a SpatRaster-object of the terra package.
To load terra and show the help pages use
library(terra)
?terra
In order to transform a magclass object to a SpatRaster-object the magclass package provides the function as.SpatRaster.
r <- as.SpatRaster(x)
> print(r)
class : SpatRaster
dimensions : 280, 720, 133 (nrow, ncol, nlyr)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -56, 84 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
source(s) : memory
names : y1985..crop, y1995..crop, y2000..crop, y2005..crop, y2010..crop, y2015..crop, ...
min values : 0.00000, 0.0000000, 0.0000000, 0.0000000, 0.000000, 0.0000000, ...
max values : 0.24699, 0.2469919, 0.2469941, 0.2470168, 0.247136, 0.2472784, ...
time (years): 1985 to 2100
The function also merges the names in the time and data dimensions in the form "<year>..<data>" in order to create unique layer names within the SpatRaster-object.
> names(r)
[1] "y1985..crop" "y1995..crop" "y2000..crop"
[4] "y2005..crop" "y2010..crop" "y2015..crop"
[7] "y2020..crop" "y2025..crop" "y2030..crop"
[10] "y2035..crop" "y2040..crop" "y2045..crop"
[13] "y2050..crop" "y2055..crop" "y2060..crop"
[16] "y2070..crop" "y2080..crop" "y2090..crop"
[19] "y2100..crop" "y1985..past" "y1995..past"
[22] "y2000..past" "y2005..past" "y2010..past"
...
The terra also provides some functionality to plot raster data. In order to plot the cropland area in 2050 use
plot(r[["y2050..crop"]])
Exercises (click on the arrows to uncover the solution)
The file `cell.bii_0.5.mz` in the MAgPIE output folder reports the Biodiversity Intactness Index (BII) in each grid cell. Read it into R and assign it to an object `bii`.
bii <- read.magpie('<path-to-output-subfolder>/cell.bii_0.5.mz')How many 0.5 degree grid cells and time steps does the object `bii` contain? Also use `magclass::getComment()` to find out in which unit the data is stored.
Find out how many 0.5 degree grid cells and time steps
biicontains by usingstr(bii). You will find out thatbiihas 67420 grid cells, 18 time steps and only one data layer.magclass::getComment(x)gives'unitless'Calculate the difference in BII per grid cell between 2015 and 2050 and store the difference in an object called `biiDiff`.
biiDiff <- bii[, 2050,] - bii[, 2015, ]Plot the difference in BII per grid cell between 2015 and 2050 using the `terra::plot()` function.
terra::plot(as.SpatRaster(biiDiff), range = c(-0.25, 0.25))
< Regional output analysis | Overview | Creating a new realization >