Merge pull request #25 from MagicForrest/dev-v0.8.0

Dev v0.8.0

Author: MagicForrest

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: DGVMTools
 Type: Package
-Version: 0.7.1
+Version: 0.8.0
 Title: DGVM Processing, Analysis and Plotting Tools
 Date: 2019-02-27
 Authors@R: c(person("Matthew", "Forrest", role=c("aut", "cre"), email="[email protected]"),
@@ -35,7 +35,6 @@ Suggests:
     rmarkdown,
     ncdf4,
     rgeos,
-    udunits2,
     testthat
 License: GPL | file LICENSE
 URL: https://www.bik-f.de
@@ -65,7 +64,6 @@ Collate:
     'commonSTAInfo.R'
     'compareLayers.R'
     'copyLayers.R'
-    'countCategoricalData.R'
     'crop-methods.R'
     'defineSource.R'
     'determinePFTs.R'
@@ -89,6 +87,7 @@ Collate:
     'plotSeasonal.R'
     'plotSpatial.R'
     'plotSpatialComparison.R'
+    'plotSubannual.R'
     'plotTemporal.R'
     'plotting-framework-functions.R'
     'renameLayers.R'

NAMESPACE

@@ -31,15 +31,13 @@ export(all.seasons)
 export(annual)
 export(availableQuantities)
 export(averageFields)
-export(byIDfromList)
 export(calcNewField)
 export(categoricalComparison)
 export(cmass.palette)
 export(commonSTAInfo)
 export(compareLayers)
 export(continuousComparison)
 export(copyLayers)
-export(countCategoricalData)
 export(defineSource)
 export(determinePFTs)
 export(difference.palette)
@@ -60,7 +58,6 @@ export(is.STAInfo)
 export(is.Source)
 export(is.Statistics)
 export(is.leapyear)
-export(lai.palette)
 export(layerOp)
 export(listPFTs)
 export(lookupQuantity)
@@ -75,6 +72,7 @@ export(plotScatterComparison)
 export(plotSeasonal)
 export(plotSpatial)
 export(plotSpatialComparison)
+export(plotSubannual)
 export(plotTemporal)
 export(promoteToRaster)
 export(proportionsComparison)
@@ -85,6 +83,7 @@ export(reversed.magma)
 export(reversed.plasma)
 export(reversed.tim.colors)
 export(reversed.viridis)
+export(seasonalComparison)
 export(selectDays)
 export(selectGridcells)
 export(selectLayers)

NEWS.md

@@ -0,0 +1,95 @@
+#  DGVMTools v0.8.0 (2019-03-18) 
+
+This release has a couple of new features, few minor bug fixes and some tidying and cleaning.  See also new repository for example scripts at https://github.com/MagicForrest/DGVMTools_Scripts
+
+## Feature improvements:
+
+* The getBiomes() function can now take a list of Source objects and average over them before calculating the biomes (#19)
+* Seasonal comparison metrics, including Mean Phase Difference, now included.  plotSpatial also works for seasonal comparisons.
+* The compareLayers() function can now calculate user-defined metrics.
+* Function plotSubannual() included as a replacement for plotSeasonal(). plotSeasonal() had an inconsistent interface with the other plot functions and is now deprecated.  In contrast, plotSubannual() has consistent arguments and also re-uses more code internally, however **facetting is not yet fully implemented**.
+* In selectGridcells(), the @spatial.extent now stores the gridcells or shapefile that was used.
+* Comparison objects now has an @type characters slot to say if they are "continuous", "categorical", "seasonal" or "relative.abundance" (this is also a potential break).
+
+## Potential breaks:
+
+* Fields are now saved to disk using the .RData file extension instead of .DGVMField to be more consistent with the R universe. You can simply rename your files of necessary.  
+* Ugly LAI colour palette removed.
+* Function countCategoricalData() removed.
+* Function byIDfromList() now internal.
+* Function addArea() makes layers called "Area" instead of "area".
+* The @stats slot of a Comparison object is now a simple R list instead of a class, giving much more flexibility and the possibility for user-defined metrics.
+
+## Bug fixes:
+
+* Fields ids now done more consistently.
+* Function renameLayers() now actually returns the Field
+
+## Behind the scenes:
+
+* More code from plotSpatial() factored out into plotting-framework-functions.R
+* Benchmarking code re-organised
+* getBiomes() much slimmed down by removing code for cropping/aggregating (this was redundant since it was done in the getField calls)
+* udunits2 dependency removed
+
+
+#  DGVMTools v0.7.1 (2019-02-27) 
+
+Bug fix due to a small regression in listPFTs().  Also took the opportunity to include examples for listPFTs, unit tests for listPFTs, and changed the 'print' functions to 'show' functions as is appropriate for S4 methods (but this does not effect users).
+
+
+#  DGVMTools v0.7.0 (2019-02-27)
+
+Changes are mostly bug fixes and behind the scenes improvements. However there have been a couple of user-facing changes which could break stuff (hence the minor version number increment from 0.6 to 0.7):
+
+
+## Potential breaks:
+
+* Statistics class has been removed in favour of a simpler and more flexible list
+* determineQuantities() and determinePFTs() are not availableQuantities() and availablePFTs()
+* The 'sources' argument in plotSpatial() is not more accurately called 'fields'
+
+## Feature improvements:
+
+* Shapefiles can be now used to select gridcells by using them as an an argument to selectGridcells(), this can also be done in getField() which allows lower memory footprint when selecting regions from large files.    
+* aDGVM support improved (thanks Simon!)
+* New benchmarking metrics included, also enhanced flexibility allowing custom metrics
+* plotSpatial() can call geom_tile() instead of geom_raster() which allows reprojections (ie polar views)
+
+## Bug fixes:
+
+* Many! (see commits, but none of them resulted in wrong results, but rather failing code)
+    
+## Behind the scenes:     
+
+* Unit tests!
+* Plotting code refactorised for consistency and efficient re-use within the package.
+* Benchmarking code re-worked
+* Day/Month/Year are now generally stored as integers (not numerics) to save a little memory
+    
+Also, from now on the master branch will be consider stable and usable.  Bug fixes and small feature improvements will be pulled directly into master.  For developers, there will be a dev-v0.8 branch to collect larger features and potentially breaking changes.
+
+
+# DGVMTools v0.6.2 (2018-10-29)
+
+Removed dependency on external gzip on windows machines for better portability.
+
+
+# DGVMTools v0.6.1 (2018-10-26)
+
+Featuring aggregation and time series plotting of daily data, use of gzip for compressed DGVMData (also on windows systems), extended vignette and bug fixes.
+
+
+# DGVMTools v0.6.0 (2018-10-23)
+
+For FireMIP workshop.
+
+
+# DGVMTools v0.5.2 (2018-10-12)
+
+For visiting collaborators.  Release v0.6 still planned in the next two weeks.
+
+
+# DGVMTools v0.5.0 (2018-10-01) 
+
+DGVMTools is now publicly available!  Whilst package is pretty much fully featured, there may be bugs, so it should be considered in the late-beta phase.  New v0.6 release to come this month.

R/Format-DGVMData.R

@@ -297,9 +297,11 @@ getField_DGVMData <- function(source,
 
   # simple ones
   spatial.extent.id  <- getGlobalAttribute("spatial.extent.id", global.attributes)
+  if(!is.null(spatial.extent.id)) { sta.info@spatial.extent.id <- spatial.extent.id  }
+  else { sta.info@spatial.extent.id <- "Full"  }
   spatial.aggregate.method  <- getGlobalAttribute("spatial.aggregate.method", global.attributes)
-  if(!is.null(spatial.extent.id)) sta.info@spatial.extent.id <- spatial.extent.id
   if(!is.null(spatial.aggregate.method)) sta.info@spatial.aggregate.method <- spatial.aggregate.method
+  else { sta.info@spatial.aggregate.method <- "none"  }
 
   # first attempt to use the attributes from the NetCDF file
   xmin  <- getGlobalAttribute("xmin", global.attributes)

R/Format-GUESS.R

@@ -167,13 +167,17 @@ openLPJOutputFile <- function(run,
     if(spatial.extent.class == "SpatialPolygonsDataFrame" || spatial.extent.class == "numeric" || is.data.frame(target.sta@spatial.extent) || is.data.table(target.sta@spatial.extent)) {
       dt <- selectGridcells(x = dt, gridcells = target.sta@spatial.extent, spatial.extent.id = target.sta@spatial.extent.id, ...)
       new.extent <- target.sta@spatial.extent
+      # if new.extent is a data.frame, convery it to a data.table for consistency
+      if(is.data.frame(new.extent) & !is.data.table(new.extent)) new.extent <- as.data.table(new.extent)
     }
 
     else {
       dt <- crop(x = dt, y = target.sta@spatial.extent, spatial.extent.id = target.sta@spatial.extent.id)
       new.extent <- extent(dt)
     } 
 
+   
+    
   }
 
   gc()
@@ -233,10 +237,16 @@ openLPJOutputFile <- function(run,
                  subannual.resolution = subannual,
                  subannual.original = subannual)
 
+  # if cropping has been done, set the new spatial.extent and spatial.extent.id
   if(!is.null(new.extent))  {
-    sta.info@spatial.extent = extent(dt)
+    sta.info@spatial.extent = new.extent
     sta.info@spatial.extent.id <- target.sta@spatial.extent.id
   }
+  # otherwise set
+  else {
+    sta.info@spatial.extent = extent(dt)
+    sta.info@spatial.extent.id <- "Full"
+  }
 
   gc()
 

R/Format-aDGVM.R

@@ -293,6 +293,13 @@ getQuantity_aDGVM_Scheme1 <- function(run, variable, target.sta, adgvm.daily)
   # Now that we have the data we can set a spatial.extent
   actual.sta.info@spatial.extent <- extent(out.all)
 
+  # Set the spatial extent to "Full" if one not provided
+  if(length(target.sta@spatial.extent.id) == 0) {
+    actual.sta.info@spatial.extent.id <- "Full"
+  }
+  else {
+    actual.sta.info@spatial.extent.id <- target.sta@spatial.extent.id
+  }
 
   # make the ID and then make and return Field
   field.id <- makeFieldID(source = run, var.string = variable@id, sta.info = actual.sta.info)
@@ -663,6 +670,14 @@ getQuantity_aDGVM_Scheme2 <- function(run, variable, target.sta)
   # Now that we have the data we can set a spatial.extent
   actual.sta.info@spatial.extent <- extent(out.all)
 
+  # Set the spatial extent to "Full" if one not provided
+  if(length(target.sta@spatial.extent.id) == 0) {
+    actual.sta.info@spatial.extent.id <- "Full"
+  }
+  else {
+    actual.sta.info@spatial.extent.id <- target.sta@spatial.extent.id
+  }
+  
   # make the ID and then make and return Field
   field.id <- makeFieldID(source = run, var.string = variable@id, sta.info = actual.sta.info)
 

R/addArea.R

@@ -78,7 +78,7 @@ gridarea2d <- function(lon, lat, scale=1.0, ellipse=FALSE) {
   area2d <- array(rep(area1d, each=length(lon)), c(length(lon), length(lat)))
   area   <- data.frame(Lon=as.vector(lon2d),
                        Lat=as.vector(lat2d),
-                       area=as.vector(area2d))
+                       Area=as.vector(area2d))
   area <- data.table(area, key=c("Lon", "Lat"))
   return(area)
 }
@@ -115,21 +115,62 @@ extract.seq <- function(x, force.regular=FALSE, descending=FALSE) {
 
 #' Adds the gridcell area to a spatial Field or data.table/data.frame
 #' 
-#' Adds the gridcell area to a spatial Field or data.table/data.frame.
-#' For unit conversion it makes use of \code{\link[udunits2]{ud.convert}}, if installed. If udunits2 is not installed the option is ignored and the returned unit is m^2.
+#' Adds the gridcell area to a spatial Field or data.table/data.frame.  Makes a new layer/column called "area".  
+#' Unit conversion from "m^2" (default) to "km^2" and "ha" is supported. 
 #' 
 #' @param input a spatial Field or a data.frame/data.table with at least the columns Lon and Lat.
-#' @param unit area unit. Default m^2, if something else is spefified udunits2 must be installed. If udunits2 is not installed this option is ignored.
+#' @param unit area unit. Default "m^2", can also be "km^2" and "ha"
 #' @param ellipse If the eath should be assumed to be a ellipsoid instead of a sphere.
 #' @param digits Numeric, number of digits to which to truncate the coordinates when merging the area data.table with the input. This is a technical detail,
 #' you only need to use it if you have troubles because of coordinates with a few decimal places.   
 #' @param verbose print some information.
+#' 
+#' The main use of this function is to calculate gridcell areas internally for gridcells weighted sums and averages in \code{aggregateSpatial} but it can
+#' be utilised by the user for any other purpose.
+#' 
+#' 
 #' @export
 #' @return same class as input
 #' @author Joerg Steinkamp \email{joerg.steinkamp@@senckenberg.de}
+#' 
+#' @examples 
+#' 
+#' \donttest{
+#'  
+#' # Get an example Field
+#' africa.dir <- system.file("extdata", "LPJ-GUESS_Runs", "CentralAfrica", package = "DGVMTools")
+#' africa.Source <- defineSource(name = "LPJ-GUESS", dir = africa.dir,  format = GUESS)
+#' field <- getField(source = africa.Source, var = "cmass", year.aggregate.method="mean")
+#' 
+#' # add area in m^2 - note "Area" column follows immediately after Lon and Lat
+#' field.m2 <- addArea(input = field, unit = "m^2")
+#' print(field.m2@data)
+#' 
+#' # add area in km^2
+#' field.km2 <- addArea(input = field, unit = "km^2")
+#' print(field.km2@data)
+#' 
+#' # plot area (also modify plot and legend title to be more meaningful)
+#' p <- plotSpatial(field.km2, "Area", title = "Gridcell Area (km^2)", subtitle = NULL)
+#' p <- p + guides(fill = guide_colourbar(title = "km^2"))
+#' print(p)
+#' 
+#' # add area in km^2 using an ellipse
+#' field.km2.ellipse <- addArea(input = field, unit = "km^2", ellipse = TRUE)
+#' field.km2.ellipse <- renameLayers(field.km2.ellipse, "Area", "Area_ellipse")
+#' 
+#' # compare areas just for fun
+#' comp.layer <- compareLayers(field1=field.km2, field2=field.km2.ellipse, 
+#'                             layers1="Area", layers2="Area_ellipse")
+#' plot.title <- "Difference in gridcell area for spherical vs ellipsoid Earth"
+#' p <- plotSpatialComparison(comp.layer, title= plot.title, subtitle = NULL)
+#' p <- p + guides(fill = guide_colourbar(title = "km^2"))
+#' print(p)
+#' 
+#' }
 addArea <- function(input, unit="m^2", ellipse=FALSE, verbose=TRUE, digits = 10) {
   ## to avoid "no visible binding for global variable" during check
-  Lat = Lon = NULL
+  Lat = Lon = Area = NULL
   if (is.na(unit))
     unit="m^2"
 
@@ -164,22 +205,11 @@ addArea <- function(input, unit="m^2", ellipse=FALSE, verbose=TRUE, digits = 10)
   }
 
   if (unit!="m^2") {
-    if (requireNamespace("udunits2", quietly=TRUE)) {
-      if (udunits2::ud.is.parseable(unit)) {
-        if (udunits2::ud.are.convertible("m^2", unit)) {
-          area$area = udunits2::ud.convert(area$area, "m^2", unit)
-        } else {
-          warning(paste("m^2 not convertible to '", unit, "'. Using m^2 instead.", sep=""))
-          unit="m^2"
-        }
-      } else {
-        warning(paste("Unit '", unit, "' not parseable! Using m^2 instead.", sep=""))
-        unit="m^2"
-      }
-    } else {
-      message("Package 'udunits2' not installed! Using m^2 instead.")
-      warning("Package 'udunits2' not installed! Using m^2 instead.")
-    }  
+    
+    if(unit == "km^2") area[, Area := Area / 10^6]
+    else if(unit == "ha") area[, Area := Area / 10^4]
+    else stop(paste("Unsupported unit string in addArea", unit))
+    
   }
 
 

R/aggregateSpatial.R

@@ -27,7 +27,7 @@ aggregateSpatial.uncompiled <- function(input.obj,
                                         ...){
 
   # Messy solution to stop "notes" about undeclared global variables stemming from data.table syntax 
-  Year = Lat = Lon = area = NULL
+  Year = Lat = Lon = Area = NULL
 
   ### SET UP THE AGGREGATE METHOD 
 
@@ -71,19 +71,18 @@ aggregateSpatial.uncompiled <- function(input.obj,
 
   # Do the aggregating
   if (method == "w.mean") {
-    if (!any(colnames(input.dt)=="area")) {
+    if (!any(colnames(input.dt)=="Area")) {
       if (verbose) message("Add column area.")
       input.dt <- addArea(input.dt, verbose=verbose, ...)
     }
     if(verbose) message(paste("Spatially averaging (with area weighting) ...", sep = ""))
 
-    
-    output.dt <- input.dt[,lapply(.SD, method.function, w=area), by=by.dims]
-    output.dt[,area:=NULL]
+    output.dt <- input.dt[,lapply(.SD, method.function, w=Area), by=by.dims]
+    output.dt[,Area:=NULL]
 
   } 
   else if (method == "w.sum") {
-    if (!any(colnames(input.dt)=="area")) {
+    if (!any(colnames(input.dt)=="Area")) {
       if (verbose) message("Add column area.")
       input.dt <- addArea(input.dt, verbose=verbose)
     }
@@ -97,10 +96,10 @@ aggregateSpatial.uncompiled <- function(input.obj,
     col.names <- col.names[!col.names %in% remove]
 
     # check to see if Year is still a column name (it might have been averaged away)
-    input.dt[, (col.names) := lapply(.SD, function(x) x * input.dt[['area']] ), .SDcols = col.names]
+    input.dt[, (col.names) := lapply(.SD, function(x) x * input.dt[['Area']] ), .SDcols = col.names]
 
     output.dt <- input.dt[, lapply(.SD, method.function), by=by.dims]
-    output.dt[,area:=NULL]
+    output.dt[,Area:=NULL]
 
   }