diff --git a/DESCRIPTION b/DESCRIPTION index c2dbf9d5120f8e134ee8b7d3dd91872960c618ff..0b45aff93c7ea89ce8ad247c16d478727ddc7013 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -26,8 +26,6 @@ Depends: Imports: abind, bigmemory, - GEOmap, - geomapdata, graphics, grDevices, mapproj, diff --git a/NAMESPACE b/NAMESPACE index c02f502122458ebeec850a5c93b498ab9de625f7..42ebf32aa133f32d15189d4a9f26abe9df0ba500 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -71,11 +71,9 @@ export(Trend) export(UltimateBrier) export(clim.colors) export(clim.palette) -import(GEOmap) import(NbClust) import(SpecsVerification) import(bigmemory) -import(geomapdata) import(graphics) import(mapproj) import(maps) diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R index f6a85380bedaa0649ae56026c2721a6e308f4f2f..ee9945f3dacf75eb67d237e37cc3f41c8b0a4ac4 100644 --- a/R/PlotEquiMap.R +++ b/R/PlotEquiMap.R @@ -62,14 +62,22 @@ #' Takes the value gray(0.5) by default or, if 'square = FALSE', takes the #' value FALSE. If set to FALSE, continents are not filled in. #'@param filled.oceans A logical value or the color name to fill in drawn -#' projected oceans. The default value is FALSE. If it is TRUE, the default +#' projected oceans. The default value is FALSE. If it is TRUE, the default #' colour is "light blue". +#'@param country.borders A logical value indicating if the country borders +#' should be plotted (TRUE) or not (FALSE). It only works when +#' 'filled.continents' is FALSE. The default value is FALSE. #'@param coast_color Colour of the coast line of the drawn projected continents. #' Takes the value gray(0.5) by default. #'@param coast_width Line width of the coast line of the drawn projected #' continents. Takes the value 1 by default. #'@param lake_color Colour of the lake or other water body inside continents. -#' The default value is NULL. +#' The default value is NULL. +#'@param shapefile A character string of the path to a .rds file or a list +#' object containinig shape file data. If it is a .rds file, it should contain +#' a list. The list should contains 'x' and 'y' at least, which indicate the +#' location of the shape. If it is specified, 'country.borders' won't be +#' used. The default value is NULL. #'@param contours Array of same dimensions as 'var' to be added to the plot #' and displayed with contours. Parameter 'brks2' is required to define the #' magnitude breaks for each contour curve. Disregarded if 'square = FALSE'. @@ -218,7 +226,7 @@ #'PlotEquiMap(sampleData$mod[1, 1, 1, 1, , ], sampleData$lon, sampleData$lat, #' toptitle = 'Predicted sea surface temperature for Nov 1960 from 1st Nov', #' sizetit = 0.5) -#'@import graphics GEOmap geomapdata maps +#'@import graphics maps #'@importFrom grDevices dev.cur dev.new dev.off gray #'@importFrom stats cor #'@export @@ -228,8 +236,9 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, triangle_ends = NULL, col_inf = NULL, col_sup = NULL, colNA = NULL, color_fun = clim.palette(), square = TRUE, filled.continents = NULL, - filled.oceans = FALSE, + filled.oceans = FALSE, country.borders = FALSE, coast_color = NULL, coast_width = 1, lake_color = NULL, + shapefile = NULL, contours = NULL, brks2 = NULL, contour_lwd = 0.5, contour_color = 'black', contour_lty = 1, contour_draw_label = TRUE, contour_label_scale = 1, @@ -434,6 +443,11 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, ocean_color <- "light blue" } + # Check country.borders + if (!is.logical(country.borders)) { + stop("Parameter 'country.borders' must be logical.") + } + # Check coast_color if (is.null(coast_color)) { if (filled.continents) { @@ -453,15 +467,45 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, # Check lake_color if (!is.null(lake_color)) { -# if (filled.continents) { -# lake_color <- 'white' -# } -# } else { if (!.IsColor(lake_color)) { stop("Parameter 'lake_color' must be a valid colour identifier.") } } + # shapefile + if (!is.null(shapefile)) { + if (is.list(shapefile)) { + shape <- shapefile + if (any(!c('x', 'y') %in% names(shape))) { + stop("The list names of the object in 'shapefile' .rds file should ", + "have at least 'x' and 'y'.") + } + if (length(shape$x) != length(shape$y)) { + stop("The length of x and y in 'shapefile' list should be equal.") + } + } else if (!is.character(shapefile)) { + stop("Parameter 'shapefile' must be a .rds file or a list.") + } else { # .rds file + if (!file.exists(shapefile)) { + stop("Parameter 'shapefile' is not a valid file.") + } + if (!grepl("\\.rds$", shapefile)) { + stop("Parameter 'shapefile' must be a .rds file or a list.") + } + shape <- readRDS(file = shapefile) + if (!is.list(shape)) { + stop("Parameter 'shapefile' should be a .rds file of a list object.") + } + if (any(!c('x', 'y') %in% names(shape))) { + stop("The list names of the object in 'shapefile' .rds file should ", + "have at least 'x' and 'y'.") + } + if (length(shape$x) != length(shape$y)) { + stop("The length of x and y in 'shapefile' list should be equal.") + } + } + } + # Check contours if (!is.null(contours)) { if (dim(contours)[1] != dims[1] || dim(contours)[2] != dims[2]) { @@ -656,10 +700,6 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, } } - #library(GEOmap) - #library(geomapdata) - #library(maps) - utils::data(coastmap, package = 'GEOmap', envir = environment()) # # Input arguments # ~~~~~~~~~~~~~~~~~ @@ -844,47 +884,26 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, } old_lwd <- par('lwd') par(lwd = coast_width) - # If [0, 360], use GEOmap; if [-180, 180], use maps - if (min(lon) >= 0) { - ylat <- latmin:latmax - xlon <- lonmin:lonmax - proj <- GEOmap::setPROJ(1, LON0 = mean(xlon), - LAT0 = mean(ylat), LATS = ylat, LONS = xlon) - lakes <- which(coastmap$STROKES$col == "blue") - par(new = TRUE) - if (filled.continents) { - coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- continent_color - if (is.null(lake_color)) { - coastmap$STROKES$col[lakes] <- continent_color - } else { - coastmap$STROKES$col[lakes] <- lake_color #"white" - } - GEOmap::plotGEOmap(coastmap, PROJ = proj, border = coast_color, - add = TRUE, lwd = coast_width) - } else { - coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- coast_color - if (is.null(lake_color)) { - coastmap$STROKES$col[lakes] <- coast_color - } else { - coastmap$STROKES$col[lakes] <- lake_color #"white" - } - GEOmap::plotGEOmap(coastmap, PROJ = proj, #MAPcol = coast_color, - add = TRUE, lwd = coast_width, MAPstyle = 2) - } - + # If [0, 360], use GEOmap; if [-180, 180], use maps::map + # UPDATE: Use maps::map for both cases. The difference between GEOmap and + # maps is trivial. The only thing we can see for now is that + # GEOmap has better lakes. + if (!is.null(shapefile)) { + coast <- maps::map(shape, interior = country.borders, wrap = TRUE, + xlim = xlim_conti, ylim = c(-89.99, 89.99), + fill = filled.continents, add = TRUE, plot = FALSE) } else { - # [-180, 180] - coast <- map(continents, interior = FALSE, wrap = TRUE, + coast <- maps::map(continents, interior = country.borders, wrap = TRUE, xlim = xlim_conti, ylim = c(-89.99, 89.99), fill = filled.continents, add = TRUE, plot = FALSE) - if (filled.continents) { - polygon(coast, col = continent_color, border = coast_color, lwd = coast_width) - } else { - lines(coast, col = coast_color, lwd = coast_width) - } - if (!is.null(lake_color)) { - map('lakes', add = TRUE, fill = filled.continents, col = lake_color) - } + } + if (filled.continents) { + polygon(coast, col = continent_color, border = coast_color, lwd = coast_width) + } else { + lines(coast, col = coast_color, lwd = coast_width) + } + if (!is.null(lake_color)) { + maps::map('lakes', add = TRUE, fill = filled.continents, col = lake_color) } par(lwd = old_lwd) @@ -893,7 +912,7 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, old_lwd <- par('lwd') par(lwd = coast_width) - outline <- map(continents, fill = T, plot = FALSE) # must be fill = T + outline <- maps::map(continents, fill = T, plot = FALSE) # must be fill = T xbox <- xlim_conti + c(-2, 2) ybox <- c(-92, 92) outline$x <- c(outline$x, NA, c(xbox, rev(xbox), xbox[1])) diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd index b574904cedbc13e914b8059e21ef3e8e5119b7a1..52637102ce25fbe027df54bd1d3da9fa77b6cc4e 100644 --- a/man/PlotEquiMap.Rd +++ b/man/PlotEquiMap.Rd @@ -24,9 +24,11 @@ PlotEquiMap( square = TRUE, filled.continents = NULL, filled.oceans = FALSE, + country.borders = FALSE, coast_color = NULL, coast_width = 1, lake_color = NULL, + shapefile = NULL, contours = NULL, brks2 = NULL, contour_lwd = 0.5, @@ -143,9 +145,13 @@ Takes the value gray(0.5) by default or, if 'square = FALSE', takes the value FALSE. If set to FALSE, continents are not filled in.} \item{filled.oceans}{A logical value or the color name to fill in drawn -projected oceans. The default value is FALSE. If it is TRUE, the default +projected oceans. The default value is FALSE. If it is TRUE, the default colour is "light blue".} +\item{country.borders}{A logical value indicating if the country borders +should be plotted (TRUE) or not (FALSE). It only works when +'filled.continents' is FALSE. The default value is FALSE.} + \item{coast_color}{Colour of the coast line of the drawn projected continents. Takes the value gray(0.5) by default.} @@ -155,6 +161,12 @@ continents. Takes the value 1 by default.} \item{lake_color}{Colour of the lake or other water body inside continents. The default value is NULL.} +\item{shapefile}{A character string of the path to a .rds file or a list +object containinig shape file data. If it is a .rds file, it should contain +a list. The list should contains 'x' and 'y' at least, which indicate the +location of the shape. If it is specified, 'country.borders' won't be +used. The default value is NULL.} + \item{contours}{Array of same dimensions as 'var' to be added to the plot and displayed with contours. Parameter 'brks2' is required to define the magnitude breaks for each contour curve. Disregarded if 'square = FALSE'.} diff --git a/tests/testdata/testIrregularGrid_1.rds b/tests/testdata/testIrregularGrid_1.rds new file mode 100644 index 0000000000000000000000000000000000000000..21a46a726737566eee7880c86aacb45d7fa37ba2 Binary files /dev/null and b/tests/testdata/testIrregularGrid_1.rds differ diff --git a/tests/testdata/testIrregularGrid_2.rds b/tests/testdata/testIrregularGrid_2.rds new file mode 100644 index 0000000000000000000000000000000000000000..d014eafdf520077459c3bd6e49948ed7608c2eba Binary files /dev/null and b/tests/testdata/testIrregularGrid_2.rds differ diff --git a/tests/testdata/testIrregularGrid_3.rds b/tests/testdata/testIrregularGrid_3.rds new file mode 100644 index 0000000000000000000000000000000000000000..3dc6acb70c4f69ece6d812992335391289a1c3f7 Binary files /dev/null and b/tests/testdata/testIrregularGrid_3.rds differ diff --git a/tests/testthat/test-CDORemap.R b/tests/testthat/test-CDORemap.R index 4d70ca7f39eed917fe0f3ade387116a4e759d803..16da1b090a39a2dfa5659cfabb5c8c2af57aa62e 100644 --- a/tests/testthat/test-CDORemap.R +++ b/tests/testthat/test-CDORemap.R @@ -4,23 +4,12 @@ context("s2dv::CDORemap tests") data1 <- array(1:360*181*2, dim = c(lon = 360, lat = 181, time = 2)) lons1 <- seq(0, 359) lats1 <- seq(-90, 90) -data1_1 <- s2dv:::.aperm2(data1, c(3, 1, 2)) -data1_2 <- s2dv::InsertDim(data1, 1, 1, name = 'var') -data1_3 <- s2dv::InsertDim(data1_2, 5, 1, name = 'dat') +data1_1 <- .aperm2(data1, c(3, 1, 2)) +data1_2 <- InsertDim(data1, 1, 1, name = 'var') +data1_3 <- InsertDim(data1_2, 5, 1, name = 'dat') # data2: irregular grid -path2 <- '/esarchive/exp/ipsl-cm6a-lr/cmip6-dcppA-hindcast_i1p1/original_files/cmorfiles/DCPP/IPSL/IPSL-CM6A-LR/dcppA-hindcast/r1i1p1f1/SImon/siconc/gn/v20200101/$var$_SImon_IPSL-CM6A-LR_dcppA-hindcast_s2016-r1i1p1f1_gn_201701-202612.nc' -library(startR) -suppressWarnings( -data2 <- Start(dat = path2, - var = 'siconc', - x = 'all', #indices(2:361), - y = 'all', #indices(2:331), - time = indices(1), #'all', - return_vars = list(#x = NULL, y = NULL, - nav_lat = NULL, nav_lon = NULL), - retrieve = TRUE, silent = T) -) +data2 <- readRDS('../testdata/testIrregularGrid_1.rds') lons2 <- attributes(data2)$Variables$common$nav_lon#[2:361, 2:331] lats2 <- attributes(data2)$Variables$common$nav_lat#[2:361, 2:331] @@ -28,57 +17,27 @@ data2_1 <- ClimProjDiags::Subset(data2, 1, 1, drop = 'selected') data2_2 <- ClimProjDiags::Subset(data2, 5, 1, drop = 'selected') data2_3 <- ClimProjDiags::Subset(data2, c(1, 2), list(1, 1), drop = 'selected') data2_4 <- ClimProjDiags::Subset(data2, c(1, 2, 5), list(1, 1, 1), drop = 'selected') -data2_5 <- s2dv:::.aperm2(data2, c(3,4,1,2,5)) -data2_6 <- s2dv:::.aperm2(data2, c(3,4,5,2,1)) -data2_7 <- s2dv:::.aperm2(data2, c(4,3,1,2,5)) -data2_8 <- s2dv:::.aperm2(data2_4, c(2,1)) +data2_5 <- .aperm2(data2, c(3,4,1,2,5)) +data2_6 <- .aperm2(data2, c(3,4,5,2,1)) +data2_7 <- .aperm2(data2, c(4,3,1,2,5)) +data2_8 <- .aperm2(data2_4, c(2,1)) # data3: irregular grid -path3 <- '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/r1i1p1f1/Omon/$var$/gn/v20190713/$var$_*_s$sdate$-$member$_gn_$aux$.nc' -suppressWarnings( -data3 <- Start(dataset = path3, - var = 'tos', - sdate = paste0(1960), - aux = indices(1), - aux_depends = 'sdate', - j = 'all', - i = 'all', - time = indices(1), - member = 'r1i1p1f1', - return_vars = list(j = NULL, i = NULL, - latitude = NULL, longitude = NULL), - retrieve = T, silent = T) -) +data3 <- readRDS('../testdata/testIrregularGrid_2.rds') lons3 <- attributes(data3)$Variables$common$longitude lats3 <- attributes(data3)$Variables$common$latitude data3_1 <- drop(data3) -data3_2 <- s2dv:::.aperm2(data3_1, c(2, 1)) +data3_2 <- .aperm2(data3_1, c(2, 1)) # data4: irregular grid -path4 <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/cmcc-cm2-sr5/cmip6-dcppA-hindcast_i1p1/', - 'DCPP/CMCC/CMCC-CM2-SR5/dcppA-hindcast/$member$/Omon/$var$/gn/v20210312/', - '$var$_*_s$sdate$-$member$_gn_$aux$.nc') -suppressWarnings( -data4 <- Start(dataset = path4, - var = 'tos', - sdate = '1960', - aux = 'all', - aux_depends = 'sdate', - j = 'all', - i = 'all', - time = indices(1), - member = 'r1i1p1f1', - return_vars = list(j = NULL, i = NULL, - latitude = NULL, longitude = NULL), - retrieve = T, silent = T) -) +data4 <- readRDS('../testdata/testIrregularGrid_3.rds') lons4 <- attributes(data4)$Variables$common$longitude lats4 <- attributes(data4)$Variables$common$latitude data4_1 <- drop(data4) data4_2 <- ClimProjDiags::Subset(data4, c(1,2,3,4,8), list(1,1,1,1,1), drop = 'selected') -data4_3 <- s2dv:::.aperm2(data4_2, c(1, 3, 2)) +data4_3 <- .aperm2(data4_2, c(1, 3, 2)) ############################################## @@ -115,7 +74,7 @@ res3 <- CDORemap(data1_3, lons1, lats1, grid = 'r100x50', method = 'bil', crop = expect_equal( as.vector(res), -as.vector(s2dv:::.aperm2(res1, c(2,3,1))) +as.vector(.aperm2(res1, c(2,3,1))) ) expect_equal( as.vector(res),