initiated package and added functions

DESCRIPTION

+Package: clews_biosphere
+Type: Package
+Title: Functions for the CLEWS - Intro to the terrestrial biosphere - Exercise
+Version: 1.0.0
+Author: Fabian Stenzel
+Maintainer: Fabian Stenzel <stenzel@pik-potsdam.de>
+Description: Functions to read and plot netcdf climate data, and compare NPP 
+    between the Miami model and measurements by Zhang et al. 2018
+License: GPL-3
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 7.2.3
+Depends:
+    RColorBrewer,
+    ncdf4

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+export(LPJ_datatype)
+export(LpjmlMetaData)
+export(aggregateLPJmLdata2Country)
+export(asub)
+export(autoReadInput)
+export(autoReadMetaOutput)
+export(autoReadclmOutput)
+export(auto_read_nc)
+export(calcCellarea)
+export(calcEFRs)
+export(create_header)
+export(daily2monthly)
+export(evenGrid2lpjGrid)
+export(exportVars2GlobalEnv)
+export(flowSeason)
+export(freadheader)
+export(fwriteheader)
+export(get_datatype)
+export(headersize)
+export(lpj2raster)
+export(lpjCOW)
+export(lpj_datatype)
+export(mapBetweenLPJgrids)
+export(monthly2yearly)
+export(netcdfCFT2lpjarray)
+export(netcdfDaily2lpjarray)
+export(netcdfMonthly2lpjarray)
+export(netcdfYearly2lpjarray)
+export(new_header)
+export(pie2)
+export(plotBasin2screen)
+export(plotCountryData)
+export(plotGlobal)
+export(plotGlobalFlex)
+export(plotGlobalMan)
+export(plotGlobalManToScreen)
+export(plotGlobalToScreen)
+export(plotLUshares)
+export(plotLUsharesToScreen)
+export(plotLocations)
+export(plotLocationsToScreen)
+export(plotMonth)
+export(plotRegionalW)
+export(plotTimelines)
+export(plotToScreen)
+export(plot_relDiff)
+export(printLegend)
+export(readCFTinput)
+export(readCFToutput)
+export(readDaily)
+export(readGridInputBin)
+export(readGridOutputBin)
+export(readGridOutputNC)
+export(readMonthly)
+export(readMonthlyCFToutput)
+export(readNCoutput)
+export(readReservoirInput)
+export(readYearly)
+export(readYearlyInput)
+export(read_header)
+export(read_meta)
+export(readheader2)
+export(removeHighFlows)
+export(replace_array)
+export(showRoute)
+export(subset_array)
+export(upstreamCells)
+export(writeNETCDFfromLPJ)
+export(writeWUinput)
+export(write_climate_input)
+export(write_landuse_input)
+export(writeheader)

R/clews.R

