#' Get spatial map of the input dataset.
#'
#' @param dataset A list containing different information, should be the result of reading netcdf file using
#' \code{loadNcdf}.
#' @param method A string showing different calculating method for the map. More information please refer to
#' details.
#' @param member A number showing which member is selected to get, if the dataset has a "member" dimension. Default
#' is NULL, if no member assigned, and there is a "member" in dimensions, the mean value of the members will be
#' taken.
#' @param ... several arguments including x, y, title, catchment, point, output, name, info, scale, color,
#' type in \code{?getSpatialMap_mat} for details.
#' @return A matrix representing the raster map is returned, and the map is plotted.
#' @details
#' There are following methods to be selected,
#' "meanAnnual": annual rainfall of each year is plotted.
#' "winter", "spring", "autumn", "summer": MEAN seasonal rainfall of each year is plotted.
#' Month(number 1 to 12): MEAN month rainfall of each year is plotted, e.g. MEAN march rainfall of each year.
#' "mean", "max", "min": mean daily, maximum daily, minimum daily precipitation.
#' @examples
#'
#'
#' \dontrun{
#' #gridData provided in the package is the result of \code {loadNcdf}
#' data(tgridData)
#' getSpatialMap(tgridData, method = 'meanAnnual')
#' getSpatialMap(tgridData, method = 'winter')
#'
#'
#' getSpatialMap(tgridData, method = 'winter', catchment = testCat)
#'
#' file <- system.file("extdata", "point.txt", package = "hyfo")
#' point <- read.table(file, header = TRUE, sep = ',' )
#' getSpatialMap(tgridData, method = 'winter', catchment = testCat, point = point)
#' }
#'
#'
#' # More examples can be found in the user manual on http://yuanchao-xu.github.io/hyfo/
#'
#' @export
getSpatialMap <- function(dataset, method = NULL, member = 'mean', ...) {
#check input dataset
checkHyfo(dataset)
#range of the dataset just loaded
lon <- dataset$xyCoords$x
lat <- dataset$xyCoords$y
startTime <- as.POSIXlt(dataset$Dates$start, tz = 'GMT')
yearIndex <- startTime$year + 1900
monthIndex <-startTime$mon + 1
data <- dataset$Data
# Dimension needs to be arranged. Make sure first and second dimension is lat and lon.
data <- adjustDim(data, ref = c('lon', 'lat', 'time'))
# Because in the following part, only 3 dimensions are allowed, so data has to be processed.
if (member == 'mean' & any(attributes(data)$dimensions == 'member')) {
dimIndex3 <- which(attributes(data)$dimensions != 'member')
data <- apply(data, MARGIN = dimIndex3, FUN = mean, na.rm = TRUE)
message('Mean value of the members are returned.')
} else if (member != 'mean' & any(attributes(data)$dimensions == 'member')) {
dimIndex3 <- which(attributes(data)$dimensions == 'member')
data <- chooseDim(data, dimIndex3, member, drop = TRUE)
} else if (member != 'mean' & !any(attributes(data)$dimensions == 'member')){
stop('There is no member part in the dataset, but you choose one, check the input
dataset or change your arguments.')
}
if (is.null(method)) {
warning('You should shoose a method, unless input is a matrix directly to be plotted.')
#in case the dataset is ready to plot and no need to calculate
} else if (method == 'meanAnnual') {
#mean value of the annual precipitation over the period of the data
#time <- proc.time()
if (length(unique(monthIndex)) < 12) {
warning ('There are less than 12 months in a year, the results may be inaccurate.')
}
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, method = 'annual')
#newTime <- proc.time() - time
title_d <- 'Mean Annual Precipitation (mm / year)'
} else if (method == 'winter') {
#mean value of the seasonal precipitation, in this case, winter
#time <- proc.time()
wm <- match(c(12, 1, 2), unique(monthIndex))
if (length(which(!is.na(wm))) < 3) {
stop ('Winter has less than 3 months, check data and try to calculate every month
seperately or choose another season.')
}
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, monthIndex = monthIndex,
method = 'winter')
#newTime <- proc.time() - time
title_d <- 'Mean Winter Precipitation (mm / winter)'
} else if (method == 'spring') {
wm <- match(c(3, 4, 5), unique(monthIndex))
if (length(which(!is.na(wm))) < 3) {
stop ('Spring has less than 3 months, check data and try to calculate every month
seperately or choose another season.')
}
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, monthIndex = monthIndex,
method = 'spring')
title_d <- 'Mean Spring Precipitation (mm / spring)'
} else if (method == 'summer') {
wm <- match(c(6, 7, 8), unique(monthIndex))
if (length(which(!is.na(wm))) < 3) {
stop ('Summer has less than 3 months, check data and try to calculate every month
seperately or choose another season.')
}
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, monthIndex = monthIndex,
method = 'summer')
title_d <- 'Mean Summer Precipitation (mm / summer)'
} else if (method == 'autumn') {
wm <- match(c(9, 10, 11), unique(monthIndex))
if (length(which(!is.na(wm))) < 3) {
stop ('Autumn has less than 3 months, check data and try to calculate every month
seperately or choose another season.')
}
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, monthIndex = monthIndex,
method = 'autumn')
title_d <- 'Mean Autumn Precipitation (mm / autumn)'
} else if (method == 'mean') {
#sum value of the dataset, this procedure is to get the mean value
data_new <- apply(data, MARGIN = c(2, 1), FUN = mean, na.rm = TRUE)
title_d <- 'Mean Daily Precipitation (mm / day)'
} else if (method == 'max') {
data_new <- apply(data, MARGIN = c(2, 1), FUN = suppressWarnings(max), na.rm = TRUE)
title_d <- 'Max Daily Precipitation (mm / day)'
} else if (method == 'min') {
data_new <- apply(data, MARGIN = c(2, 1), FUN = suppressWarnings(min), na.rm = TRUE)
title_d <- 'Min Daily Precipitation (mm / day)'
} else if (is.numeric(method)) {
data_new <- apply(data, MARGIN = c(2, 1), FUN = getMeanPreci, yearIndex = yearIndex, monthIndex = monthIndex,
method = method)
title_d <- paste(month.abb[method], 'Precipitation (mm / month)', sep = ' ')
} else {
wrongMethod <- method
stop(paste('no method called', wrongMethod))
}
# This is to give attributes to the matrix and better be melted in ggplot.
colnames(data_new) <- round(lon, 2)
rownames(data_new) <- round(lat, 2)
# If ... also has a title argument, this will cause conflicts. so title has to be renamed as title_d
# This has to be paid a lot of attention when use ... to pass arguments.
output <- getSpatialMap_mat(matrix = data_new, title_d = title_d, ...)
return(output)
}
#' Replot raster matrix
#'
#' replot the matrix output from \code{getSpatialMap}, when \code{output = 'data'} or output is default
#' value.
#'
#' @param matrix A matrix raster, should be the result of \code{getSpatialMap()}, output should be default
#' or 'data'
#' @param title_d A string showing the title of the plot, defaut is NULL.
#' @param catchment A catchment file geting from \code{shp2cat()} in the package, if a catchment is available for background.
#' @param point A dataframe, showing other information, e.g., location of the gauging stations. The
#' the data.frame should be with columes "name, lon, lat, z, value".
#' @param output A string showing the type of the output, if \code{output = 'ggplot'}, the returned
#' data can be used in ggplot and \code{getSpatialMap_comb()}; if \code{output = 'plot'}, the returned data is the plot containing all
#' layers' information, and can be plot directly or used in grid.arrange; if not set, the raster matrix data
#' will be returned.
#' @param name If \code{output = 'ggplot'}, name has to be assigned to your output, in order to differentiate
#' different outputs in the later multiplot using \code{getSpatialMap_comb}.
#' @param info A boolean showing whether the information of the map, e.g., max, mean ..., default is FALSE.
#' @param scale A string showing the plot scale, 'identity' or 'sqrt'.
#' @param color Most of time you don't have to set this, but if you are not satisfied with the
#' default color, you can set your own palette here. e.g., \code{color = c('red', 'blue')}, then
#' the value from lowest to highest, will have the color from red to blue. More info about color,
#' please check ?palette().
#' @param ... \code{title, x, y} showing the title and x and y axis of the plot. e.g. \code{title = 'aaa'}
#'default is about precipitation.
#' @return A matrix representing the raster map is returned, and the map is plotted.
#' @examples
#'
#' \dontrun{
#' data(tgridData)# the result of \code{loadNcdf}
#' #the output type of has to be default or 'data'.
#' a1 <- getSpatialMap(tgridData, method = 'mean')
#' a2 <- getSpatialMap(tgridData, method = 'max')
#' a3 <- getSpatialMap(tgridData, method = 'winter')
#' a4 <- getSpatialMap(tgridData, method = 'summer')
#' #For example, if we want to investigate the difference between mean value and max.
#'
#' a5 <- a2 - a1
#' getSpatialMap_mat(a4)
#'
#' #Or to investigate the difference between winter value and summer value.
#' a6 <- a3 - a4
#' getSpatialMap_mat(a6)
#'
#' }
#'
#'
#' # More examples can be found in the user manual on https://yuanchao-xu.github.io/hyfo/
#'
#' @export
#' @import ggplot2 plyr maps sp sf
#' @importFrom stats median
#' @importFrom methods is
#' @importFrom reshape2 melt
#' @references
#'
#' \itemize{
#' \item R Core Team (2015). R: A language and environment for statistical computing. R Foundation for
#' Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
#'
#' \item Hadley Wickham (2007). Reshaping Data with the reshape Package. Journal of Statistical Software,
#' 21(12), 1-20. URL http://www.jstatsoft.org/v21/i12/.
#'
#' \item Hadley Wickham (2011). The Split-Apply-Combine Strategy for Data Analysis. Journal of Statistical
#' Software, 40(1), 1-29. URL http://www.jstatsoft.org/v40/i01/.
#'
#' \item Original S code by Richard A. Becker and Allan R. Wilks. R version by Ray Brownrigg. Enhancements
#' by Thomas P Minka <tpminka at media.mit.edu> (2015). maps: Draw Geographical Maps. R package version
#' 2.3-11. https://CRAN.R-project.org/package=maps
#'
#' \item Pebesma, Edzer, and Roger Bivand. 2023a. Sp: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=sp.
#'
#' \item Roger Bivand and Colin Rundel (2015). rgeos: Interface to Geometry Engine - Open Source (GEOS). R
#' package version 0.3-11. https://CRAN.R-project.org/package=sf
#'
#' }
#'
#'
#'
#'
#'
getSpatialMap_mat <- function(matrix, title_d = NULL, catchment = NULL, point = NULL, output = 'data',
name = NULL, info = FALSE, scale = 'identity', color = NULL, ...) {
#check input
checkWord <- c('lon', 'lat', 'z', 'value')
if (is.null(attributes(matrix)$dimnames)) {
stop('Input matrix is incorrect, check help to know how to get the matrix.')
} else if (!is.null(catchment) & is(catchment)[1] != "SpatialPolygonsDataFrame") {
stop('Catchment format is incorrect, check help to get more details. ')
} else if (!is.null(point) & any(is.na(match(checkWord, attributes(point)$names)))) {
stop('point should be a dataframe with colnames "lon, lat, z, value".')
}
#ggplot
#for the aes option in ggplot, it's independent from any other command through all ggplot, and aes() function
#get data from the main dataset, in this case, data_ggplot. for other functions in ggplot, if it wants to use
#data from the main dataset as parameters, it has to use aes() function. if not, it has to use data available
#in the environment.
#in other words, all the parameters in aes(), they have to come from the main dataset. Otherwise, just put them
#outside aes() as normal parameters.
if (info == TRUE) {
plotMax <- round(max(matrix, na.rm = TRUE), 2)
plotMin <- round(min(matrix, na.rm = TRUE), 2)
plotMean <- round(mean(matrix, na.rm = TRUE), 2)
plotMedian <- round(median(matrix, na.rm = TRUE), 2)
word <- paste('\n\n', paste('Max', '=', plotMax), ',', paste('Min', '=', plotMin), ',',
paste('Mean', '=', plotMean), ',', paste('Median', '=', plotMedian))
} else {
word <- NULL
}
x_word <- paste('Longitude', word)
world_map <- map_data('world')
# For some cases, matrix has to be reshaped, because it's too fat or too slim, to make
# it shown on the map, the ratio is x : y is 4 : 3.
matrix <- reshapeMatrix(matrix)
# cannot remove NA, or the matrix shape will be changed.
data_ggplot <- melt(matrix, na.rm = FALSE)
colnames(data_ggplot) <- c('lat', 'lon', 'value')
theme_set(theme_bw())
if (is.null(color)) color <- c('yellow', 'orange', 'red')
# if (is.null(color)) color <- rev(rainbow(n = 20, end = 0.7))
mainLayer <- with(data_ggplot, {
ggplot(data = data_ggplot) +
geom_tile(aes(x = lon, y = lat, fill = value)) +
#scale_fill_discrete()+
scale_fill_gradientn(colours = color, na.value = 'transparent') +#usually scale = 'sqrt'
#guide = guide_colorbar, colorbar and legend are not the same.
guides(fill = guide_colourbar(title='Rainfall (mm)', barheight = rel(9), trans = scale)) +#usually scale = 'sqrt'
geom_map(data = world_map, map = world_map, aes(map_id = region), fill = 'transparent',
color='black') +
# guides(fill = guide_colorbar(title='Rainfall (mm)', barheight = 15))+
xlab(x_word) +
ylab('Latitude') +
ggtitle(title_d) +
labs(empty = NULL, ...) +#in order to pass "...", arguments shouldn't be empty.
theme(plot.title = element_text(size = rel(1.8), face = 'bold'),
axis.title.x = element_text(size = rel(1.7)),
axis.title.y = element_text(size = rel(1.7)),
axis.text.x = element_text(size = rel(1.9)),
axis.text.y = element_text(size = rel(1.9)),
legend.text = element_text(size = rel(1.3)),
legend.title = element_text(size = rel(1.3)))
# coord_fixed(ratio = 1, xlim = xlim, ylim = ylim)
# geom_rect(xmin=min(lon)+0.72*(max(lon)-min(lon)),
# xmax=min(lon)+0.99*(max(lon)-min(lon)),
# ymin=min(lat)+0.02*(max(lat)-min(lat)),
# ymax=min(lat)+0.28*(max(lat)-min(lat)),
# fill='white',colour='black')+
# annotate('text', x = min(lon), y = min(lat), label=word, hjust = 0, vjust = -1)
})
printLayer <- mainLayer
#catchment conversion
if (is.null(catchment) == FALSE) {
a <- catchment
a@data$id <- rownames(a@data)
b <- fortify(a, region = 'id')
c <- join(b, a@data, by = 'id')
catchmentLayer <- with(c, {
geom_polygon(data = c, aes(long, lat, group = group), color = 'black',
fill = 'transparent')
})
printLayer <- printLayer + catchmentLayer
}
#plot point
if (is.null(point) == FALSE) {
pointLayer <- with(point, {
geom_point(data = point, aes(x = lon, y = lat, size = value, colour = z),
guide = guide_legend(barheight = rel(3)))
})
printLayer <- printLayer + pointLayer
}
print(printLayer)
if (output == 'ggplot') {
if (is.null(name)) stop('"name" argument not found,
If you choose "ggplot" as output, please assign a name.')
data_ggplot$Name <- rep(name, dim(data_ggplot)[1])
return (data_ggplot)
} else if (output == 'plot') {
return(printLayer)
} else {
return(matrix)
}
}
#' Combine maps together
#' @param ... different maps generated by \code{getSpatialMap(, output = 'ggplot')}, see details.
#' @param nrow A number showing the number of rows.
#' @param list If input is a list containing different ggplot data, use \code{list = inputlist}.
#' @param x A string of x axis name.
#' @param y A string of y axis name.
#' @param title A string of the title.
#' @param output A boolean, if chosen TRUE, the output will be given.
#' @return A combined map.
#' @examples
#'
#'
#' \dontrun{
#' data(tgridData)# the result of \code{\link{loadNcdf}}
#' #The output should be 'ggplot'
#' a1 <- getSpatialMap(tgridData, method = 'summer', output = 'ggplot', name = 'a1')
#' a2 <- getSpatialMap(tgridData, method = 'winter', output = 'ggplot', name = 'a2')
#' a3 <- getSpatialMap(tgridData, method = 'mean', output = 'ggplot', name = 'a3')
#' a4 <- getSpatialMap(tgridData, method = 'max', output = 'ggplot', name = 'a4')
#' getSpatialMap_comb(a1, a2)
#'
#' # or you can put them into a list.
#' getSpatialMap_comb(list = list(a1, a2), nrow = 2)
#' }
#'
#'
#'
#' # More examples can be found in the user manual on https://yuanchao-xu.github.io/hyfo/
#'
#' @details
#' For \code{getSpatialMap_comb}, the maps to be compared should be with same size and resolution,
#' in other words, they should be fully overlapped by each other.
#'
#' If they have different resolutions, use \code{interpGridData{ecomsUDG.Raccess}} to interpolate.
#'
#' @export
#' @import ggplot2 maps
#' @importFrom data.table rbindlist
#' @importFrom methods is
#' @references
#'
#' \itemize{
#' \item H. Wickham. ggplot2: elegant graphics for data analysis. Springer New York, 2009.
#' }
getSpatialMap_comb <- function(..., list = NULL, nrow = 1, x = '', y = '', title = '',
output = FALSE) {
if (!is.null(list)) {
data_ggplot <- rbindlist(list)
} else {
maps <- list(...)
checkBind(maps, 'rbind')
data_ggplot <- rbindlist(maps)
}
if (!is(data_ggplot)[1] == 'data.frame') {
warning('Your input is probably a list, but you forget to add "list = " before it.
Try again, or check help for more information.')
} else if (is.null(data_ggplot$Name)) {
stop('No "Name" column in the input data, check the arguments in getSpatialMap(), if
output = "ggplot" is assigned, more info please check ?getSpatialMap().')
}
data_ggplot$Name <- factor(data_ggplot$Name, levels = unique(data_ggplot$Name), ordered = TRUE)
# lim <- getLim(data_ggplot$lon, data_ggplot$lat)
# xlim <- lim[[1]]
# ylim <- lim[[2]]
world_map <- map_data('world')
theme_set(theme_bw())
mainLayer <- with(data_ggplot, {
ggplot(data = data_ggplot) +
geom_tile(aes(x = lon, y = lat, fill = value)) +
#scale_fill_gradient(high = 'red', low = 'yellow')+
scale_fill_gradientn(colours = c('yellow', 'orange', 'red'), na.value = 'transparent') +#usually scale = 'sqrt'
geom_map(data = world_map, map = world_map, aes(map_id = region), fill = 'transparent', color = 'black') +
# guides(fill = guide_colourbar(title='Rainfall (mm)', barheight = rel(9), trans = scale)) +#
facet_wrap(~ Name, nrow = nrow) +
theme(plot.title = element_text(size = rel(1.8), face = 'bold'),
axis.title.x = element_text(size = rel(1.7)),
axis.title.y = element_text(size = rel(1.7)),
axis.text.x = element_text(size = rel(1.9)),
axis.text.y = element_text(size = rel(1.9)),
legend.text = element_text(size = rel(1.3)),
legend.title = element_text(size = rel(1.3))) +
# no solultion for some very fat or very slim, in facet ggplot2, so, it's not buitiful.
#coord_equal() +
labs(x = x, y = y, title = title)
})
suppressWarnings(print(mainLayer))
if (output == TRUE) return(data_ggplot)
}
reshapeMatrix <- function(matrix) {
# This is for the map plot to keep the ratio x : y == 4:3
# mainly used in map plot in ggplot2.
# So the input matrix should be reshaped, add in some NA values,
# in order to be shown appropriately in ggplot.
lon <- as.numeric(colnames(matrix))
lat <- as.numeric(rownames(matrix))
dx <- mean(diff(lon))
dy <- mean(diff(lat))
lx <- max(lon) - min(lon)
ly <- max(lat) - min(lat)
if (0.75 * lx < ly) {
# In this case, x needs to be made longer
xhalf <- 0.67 * ly
xadd <- xhalf - lx / 2
# calculate how many columns needs to be added.
nxadd <- abs(round(xadd / dx))
madd1 <- matrix(data = NA, nrow = length(lat), ncol = nxadd)
madd2 <- madd1
colnames(madd1) <- seq(to = min(lon) - dx, length = nxadd, by = dx)
colnames(madd2) <- seq(from = max(lon) + dx, length = nxadd, by = dx)
matrix_new <- cbind(madd1, matrix, madd2)
} else if (0.75 * lx > ly) {
yhalf <- 0.38 * lx
yadd <- yhalf - ly / 2
nyadd <- abs(round(yadd / dy))
madd1 <- matrix(data = NA, nrow = nyadd, ncol = length(lon))
madd2 <- madd1
rownames(madd1) <- seq(to = max(lat) + dy, length = nyadd, by = -dy)
rownames(madd2) <- seq(from = min(lat) - dx, length = nyadd, by = -dy)
matrix_new <- rbind(madd1, matrix, madd2)
} else {
matrix_new <- matrix
}
return(matrix_new)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.