MagicForrest / DGVMTools

R package for processing, analysing and visualising ouput from Dynamic Global Vegetation Models (DGVMs)
GNU General Public License v3.0
29 stars 22 forks source link

Calling addArea() with non-evenly spaced gridlist/coordinates #91

Open hol430 opened 4 months ago

hol430 commented 4 months ago

When attempting to use the addArea() function with some LPJ-Guess outputs, I get the following error:

Error in DGVMTools::addArea(agg, lon_centres = lons, lat_centres = lats) : 
  In addArea(), not all latitudes in the data are present in the supplied in the 'latitude_centres' argument.

My dataset here is a gridded lpj-guess output of LAI over Australia, and the coordinates I passed into the function were simply all unique coordinates in my gridlist. After some investigation, I think the problem comes down to the fact that my gridlist doesn't contain an evenly-spaced sequence of latitudes (note the missing -40.5° and -39.5° - these are ocean gridcells at all longitudes if your grid is a small box around Australia):

> lats
 [1] -43.5 -42.5 -41.5 -38.5 -37.5 -36.5 -35.5 -34.5 -33.5 -32.5 -31.5 -30.5
[13] -29.5 -28.5 -27.5 -26.5 -25.5 -24.5 -23.5 -22.5 -21.5 -20.5 -19.5 -18.5
[25] -17.5 -16.5 -15.5 -14.5 -13.5 -12.5 -11.5

I think the problem is caused by the use of extract.seq() in addArea(). In this gridlist, less than half of the latitudes are unevenly spaced, so extract.seq() returns a full sequence. Is there a reason to do this? If I'm reading things correctly, it seems that returning a full sequence for an unevenly-spaced grid will always result in an error in addArea(). That said, I'm not certain what the function should output in my case - should it account for the presence of oceans, or simply assume that there are some large gridcells in that area?

I tried changing the addArea() function to use sort(unique()) instead of extract.seq(), which appears to work on some level (I no longer get the error), though I'm not sure if the output should be considered "correct" around the affected area (the bass strait in my case).

If it helps, here is a reproducible example which demonstrates the problem:

mlai.out.gz

library("DGVMTools")

infile <- "/path/to/mlai.out.gz"

src <- defineSource("src", "LAI Source", dirname(infile), format = GUESS)
field <- getField(src, quant = "mlai")
agg <- aggregateYears(aggregateSubannual(field, method = "max"), "mean")

lons <- unique(sort(field@data$Lon))
lats <- unique(sort(field@data$Lat))

d <- addArea(agg, lon_centres = lons, lat_centres = lats)