+#' Get time resolution of a netcdf file
+#'
+#' Utility function to determine time resolution of a netcdf file
+#'
+#' @param file_nc netcdf file name
+#'
+#' @return character string: determined timestep (daily, monthly, yearly)
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+get_timestep <- function(file_nc){
+  tunit <- file_nc$dim$time$units
+  tvals <- file_nc$dim$time$vals
+  if (grepl("years since", tunit, ignore.case = TRUE)){
+    tres <- "annual"
+    firstyr <- as.character( as.integer(
+      unlist(strsplit(unlist(strsplit(tunit, split = ' ', 
+                                      fixed = TRUE))[3], split = '-', fixed = TRUE))[1]) + 
+        tvals[1] )
+  }else if (grepl("year", tunit, ignore.case = TRUE)){
+    tres <- "annual"
+    firstyr <- tvals[1]
+  }else if (grepl("days since", tunit, fixed = TRUE)){
+    ddiff <- tvals[2]-tvals[1]
+    firstyr <- as.character( as.integer(
+      unlist(strsplit(unlist(strsplit(tunit, split = ' ', 
+                                      fixed = TRUE))[3], split = '-', fixed = TRUE))[1]) + 
+        floor(tvals[1]/365) )
+    if (ddiff > 27 && ddiff < 32){
+      tres <- "monthly"
+    }else if (ddiff > 364 && ddiff < 367){
+      tres <- "annual"
+    }else if (ddiff == 1){
+      tres <- "daily"
+    }else{
+      stop("Automatic detection of firstyear and time resolution failed.")
+    }
+  }else{
+    stop("Automatic detection of firstyear and time resolution failed.")
+  }
+  return(list(tres = tres, firstyr = as.numeric(firstyr)))
+}
+
+#' Get main variable of a netcdf file
+#'
+#' Utility function to guess the main variable of a netcdf file
+#'
+#' @param file_nc netcdf file name
+#'
+#' @return character string: determined main variable name
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+get_main_variable <- function(file_nc){
+  for (var in names(file_nc$var)){
+    ndims <- file_nc$var[[var]]$ndims
+    dimNames <- c()
+    for (d in 1:ndims){
+      dimNames <- append(dimNames,file_nc$var[[var]]$dim[[d]]$name) 
+    }
+    #print(paste0(var,paste(dimNames,collapse = ",")))
+    if (grepl("lon",paste(dimNames,collapse = " "),ignore.case = T) && grepl("lat",paste(dimNames,collapse = " "),ignore.case = T)){
+      return(var)
+    }
+  }
+  print(paste0("None of the variables could certainly be identified as main variable, guessing the last one: ",var))
+  return(var)
+}
+
+#' Reads netcdf and returns the requested variable as an array
+#'
+#' Reads an arbitrary netcdf and returns the requested variable in the given 
+#' year range as an array
+#'
+#' @param nc_in_file netcdf file name
+#' @param var optional variable to be read, in case automatic detection does 
+#'        not work as intended or several variables are stored within the file
+#' @param get_year_start first year to be read 
+#'        (if not specified will default to first record year)
+#' @param get_year_stop final year to be read 
+#'        (if not specified will default to last record year)
+#' @param suppress_read_print whether to suppress the info "reading file XYZ" 
+#'        (default FALSE)
+#'
+#' @return array with netcdf's data, dim=c([longitude],[latitude],[bands],[months],[years])
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+read_netcdf <- function(nc_in_file,
+                        var = NULL,
+                        get_year_start = NULL,
+                        get_year_stop = NULL,
+                        suppress_read_print = FALSE) {
+  
+  file_nc <- ncdf4::nc_open(filename = nc_in_file)
+  if (is.null(var)) var <- get_main_variable(file_nc = file_nc)
+  if (!suppress_read_print) {
+    print(paste0("Reading in: ", nc_in_file))
+    print(paste0("Attempting to read variable: ",var,
+                 ". If this is not correct,",
+                 " please correct via argument var."))
+  }
+  
+  ndims <- file_nc$var[[var]]$ndims
+  tunit <- file_nc$dim$time$units
+  tvals <- file_nc$dim$time$vals
+  timing <- get_timestep(file_nc = file_nc)
+  tres <- timing$tres
+  firstyr <- timing$firstyr
+  if (ndims == 3){
+    data <- ncdf4::ncvar_get(nc = file_nc, varid = var,
+                             start = c(1,1,1), count = c(1,1,-1))
+    data_dim <- dim(data)
+    nbands <- 1
+    tsteps <- data_dim[1]
+  }else if (ndims >= 4){
+    data <- ncdf4::ncvar_get(nc = file_nc, varid = var,
+                             start = c(1,1,1,1), count = c(1,1,-1,-1))
+    data_dim <- dim(data)
+    nbands <- data_dim[1]
+    tsteps <- data_dim[2]
+  }else{
+    stop("Less than 3 dimensions for file. aborting.")
+  }
+  if (tres == "annual"){
+    nyears <- tsteps
+    nmonths <- 1
+  }else if (tres == "monthly"){
+    nyears <- tsteps/12
+    nmonths <- 12
+  }
+  if (is.null(get_year_start)) get_year_start <- firstyr
+  if (is.null(get_year_stop)) get_year_stop <- firstyr + nyears - 1
+  ngetyears <- get_year_stop - get_year_start + 1
+  
+  nlatin <- file_nc$dim$lat$len
+  nlonin <- file_nc$dim$lon$len
+  lat_values <- file_nc$dim$lat$vals
+  lon_values <- file_nc$dim$lon$vals
+  outdata <- array(0, dim = c(nlatin,nlonin,nbands,nmonths,ngetyears))
+  
+  # get spatial extent and resolution
+  # this will give a warning if the NetCDF has more than one data field, 
+  # e.g. crop bands or time axis
+  for (year in get_year_start:get_year_stop){
+    for (month in 1:nmonths){
+      if (nbands == 1){
+        data <- ncdf4::ncvar_get(nc = file_nc, varid = var, count=c(-1,-1,1),
+                                 start=c(1,1,((year - firstyr)*nmonths + month)))
+        # check whether data needs to be flipped vertically
+        # print(paste0("lat: ",lat_values[1],", ",lat_values[2]))
+        if (lat_values[1]>lat_values[2]){
+          outdata[,,1,month,(year - get_year_start + 1)] <- data[,nlatin:1]
+          lat_values <- rev(lat_values)
+        }else{
+          outdata[,,1,month,(year - get_year_start + 1)] <- data
+        }
+      }else{ #nbands>1
+        data <- ncvar_get(nc = file_nc, varid = var, count=c(-1,-1,-1,1),
+                          start=c(1,1,1,((year - firstyr)*nmonths + month)))
+        # check whether data needs to be flipped vertically
+        if (lat_values[1]<lat_values[2]){
+          outdata[,,,month,(year - get_year_start + 1)] <- data[,nlatin:1,]
+          lat_values <- rev(lat_values)
+        }else{
+          outdata[,,,month,(year - get_year_start + 1)] <- data
+        }
+      }# end if nbands == 1
+    }# end for month in 1:nmonths
+  }# end for year in ...
+  ncdf4::nc_close(file_nc)
+  dim(outdata) <- c(lon = nlonin, lat = nlatin, band = nbands, 
+                    month = nmonths, year = ngetyears)
+  dimnames(outdata) <- list(lon = lon_values, lat = lat_values, band = 1:nbands, 
+                            month = 1:nmonths, year = get_year_start:get_year_stop)
+  return(drop(outdata))
+}
+
+#' Plot an array of lon_lat data to screen or file 
+#'
+#' Plot an array of lon_lat data to screen or file 
+#' e.g. obtained from read_netcdf and averaged over time.
+#' Prints to screen if file argument is not supplied (default).
+#'
+#' @param data array to plot. Needs to be in format [longitude,latitude].
+#' @param file character path. file location to save the plot to. 
+#'        If not supplied prints to screen. Default: NULL
+#' @param title character string printed as title.
+#'
+#' @return None
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+plot_lon_lat_array <- function(data, file = NULL, title = ""){
+  di <- dim(data)
+  if (length(di)>2) stop("Too many dimensions on data object, 
+                         please reduce by picking/summing/averaging.")
+  par(mar=c(3,3,0,0)) #bltr
+  ra <- terra::rast(t(data[,di[2]:1]))
+  range <- range(data)
+  extent <- terra::ext(c(0, di[1], 30, di[2]))
+  if (!is.null(file)) png(file, width=7.25, height=3.5, units="in", res=300, pointsize=6,type="cairo")
+  terra::plot(ra, main = title,ext = extent)
+  if (!is.null(file)) dev.off()
+}
+
+#' Returns cellarea of given netcdf in m
+#'
+#' Returns approximate cellarea of given netcdf file in m
+#'
+#' @param nc_in_file netcdf file name
+#'
+#' @return cellarea array with dimensions [longitude, latitude] in m
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+get_netcdf_cellsize <- function(nc_in_file = tmp_file){
+  file_nc <- ncdf4::nc_open(filename = nc_in_file)
+  lat_values <- file_nc$dim$lat$vals
+  lon_values <- file_nc$dim$lat$vals
+  res.lat <- abs(lat_values[1] - lat_values[2])
+  res.lon <- abs(lon_values[1] - lon_values[2])
+  nc_close(file_nc)
+  
+  earthradius <- 6371000.785 # in m
+  cellwidth <- earthradius * 2 * pi / 360 # cellwidth per degree at equator in m
+  # cells are approx. rectangular at the equator - they get smaller towards 
+  # the poles (scaling with the cos of the midpoint latitude - here in radians)
+  cellarea_raw <- (cellwidth * res.lon) * (cellwidth * res.lat) *
+    cos(lat_values/180 * pi) # in m
+  # we want the cellarea in the same array format as tmp and pre
+  cellarea <- rep(x = cellarea_raw, each = length(tmp_lon))
+  dim(cellarea) <- c(length(tmp_lon), length(tmp_lat))
+  return(cellarea)
+}
+
+#' Plot the climate of a given location
+#'
+#' Plot an array of lon_lat data to screen or file 
+#' e.g. obtained from read_netcdf and averaged over time.
+#' Prints to screen if file argument is not supplied (default).
+#'
+#' @param tmp temperature array with 12 values for each month in °C
+#' @param pre precipitation array with 12 values for each month in mm
+#' @param file character path. file location to save the plot to. 
+#'        If not supplied prints to screen. Default: NULL
+#' @param title character string printed as title.
+#'
+#' @return None
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+climate_plot <- function(temp, prec, file = NULL, title = ""){
+  par(mar=c(4,4,1,4),oma=c(0,0,0,0)) #bltr
+  if (!is.null(file)) png(file, width=7.25, height=3.5, 
+                          units="in", res=300, pointsize=6,type="cairo")
+  # plot average monthly precipitation for past 30 years
+  x_pos <- barplot(height = prec,col = "blue", ylab = "prec in mm",
+                   names.arg = c("J","F","M","A","M","J","J","A","S","O","N","D"))
+  # plot the average temp. of the past 30 years
+  par(fig=c(0,1,0,1))
+  lines(x = x_pos, y = temp, ylim = range(temp),xlab = "Month", ylab = "",col = "red", lwd = 2)
+  axis(side = 4, col = "red")
+  mtext(side = 4, line = 2.5, text = "temp in °C")
+  if (!is.null(file)) dev.off()
+}
+
+#' Read in the Zhang et al. csv file with NPP measurements
+#'
+#' Read in the Zhang et al. csv file with NPP measurements and return a 
+#' list of relevant records
+#'
+#' @param csv_in_file character string. full csv file path
+#'
+#' @return list object with dimensions [lat, lon, totnpp, year] per record
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+read_in_zhang_file <- function(csv_in_file){
+  npp_measurements <- read.csv(file = csv_in_file)
+  return(npp_measurements[,c("latitude","longitude","totnpp","year")])
+}
+
+#' Add modelled NPP to measured zhang data
+#'
+#' Add the best fitting modelled NPP to the measured zhang data list object 
+#' (e.g. obtained from read_in_zhang_file() function)
+#'
+#' @param zhang_data zhang data list object 
+#'        (e.g. obtained from read_in_zhang_file() function)
+#' @param modelled_npp array with modelled NPP (dimensions [lat,lon,year])
+#'
+#' @return modified zhang_data input list, with added colums: measured, modelled
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+add_modelled_npp_to_zhang_npp_measurements <- function(zhang_data, modelled_npp) {
+  min_model_year <- min(as.integer(dimnames(modelled_npp)$year))
+  years <- unique(zhang_data$year)
+  #      Average,1991,1959-98,1988-90,Potential,1993,1986-93,1982-89,1983-90,1992,1982-96,1988,latest interval,1990-93,1987-97
+  from <- c(  51,  91,     59,     88,       NA,  93,     86,     82,     83,  92,     82, 88,              NA,     90,     87) - (2000-min_model_year)
+  to <-   c( 122,  91,     98,     90,       NA,  93,     93,     89,     90,  92,     96, 88,              NA,     93,     97) - (2000-min_model_year)
+  npp_measured <- zhang_data$totnpp
+  npp_modelled <- array(NA,dim=length(npp_measured))
+  for (i in 1:length(npp_measured)){
+    if (zhang_data$year[i] %in% c("Potential","latest interval") 
+        | zhang_data$totnpp[i]<0 
+        | is.na(zhang_data$totnpp[i])) {
+      npp_measured[i] <- NA
+      npp_modelled[i] <- NA
+      next
+    }
+    latPick <- round((zhang_data$latitude[i] - 0.25) * 2) / 2 + 0.25
+    lonPick <- round((zhang_data$longitude[i] - 0.25) * 2) / 2 + 0.25
+    ind <- which(years == zhang_data$year[i])
+    npp_modelled[i] <- mean(NPP[as.character(lonPick),as.character(latPick),from[ind]:to[ind]])
+  }
+  zhang_data$measured <- npp_measured
+  zhang_data$modelled <- npp_modelled
+  return(zhang_data)
+}
+
+#' Create a scatterplot
+#'
+#' Create a scatterplot of two data series e.g. measured vs. modelled
+#'
+#' @param x data to plot on x axis
+#' @param y data to plot on y axis
+#' @param file character path. file location to save the plot to. 
+#'        If not supplied prints to screen. Default: NULL
+#' @param title character string printed as title.
+#' @param xlab label for x axis
+#' @param ylab label for y axis
+#'
+#' @return None
+#'
+#' @examples
+#' \dontrun{
+#' }
+#'
+#' @export
+scatter_plot <- function(x, y, file = NULL, title = "", xlab = "", ylab = ""){
+  if (!is.null(file)) png(file, width=7.25, height=3.5, 
+                          units="in", res=300, pointsize=6,type="cairo")
+  par(mar=c(4.5,4.5,1,1),bty = "o")
+  max = max(y,x, na.rm = T)
+  plot(x = x, y = y, ylim = c(0,max), xlim = c(0,max), xlab = xlab, 
+       ylab = ylab, asp=1, main = title)
+  abline(a = 0,b = 1)
+  fit <- summary(lm(x~y))
+  text(x = 0, y = max*0.9,paste("R^2 =",round(fit$r.squared,4)),adj=c(0,1))
+  if (!is.null(file)) dev.off()
+}
\ No newline at end of file

README.md

+## Installation
+
+You can install `clews_biosphere` by git cloning this repository:
+
+```bash
+git clone https://gitlab.pik-potsdam.de/stenzel/clews_biosphere.git <path_to_clews_biosphere>
+```
+(replace <path_to_clews_biosphere> with the location you want to save it to)
+
+and install via  [`devtools`](https://rawgit.com/rstudio/cheatsheets/master/package-development.pdf):
+
+```R
+devtools::install("<path_to_clews_biosphere>")
+library("clews_biosphere")
+```
+
+alternatively, you can also load it from source without installing it:
+
+```R
+devtools::load_all("<path_to_clews_biosphere>")
+```
+
+## Examples

clews.Rproj

+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace,vignette