Compare LPJ-Guess outputs against measurement points from external dataset

[antonycastro]

Hi,

I want to compare my LPL-Guess outputs against Brienen et al 2015 dataset (attached; doi:10.1038/nature14283). Brienen dataset are points from measurement stations and I need to analyze the yearly differences between the LPJ-Guess outputs and these at the given points' locations. Any suggestions on the best way to proceed with DGVM tools??

All the best,

briennen_for_comparison.zip

[MagicForrest]

Hi Antony,

Did you already solve this yourself?

I was on holiday for a few days and hadn't yet decided how to proceed with it.

[MagicForrest]

Okay. I can offer some tips if needs be. Basically you want to:

1. Read your data and format it into a data.table with the structure as DGVMTools needs it (basically with columns "Lon Lat Year ". You probably also want to split it by site, with one data.table per site.

2. Define a Source object using a call to defineSource(). This basically only metadata, but you need this in place for step 3.

3. Make one Field object for each of the data.tables you made in step 1. These are the final objects that you can compare to the LPJ-GUESS data. You will need to do this using the 'new()' function. The 'source' and 'data' slots you will be filled with the output of steps 1 and 2. You will also need to to define an "STAInfo" object and id for each Field. This is a little bit tedious (under normal circumstances the package would do this automatically), but not too tricky.

But, yeah, ask if you have questions!

[antonycastro]

I have managed to get to point 3 where I make a field object for one of the data.tables. But don´t know how could I define the "STAInfo" and ID objects for each field. Could you detail a bit more this part, please?

[MagicForrest]

Sure. First define the STAInfo. This contain at the Spatial Temporal and Annual Information about a Field. Something like:

`

this.STAInfo <- new("STAInfo", first.year = the first year of data, last.year = the last year of data, year.aggregate.method = "none", spatial.extent = extent(the data.table), spatial.extent.id = the name of the site, spatial.aggregate.method = "none", subannual.resolution ="Annual", subannual.aggregate.method = "none", subannual.original = "Annual")

` (stuff in italics you need to specify, the rest is verbatim)

Then you also need to define a Quantity (this will the same for all your Fields):

`

AGBinc.Quantity <- new("Quantity", id = "AGBnetchange", name = "Above ground biomass net chance", units = "kg C /m^2", ## ??? colours = viridis::viridis, format = "DGVMData") `

Then create the id:

`

this.field.id <- makeFieldID(source = your source object you made in step 2, var.string = AGBinc.Quantity@id, sta.info = this.STAInfo)

`

And then build the Field:

`

this.field <- new("Field", id = this.field.id, data = the data.table from step 1 , quant = AGBinc.Quantity, source = the Source from step 2, this.STAInfo) `

What I write above might not be perfect, so just ask again if something doesn't work out ;-). But in principle, when it works, you can manipulate and plot the resulting Field exactly as you would do for LPJ-GUESS output.

[antonycastro]

Thank you very much! I have made progress.

Now I am getting an error when intending to compare a Brienen measurement station against the LPJ-Guess output:

vegC.comparison <- compareLayers(field1 = field, field2 = cmass.full, layers1 = "AGB_fin", layers2 = "Total") Error in merge.data.table(x = to.dt, y = layers.to.add.dt, all.y = keep.all.from, : A non-empty vector of column names for by is required.

I wonder what I am doing wrong. Please find attached all the files needed to reproduce what I have done.

Looking forward to your guidance, best regards

Magic_Forest_files.zip

[MagicForrest]

Hi Antony,

Good job on all that. Thanks for the script and files, it made troubleshooting very straightforward. And you did everything right! The problem was a small bug in the DGVMTools, which has now been fixed in the master branch so go ahead and re-install the package and it should work.

There is another hitch though, again on the package side, not your fault. The plotTemporal() function only works on Field objects, not Comparison objects. There should be a plotTemporalComparison() function, which I had completely forgotten to do! Sorry, I am not sure when I will get to that. Maybe soon, I'm not sure. In the mean time there are two things you can:

1. Visualise your data as a scatter plot using plotScatterComaprison().

2. Hack your Comparison object into a Field object, based on everything you have doen so far this should be too hard...

[MagicForrest]

Hi Antony,

I forgot something, you can use plotTemporal() toplot multiple time series (from different Fields) side-by-side. I just made some tweaks to the package so update from master branch, and then you can do this:

print(plotTemporal(list(field, cmass.full), layers = c("AGB_fin", "Total"), gridcells = c(-60, -2.5), facet = FALSE))

[you can also try facet = TRUE]

This will put the simulation and the observed on the same line graph. It needs some refinement, but it will get you visualising you results.

I will work on improving this (and maybe implementing plotTemporal() function) next week.

[antonycastro]

Everything worked. I managed to include Brienen´s points dataset into the DGVMTools format successfully and compare these against my model´s output.

Thank you very much, all the best Antony