#'@name surfplot
#'@title Plotting Interpolated Surfaces
#'@description accounts for corner case parameter spellings, variable specific contour breaks are (need to be) defined
#'@author Jemma Stachelek
#'@param rnge numeric string of no more than two dates in yyyymm format
#'@param params character. string of parameter fields to plot
#'@param fdir character file path to local data directory
#'@param yext numeric length 2 y extent
#'@param xext numeric length 2 x extent
#'@return output plots to plot window
#'@importFrom rgdal readOGR
#'@importFrom rasterVis levelplot
#'@importFrom latticeExtra layer
#'@importFrom sp spplot
#'@importFrom raster raster stack reclassify calc writeRaster
#'@export
#'@examples \dontrun{
#'surfplot(rnge = c(200707), params = c("cdom"))
#'surfplot(201513, c("ph", "c6turbidity", "c6chl", "c6cdom"),
#' yext = c(2786102, 2797996), xext = c(557217, 567415))
#'}
surfplot <- function(rnge = c(201402, 201404), params = c("c6chl", "sal"), fdir = getOption("fdir"), yext = c(2772256, 2798000), xext = c(518000.2, 566000)){
#print(fdir)
#fbcoast <- rgdal::readOGR(dsn = file.path(fdir, "DF_Basefile", "FBcoast_big.shp"), layer = "FBcoast_big", verbose = TRUE)
if(length(rnge) == 1){
rnge <- c(rnge, rnge)
}
namesalias <- read.table(text = "
chlorophyll.a c6chl
c6chla c6chl
")
#define breaks
brks <- read.table(text = "
sal list(seq(0,40,2))
salinity.pss list(seq(0,40,2))
salpsu list(seq(0,40,2))
c6chl list(seq(50,200,10))
chlext list(seq(0,5,0.5),seq(10,30,5))
chlext_hi list(seq(0,5,0.5),seq(10,30,5))
chlext_low list(seq(0,5,0.5),seq(10,30,5))
temp list(seq(14,36,2))
c6temp list(seq(14,36,2))
c6cdom list(seq(80,360,40))
ph list(seq(7.3,8.1,0.1))
c6turbidity list(seq(0,45,5))")
dirlist <- list.dirs(file.path(fdir, "DF_Surfaces"), recursive = F)
minrnge <- min(which(substring(basename(dirlist), 1, 6) >= rnge[1]))
maxrnge <- max(which(substring(basename(dirlist), 1, 6) <= rnge[2]))
rlist <- list.files(dirlist[minrnge:maxrnge], full.names = T, include.dirs = T, pattern = "\\.tif$")
plist <- tolower(sub("[.][^.]*$", "", basename(rlist)))
for(n in 1:length(plist)){
if(any(plist[n] == namesalias[,1])){
plist[n] <- as.character(namesalias[which(plist[n] == namesalias[,1]), 2])
#print(names(dt)[n])
}
}
rlist <- rlist[which(!is.na(
match(plist, params)
))]
plist <- plist[which(!is.na(match(plist, params)))]
for(i in 1:length(rlist)){
#i<-1
my.at <- unlist(eval(parse(text = as.character(brks[which(plist[i] == brks[,1]), 2]))))
print(rasterVis::levelplot(raster::raster(rlist[i]),ylim =yext, xlim = xext,par.settings=rasterVis::PuOrTheme(),at=my.at,margin=FALSE,auto.key=FALSE,scales=list(draw=FALSE),main=paste(as.character(plist[i]),unlist(strsplit(rlist[i],"/"))[length(unlist(strsplit(rlist[i],"/")))-1]))+latticeExtra::layer({sp::SpatialPolygonsRescale(sp::layout.north.arrow(),offset=c(563000,2775000),scale=4400)})+ latticeExtra::layer(sp::sp.polygons(rgdal::readOGR(dsn = file.path(getOption("fdir"), "DF_Basefile", "FBcoast_big.shp"), layer = "FBcoast_big", verbose = TRUE), fill="green",alpha=0.6)))
}
}
#'@name avmap
#'@title Create a difference map compared to average
#'@param yearmon survey of interest to compare against average
#'@param params variable name generally choice of "sal" or "chlext"
#'@param diffpath file.path to write output difference surface
#'@param avpath file.path to write output average surface
#'@param percentcov numeric account for the different raster extents by setting the percent of all surveys required before a pixel is included in difference from average computations
#'@param tolerance numeric number of months on either side of yearmon to include in the set of surfaces averaged . defaults to 1.
#'@param fdir character file path to local data directory
#'@description takes a survey date as input and searches the DF_Surfaces folder for maps of the same parameter within the range of 1-2 months from yearmon for each year. These surfaces are averaged and compared to the surface from yearmon.
#'@details resulting surfaces are written to disk if the avpath or diffpath arguments are specified
#'@export
#'@importFrom raster raster stack reclassify calc writeRaster
#'@examples \dontrun{
#'avmap(yearmon = 201505, params = "sal", diffpath = file.path(getOption("fdir"),
#' "DF_Surfaces", 201505, paste0("diff", "sal", ".tif")),
#' percentcov = 0.6, tolerance = 1)
#'avmap(yearmon = 201505, params = "sal")
#'}
avmap <- function(yearmon, params, diffpath = NULL, avpath = NULL, percentcov = 0.6, tolerance = 1, fdir = getOption("fdir")){
#generate-flist=========================================================#
flist.full <- list.files(file.path(fdir, "DF_Surfaces"), pattern = "*.grd", recursive = T, include.dirs = T, full.names = T)
flist <- flist.full[basename(flist.full) == paste(toupper(params), ".grd", sep = "") | basename(flist.full) == paste(tolower(params), ".grd", sep = "")]
sdates <- data.frame(matrix(unlist(strsplit(dirname(flist), "/")), nrow = length(flist), byrow = T))
sdates <- substring(sdates[,ncol(sdates)], 1, 6)
cursurf <- raster::raster(flist[which(sdates == yearmon)])
flist <- flist[-which(sdates == yearmon)]
sdates <- sdates[-which(sdates == yearmon)]
curmon <-as.numeric(substring(yearmon, 5, 6))
flist <- flist[as.numeric(substring(sdates, 5, 6)) <= curmon + tolerance & as.numeric(substring(sdates, 5, 6)) >= curmon - tolerance]
remlist <- which(is.na(flist))
if(length(remlist) > 0){
flist <- flist[-remlist]
}
#stack-and-calculate====================================================#
rmean <- average_rlist(flist, percentcov)
rmean_diff <- cursurf - rmean
#plotting===============================================================#
sdates <- data.frame(matrix(unlist(strsplit(dirname(flist), "/")), nrow = length(flist), byrow = T))
sdates <- substring(sdates[,ncol(sdates)], 1, 6)
sp::plot(rmean, main = paste("Average", params, sdates[1], "-", sdates[length(sdates)], sep = " "))
sp::plot(rmean_diff, main = "Difference from Average")
#save-to-file===========================================================#
if(length(avpath) > 0){
raster::writeRaster(rmean, avpath, format = "GTiff", overwrite = T)
}
if(length(diffpath) > 0){
raster::writeRaster(rmean_diff, diffpath, format = "GTiff", overwrite = T)
}
rmean_diff
}
#'@name create_rlist
#'@title Create file listing of rasters from a date range and parameter name
#'@description Create file listing of rasters from a date range and parameter name
#'@param rnge numeric string of 1, 2, or more dates dates in yyyymm format. A length 1 rnge will produce a single plot, a length 2 rnge will produce a series of plots bookended by the two dates, a rnge object with more than 2 dates will produce a series of plots exactly corresponding to the dates provided.
#'@param params character vector of parameter fields to plot legends and color ramps are defined for sal, chlext, and diffsal
#'@export
#'@examples \dontrun{
#'create_rlist(rnge = 201509, params = "sal")
#'create_rlist(rnge = c(201509, 201605), params = c("sal", "salinity.pss", "salpsu"))
#'create_rlist(rnge = c(200808, 200910, 201002, 201004, 201007, 201102, 201105,
#' 201206, 201209, 201212, 201305, 201308, 201311, 201404, 201407, 201410,
#' 201502, 201505, 201507, 201509), params = 'chlext')
#'
# rlist <- create_rlist(good_years, "chlext")
#'}
create_rlist <- function(rnge, params){
fdir <- getOption("fdir")
if(length(rnge) == 1){
rnge <- c(rnge, rnge)
}
namesalias <- read.table(text = "
chlorophyll.a c6chl
c6chla c6chl
")
dirlist <- list.dirs(file.path(fdir, "DF_Surfaces"), recursive = F)
minrnge <- min(which(substring(basename(dirlist), 1, 6) >= rnge[1]))
maxrnge <- max(which(substring(basename(dirlist), 1, 6) <= rnge[2]))
datelist <- sapply(list.dirs(dirlist[minrnge:maxrnge]), function(x) substring(x, nchar(x) - 5))
rlist <- list.files(dirlist[minrnge:maxrnge], full.names = T, include.dirs = T, pattern = "\\.tif$")
plist <- tolower(sub("[.][^.]*$", "", basename(rlist)))
for(n in 1:length(plist)){
if(any(plist[n] == namesalias[,1])){
plist[n] <- as.character(namesalias[which(plist[n] == namesalias[,1]), 2])
}
}
if(length(rnge) > 2){
rnge <- rnge[order(rnge)]
rlist <- list.files(dirlist[which(substring(basename(dirlist), 1, 6) %in% rnge)], full.names = T, include.dirs = T, pattern = "\\.tif$")
plist <- tolower(sub("[.][^.]*$", "", basename(rlist)))
}
rlist <- rlist[which(!is.na(match(plist, params)))]
plist <- plist[which(!is.na(match(plist, params)))]
list(rlist = rlist, plist = plist, datelist = datelist)
}
#'@name grassmap
#'@title Publication quality maps with GRASS
#'@description Create publication quality maps with GRASS GIS. Only runs on Linux.
#'@author Jemma Stachelek
#'@param fpath file.path to geotiff file
#'@param rnge numeric string of 1, 2, or more dates dates in yyyymm format. A length 1 rnge will produce a single plot, a length 2 rnge will produce a series of plots bookended by the two dates, a rnge object with more than 2 dates will produce a series of plots exactly corresponding to the dates provided.
#'@param params character vector of parameter fields to plot legends and color ramps are defined for sal, chlext, phycoc, and diffsal
#'@param fdir character file path to local data directory
#'@param cleanup logical remove intermediate rasters and shapefiles?
#'@param rotated logical rotate canvas to fit Florida Bay more squarely? This requires the i.rotate extension to be installed and addons configured (not working).
#'@param labelling logical lablel output with yearmon?
#'@param numrow numeric number of rows in a multipanel
#'@param numcol numeric number of columns in a multipanel
#'@param mapextent numeric vector of length 4. Format is bottomleft-x, topright-x, bottomleft-y, topright-y.
#'@param basin character basin name from fathom_basins_proj.shp
#'@param label_string character label text
#'@param label_size int label size
#'@param print_track logical print dataflow track?
#'@return output plots to the QGIS_plotting folder
#'@details Probably need to implement this as a seperate package to improve portability. Set param to "diffsal" to plot outpot of avmap function. Will output an imagemagick plot to the working directory and a pdf plot to the file.path(getOption("fdir"), "QGIS_Plotting") folder. The optional fpath argument only supports pointing to a single geotiff. This function can be called from the commandline using Rscript by loading the methods package and creating/loading a python 2 environment using the commands conda create --name python2 python=2 and source activate python2.
#'@import rgrass7
#'@import maptools
#'@import rgeos
#'@import methods
#'@export
#'@examples \dontrun{
#'#list supported parameters
#'grassmap(rnge = c(201512), params = c("sal"))
#'grassmap(rnge = c(201512), params = c("diffsal"))
#'grassmap(rnge = c(201407), params = c("chlext"))
#'
#'#multiple survey dates
#'grassmap(rnge = c(201509, 201512), params = c("sal"))
#'
#'#exclude Card Sound
#'grassmap(rnge = c(201507), params = c("sal"), mapextent = c(494952.6, 564517.2, 2758908, 2799640))
#'
#'#print survey track and zoom to Manatee + Barnes
#'grassmap(201513, "chlext", mapextent = c(557217, 567415, 2786102, 2797996), print_track = TRUE)
#'
#'grassmap(rnge = 201512, params = "sal", mapextent = c(557217,
#' 567415, 2786102, 2797996), label_string = "2015-12-01")
#'
#'grassmap(rnge = c(201512), params = c("sal"), basin = "Manatee Bay")
#'
#'#specify raster file directly
#'grassmap(fpath = file.path(getOption("fdir"), "DF_Surfaces", "200904", "sal.tif"), params = "sal")
#'
#'#specify label string
#'grassmap(fpath = file.path("salpre.proj_mean.tif"), params = "sal", label_string = "Pre C-111")
#'
#'#create a new color ramp by editing DF_Basefile/*.file and update figure makefile
#'logramp(n = 8, maxrange = 20) #chlext
#'scales::show_col(viridis::viridis_pal()(9))
#'logramp(n = 8, minrange = 10, maxrange = 65) #phycoc
#'
#'grassmap(rnge = c(201509, 201512), params = "sal", numrow = 2, numcol = 1, labelling = FALSE)
#'
#'grassmap(rnge = 201509, params = "phycoc", label_string = "phycoc test")
#'
#'grassmap(rnge = c(200808, 200910), params = 'chlext', numrow = 2, numcol = 1,
#' labelling = TRUE, label_string = c("Aug 2008", "Oct 2009"))
#'
#'
# grassmap(rnge = 201509, params = 'chlext', labelling = TRUE, label_size = 35, label_string = "Sep 2015")
#'}
grassmap <- function(fpath = NULL, rnge = NULL , params, mapextent = NA, numrow = NULL, numcol = NULL, fdir = getOption("fdir"), basin = "full", label_string = NULL, label_size = 21, labelling = TRUE, print_track = FALSE, cleanup = TRUE, rotated = TRUE){
if(as.character(Sys.info()["sysname"]) != "Linux"){
stop("This function only works with Linux!")
}
fathombasins <- rgdal::readOGR(file.path(fdir, "DF_Basefile/fathom_basins_proj.shp"), layer = "fathom_basins_proj", verbose = FALSE)
fboutline <- rgdal::readOGR(dsn = file.path(getOption("fdir"), "DF_Basefile/FBcoast_big.shp"), layer = "FBcoast_big", verbose = FALSE)
paramkey <- read.table(text = "
sal,salrules.file
salpsu,salrules.file
salinity.pss,salrules.file
chlext,chlextrules.file
chlext_low,chlextrules.file
chlext_hi,chlextrules.file
phycoc,phycocrules.file
diffsal,diffsalrules.file",
sep = ",", stringsAsFactors = FALSE)
rulesfile <- paramkey[which(params == paramkey[,1]), 2]
if(length(fpath) == 0){
rlist <- create_rlist(rnge = rnge, params = params)$rlist
plist <- create_rlist(rnge = rnge, params = params)$plist
}else{
rlist <- fpath
plist <- rep(params, each = length(rlist))
}
if(print_track == TRUE){
surveytrack <- coordinatize(streamget(rnge[1]), latname = "lat_dd", lonname = "lon_dd")
#interp_pnts <- coordinatize(read.csv(file.path(fdir, "DF_Subsets", paste0(rnge[1], "s.csv"))), latname = "lat_dd", lonname = "lon_dd")
}
print(rlist)
for(i in 1:length(rlist)){
#set extent============================================================#
if(basin != "full"){
firstras <- raster::raster(rlist[1])
firstras <- raster::crop(firstras, fathombasins[fathombasins$NAME == basin,])
}else{
firstras <- raster::raster(rlist[i])
}
if(basin != "full"){
tempras <- raster::raster(rlist[i])
tempras <- raster::crop(tempras, fathombasins[fathombasins$NAME == basin,])
}else{
tempras <- raster::raster(rlist[i])
}
# if(!is.na(mapextent)){
# browser()
# firstras <- raster::crop(firstras, mapextent)
# }
# if(!is.na(mapextent)){
# tempras <- raster::crop(tempras, mapextent)
# }
#create raster outline====================================================#
#no fpath
if(length(fpath) == 0){
#no fpath but label_string
if(length(label_string) == 0){
label_string_single <- rasname <- paste(substring(dirname(rlist[i]), nchar(dirname(rlist[i])) - 5, nchar(dirname(rlist[i]))))
}else{
label_string_single <- label_string[i]
rasname <- paste(substring(dirname(rlist[i]), nchar(dirname(rlist[i])) - 5, nchar(dirname(rlist[i]))))
}
raspath <- file.path(paste(fdir, "/QGIS_plotting", sep = ""), paste(rasname, ".tif", sep = ""))
outpath <- file.path(paste(fdir, "/QGIS_plotting", sep = ""), paste(rasname, "poly.shp" ,sep = ""))
}else{
#fpath & label_string
label_string_single <- label_string[i]
rasname <- paste0(strsplit(basename(fpath), "\\.")[[1]][-(length(strsplit(basename(fpath), "\\.")[[1]]))], collapse = "")
raspath <- file.path(paste0(fdir, "/QGIS_plotting", sep = ""), basename(fpath))
outpath <- file.path(paste(fdir, "/QGIS_plotting", sep = ""), paste(rasname, "poly.shp" ,sep = ""))
}
#=============================================================#
raster::writeRaster(tempras, raspath, format = "GTiff", overwrite = TRUE)
shellcmds = paste("/usr/bin/gdal_polygonize.py", raspath, "-f","'ESRI Shapefile'", outpath)
system(shellcmds)
outpoly <- rgdal::readOGR(dsn = outpath, layer = paste(rasname, "poly", sep = ""), verbose = TRUE)
print(class(outpoly))
#requireNamespace("maptools")
require("rgeos")
require("maptools")#cannot seem to execute below without call to require
maptools::gpclibPermit()
outpoly <- maptools::unionSpatialPolygons(outpoly, IDs = rep(1, length(outpoly)))
outlines <- as(outpoly, 'SpatialLines')
outlines <- sp::SpatialLinesDataFrame(outlines, data = as.data.frame(1))
#GRASS block===============================================================#
loc <- rgrass7::initGRASS("/usr/lib/grass72", home = file.path(fdir, "QGIS_plotting"), override = TRUE)
#raster
firstras <- as(firstras, "SpatialGridDataFrame")
firstras.g <- rgrass7::writeRAST(firstras, "firstras", flags = c("overwrite", "quiet"))
tempras <- as(tempras, "SpatialGridDataFrame")
tempras.g <- rgrass7::writeRAST(tempras, "tempras", flags = c("overwrite"))
rgrass7::execGRASS("g.region", raster = "tempras")
if(params %in% c("chlext", "chlext_low", "chlext_hi", "phycoc")){
if(length(grep("_log", rulesfile)) != 0){
rulesfile <- gsub("_log", "", rulesfile)
}
}
rgrass7::execGRASS("r.colors", map = "tempras", rules = file.path(fdir, "DF_Basefile", rulesfile))
rgrass7::execGRASS("r.grow", input = "tempras", output = "tempras2", radius = 1.3, flags = c("overwrite", "quiet"))
#Florida Bay outline
if(is.na(mapextent)){
fboutline <- raster::crop(fboutline, raster::extent(raster::raster(tempras)))
}else{
fboutline <- raster::crop(fboutline, mapextent)
}
fbvec.g <- rgrass7::writeVECT(fboutline, "fbvec", v.in.ogr_flags = c("o"))
rgrass7::execGRASS("g.region", vector = "fbvec")
rgrass7::execGRASS("v.colors", map = "fbvec", column = "cat", color = "grey")
#survey track
if(print_track == TRUE){
trackvec.g <- rgrass7::writeVECT(surveytrack, "trackvec", v.in.ogr_flags = c("o"))
rgrass7::execGRASS("v.colors", use = "cat", map = "trackvec", color = "grey")
}
#raster outline
outvec.g <- rgrass7::writeVECT(outlines, "outvec", v.in.ogr_flags = c("o"))
rgrass7::execGRASS("g.region", vector = "outvec")
rgrass7::execGRASS("g.region", raster = "firstras")
rgrass7::execGRASS("g.region", vector = "fbvec")
if(labelling == TRUE){
if(length(label_string) == 0 & length(label_string_single) == 0){
stop("Must specify a label_string if labelling == TRUE")
}
# #compose plotting commands here####
fileConn <- file(file.path(fdir, "QGIS_plotting", "grassplot.file"))
writeLines(c("raster tempras2",
"vlines outvec",
" color black",
" style dashed",
" end",
"vareas fbvec",
" masked y",
" end",
paste("text 20% 87% ", label_string_single, sep = ""),
paste0("fontsize ", label_size),
" background white",
" border black",
" end",
"end"), fileConn)
close(fileConn)
if(print_track == TRUE){
fileConn <- file(file.path(fdir, "QGIS_plotting", "grassplot.file"))
writeLines(c("raster tempras2",
"vpoints trackvec",
" color black",
" fcolor black",
" symbol basic/cross1",
" size 5",
" end",
"vlines outvec",
" color black",
" style dashed",
" end",
"vareas fbvec",
" masked y",
" end",
paste("text 17% 85% ", label_string_single, sep = ""),
" fontsize 35",
" background white",
" border black",
" end",
"end"),fileConn)
close(fileConn)
}
}else{
fileConn <- file(file.path(fdir,"QGIS_plotting","grassplot.file"))
writeLines(c("raster tempras2",
"vlines outvec",
" color black",
" style dashed",
" end",
"vareas fbvec",
" masked y",
" end",
"end"),fileConn)
close(fileConn)
}
# label_string <- NULL
# browser()
# system("iconv -f UTF-8 -t ISO_8859-1 grassplot.file > grassplot.file")
rgrass7::execGRASS("ps.map", input = file.path(paste(fdir, "/QGIS_plotting", sep=""), "grassplot.file"), output = file.path(paste(fdir, "/QGIS_plotting", sep = ""), paste(rasname, ".pdf", sep = "")), flags = "overwrite")
#==================================================================#
legendalias <- read.table(text = "chlext,Chlorophyll (ug/L)
chlext_low,Chlorophyll (ug/L)
chlext_hi,Chlorophyll (ug/L)
sal,Salinity
salpsu,Salinity
salinity.pss,Salinity
phycoc,Phycocyanin (RFU)
diffsal,Salinity minus average", sep = ",", stringsAsFactors = FALSE)
legendname <- legendalias[which(params == legendalias[,1]), 2]
if(params %in% c("sal", "salinity.pss", "salpsu")){
paramxcoord <- 2160
legendunits<-seq(from = 5,to = 54, by = 0.1)
legendunits_print <- "'5 10 15 20 25 30 35 40'"
legendunits_spacing <- 220
legend_xlim <- 270
legend_crop_extent <- 2404
}
if(params %in% c("chlext", "chlext_low", "chlext_hi")){
paramxcoord <- 1980
legendunits <- log(seq(from = 0, to = 13.5, by = 0.1) + 1)
legendunits_print <- "'0.0 0.7 2.0 4.0 7.0 13.0'"
legendunits_spacing <- 275
if(length(grep("_log", rulesfile)) == 0){
rulesfile <- paste0(rulesfile,"_log")
}
legend_xlim <- 300
legend_crop_extent <- 2404
}
if(params %in% c("phycoc")){
paramxcoord <- 1980
legendunits <- log(seq(from = 9, to = 65, by = 0.1) + 1)
legendunits_print <- "'9.0 12.1 16.2 21.5 28.4 37.5'"
legendunits_spacing <- 275
if(length(grep("_log", rulesfile)) == 0){
rulesfile <- paste0(rulesfile,"_log")
}
legend_xlim <- 300
legend_crop_extent <- 2404
}
if(params == "diffsal"){
paramxcoord <- 1960
legendunits <- seq(from = -30, to = 35, by = 1)
legendunits_print <- "'-30 -25 -20 -15 -10 -5 0 5 10 15 20'"
legendunits_spacing <- 120
legend_xlim <- 270
legend_crop_extent <- 2404
}
#print legend========================================================#
legras <- raster::raster(tempras)
legras[1:length(legras)] <- legendunits
tempras <- as(legras, "SpatialGridDataFrame")
tempras.g <- rgrass7::writeRAST(tempras, "tempras", flags = c("overwrite"))
rgrass7::execGRASS("r.support", map = "tempras", units = legendname)
rgrass7::execGRASS("g.region", raster = "tempras")
rgrass7::execGRASS("r.colors", map = "tempras", rules = file.path(fdir, "DF_Basefile", rulesfile))
rgrass7::execGRASS("r.to.vect", input = "tempras", output = "outvec", type = "area", flags = "overwrite")
rgrass7::execGRASS("v.hull", input = "outvec", output = "outvec2", flags = "overwrite")
rgrass7::execGRASS("ps.map", input = file.path(paste(fdir, "/QGIS_plotting", sep = ""), "legendplot.file"), output = file.path(paste(fdir, "/QGIS_plotting", sep = ""), "legend", paste("legend", ".pdf", sep = "")), flags = "overwrite")
#make files================================================================#
#system(paste("echo", "'", legendname, substring(legendunits_print, 2, nchar(legendunits_print) - 1), legendunits_spacing, legend_xlim, legend_crop_extent, "'", ">> 'single.txt'"))
if(length(rlist) == 1){
makefile <- system.file("grass-image_makefiles/Makefile_single", package = "DataflowR")
system(paste0("make -f ", makefile,
" testpanel.png BASEDIR=", fdir,
" YEARMON=", rasname,
" PARAM=", shQuote(legendname),
" LEGENDUNITS=", legendunits_print,
" LEGENDUNITSSPACING=", legendunits_spacing,
" LEGEND_XLIM=", legend_xlim,
" PARAMXCOORD=", paramxcoord,
" LEGEND_CROP_EXTENT=", legend_crop_extent
))
if(cleanup == TRUE){
system(paste0("make -f ", makefile, " clean"))
}
}
#==================================================================#
if(cleanup == TRUE){
rmlist <- list.files(file.path(paste(fdir, "/QGIS_plotting", sep = "")), pattern = paste(rasname, "*", sep = ""), include.dirs = TRUE, full.names = TRUE)
rmlist <- rmlist[-grep("*.pdf", rmlist)]
file.remove(rmlist)
}
}
#==================================================================#
# system(paste(
# "echo", "'", legendname, substring(legendunits_print, 2, nchar(legendunits_print) - 1), legendunits_spacing, legend_xlim, legend_crop_extent, "'", ">> 'multi.txt'"))
#assumes that all pdfs in QGIS_plotting are to be part of panel
if(length(rlist) > 1){ # & !is.na(panel.dim)
makefile <- system.file("grass-image_makefiles/Makefile_multi", package = "DataflowR")
system(paste0("make -f ", makefile,
" multipanel.png BASEDIR=", fdir,
" PARAM=", shQuote(legendname),
" LEGENDUNITS=", legendunits_print,
" LEGENDUNITSSPACING=", legendunits_spacing,
" LEGEND_XLIM=", legend_xlim,
" PARAMXCOORD=", paramxcoord,
" NROW=", numrow,
" NCOL=", numcol,
" LEGEND_CROP_EXTENT=", legend_crop_extent
))
if(cleanup == TRUE){
system(paste0("make -f ", makefile, " clean"))
}
}
if(cleanup == TRUE){
rmlist <- list.files(file.path(paste(fdir,"/QGIS_plotting", sep = "")), pattern = paste("out", "*", sep = ""), include.dirs = TRUE, full.names = TRUE)
file.remove(rmlist)
}
}
#'@name average_rlist
#'@title Create an average raster from a raster list
#'@param flist character file.list
#'@param percentcov numeric tolerance to include a stacked pixel based on percent NA
#'@description Create an average raster from a raster list
#'@export
average_rlist <- function(flist, percentcov){
rstack <- raster::stack(flist)
rstack <- raster::reclassify(rstack, c(-Inf, 0, NA))
rmean <- raster::calc(rstack, fun = mean, na.rm = T)
rlen <- raster::calc(!is.na(rstack), fun = sum)
rmean[rlen < (percentcov * length(flist))] <- NA
rmean
#res <- cursurf - rmean
#list(rstack = rstack, rmean = rmean, rlen = rlen, res = res)
}
#'@name cv_rlist
#'@title Create a variability raster from a raster list
#'@param flist character file.list
#'@param percentcov numeric tolerance to include a stacked pixel based on percent NA
#'@description Create a variability raster from a raster list. Divide sd by mean.
#'@export
#'@examples \dontrun{
#' rlist <- create_rlist(good_years$x, "chlext")
#' cv_rlist(rlist$rlist, 0.8)
#'}
cv_rlist <- function(flist, percentcov){
rstack <- raster::stack(flist)
rstack <- raster::reclassify(rstack, c(-Inf, 0, NA))
rmean <- raster::calc(rstack, fun = mean, na.rm = T)
rsd <- raster::calc(rstack, fun = sd, na.rm = T)
rlen <- raster::calc(!is.na(rstack), fun = sum)
rmean[rlen < (percentcov * length(flist))] <- NA
rsd[rlen < (percentcov * length(flist))] <- NA
rsd / rmean
}
# #create florida inset
# library(mapdata)
# data(worldHiresMapEnv)
# map("worldHires","usa",xlim=c(-86,-80),ylim=c(24,31),fill=TRUE,col="gray95")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.