changes in spatial

Dropping sp and rgdal. Changes to plot_trend: only stamen maps can be easily obtained and are now employed.

.gitignore

@@ -18,3 +18,5 @@ README_cache/*
 inst/doc
 doc
 Meta
+/doc/
+/Meta/

DESCRIPTION

@@ -23,7 +23,7 @@ Authors@R: c(person("Ilaria", "Prosdocimi",
                     introducing vectorisation")
                     )
 Maintainer: Ilaria Prosdocimi <prosdocimi.ilaria@gmail.com>
-URL: http://ilapros.github.io/rnrfa/
+URL: https://ilapros.github.io/rnrfa/
 BugReports: https://github.com/ilapros/rnrfa/issues
 Description: Utility functions to retrieve data from the UK National River Flow
   Archive (<https://nrfa.ceh.ac.uk/>, terms and conditions:
@@ -36,7 +36,6 @@ Description: Utility functions to retrieve data from the UK National River Flow
   <https://journal.r-project.org/archive/2016/RJ-2016-036/RJ-2016-036.pdf>.
 Depends: R (>= 3.5)
 Imports:
-    rgdal,
     curl,
     jsonlite,
     lubridate,
@@ -46,7 +45,7 @@ Imports:
     zoo,
     ggmap,
     ggplot2,
-    sp,
+    sf,
     parallel,
     tibble
 Suggests: 

NAMESPACE

@@ -10,7 +10,6 @@ export(plot_rain_flow)
 export(plot_trend)
 export(seasonal_averages)
 export(station_ids)
-import(rgdal)
 importFrom(curl,has_internet)
 importFrom(ggmap,get_map)
 importFrom(ggmap,ggmap)
@@ -38,10 +37,12 @@ importFrom(httr,http_error)
 importFrom(jsonlite,fromJSON)
 importFrom(lubridate,year)
 importFrom(parallel,parLapply)
-importFrom(sp,CRS)
-importFrom(sp,coordinates)
-importFrom(sp,proj4string)
-importFrom(sp,spTransform)
+importFrom(sf,st_coordinates)
+importFrom(sf,st_crs)
+importFrom(sf,st_multipoint)
+importFrom(sf,st_point)
+importFrom(sf,st_sfc)
+importFrom(sf,st_transform)
 importFrom(stats,aggregate)
 importFrom(stats,glm)
 importFrom(stats,quantile)

NEWS.md

@@ -1,11 +1,14 @@
-
 v2.0.5
 --------------------------------------
-Minor change:
+Major changes: 
+
+- The hard dependency on the rgdal package is dropped (this is in light of the announced retirement of rgdal in 2023). In the process the hard dependence to sp has been dropped in favour of sf, to avoid possible issues down the line when sp goes through the changes necessary to accommodate the rgdal retirement. The package still depends on sp via ggmap -> RgoogleMaps. 
+
+Minor changes:
 
 - To further ensure that the package fails gracefully when the NRFA api is interrogated about stations which do not exist an informative message is printed and the function(s) return a NULL object.  
 - Thanks to Andrew Duncan (@aj2duncan) osg_parse now outputs a missing value if a missing coordinate is given 
-- small fixes in documentation 
+- small fixes in documentation. 
 
 v2.0.4 and submitted to CRAN.
 --------------------------------------

R/internals.R

@@ -32,7 +32,7 @@ nrfa_api <- function(webservice, parameters, path = "") {
 
   # Parse content
   page_content <- try(httr::content(resp, "text", encoding = "UTF-8"))
-  if (class(page_content) == "try-error") {
+  if (inherits(page_content,"try-error")) {
     errs <- geterrmessage()
     message(paste("An unknwon error occurred when accessing the data",
                   "- with error message:", errs))

R/osg_parse.R

@@ -1,6 +1,6 @@
 #' Converts OS Grid Reference to BNG/WGS coordinates.
 #'
-#' @author Claudia Vitolo
+#' @author Claudia Vitolo (Ilaria Prosdocimi ported to sf)
 #'
 #' @description This function converts an Ordnance Survey (OS) grid reference to
 #' easting/northing or latitude/longitude coordinates.
@@ -99,9 +99,9 @@ osg_parse <- function(grid_refs, coord_system = c("BNG", "WGS84")) {
 
 .transform_crs <- function(x, y, from, to) {
 
-    df <- data.frame(x = as.numeric(x), y = as.numeric(y), from, to)
+    df <- data.frame(x = as.numeric(x), y = as.numeric(y), from = from, to = to)
 
-    .transform <- function(x) {
+        .transform <- function(x) {
         # transformation can only be vectorized for unique CRS
         if (length(unique(x$from)) > 1) {
           stop("Cannot handle multiple source CRS.")
@@ -111,22 +111,32 @@ osg_parse <- function(grid_refs, coord_system = c("BNG", "WGS84")) {
         }
 
         xy <- x[, c("x", "y")]
-
+       
         from <- x$from[1]
         to <- x$to[1]
+        rn <- rownames(x)
 
         # nothing to do ...
         if (from == to) return(xy)
-
-        sp::coordinates(xy) <- ~x + y
-        sp::proj4string(xy) <- sp::CRS(paste0("+init=epsg:", from))
-
-        xy_new <- sp::spTransform(xy, sp::CRS(paste0("+init=epsg:", to)))
-
-        as.data.frame(sp::coordinates(xy_new))
+        # 
+        # drop dependency on sp, include dependency on sf 
+        # this is due to the retiremnt of some packages
+        # https://r-spatial.org/r/2022/04/12/evolution.html#packages-depending-on-sp-and-raster
+        # sp::coordinates(xy) <- ~x + y
+        # sp::proj4string(xy) <- sp::CRS(paste0("+init=epsg:", from))
+        # 
+        # xy_new <- sp::spTransform(xy, sp::CRS(paste0("+init=epsg:", to)))
+        # 
+        # as.data.frame(sp::coordinates(xy_new))
+        if(nrow(xy) == 1) pointxy <- sf::st_sfc(sf::st_point(as.matrix(xy)), crs = from)
+        if(nrow(xy) > 1)  pointxy <- sf::st_sfc(sf::st_multipoint(as.matrix(xy)), crs = from)
+        xy_new <- sf::st_transform(pointxy,crs = sf::st_crs(to))
+        xy_new <- sf::st_coordinates(xy_new)[,c("X","Y"), drop = FALSE]
+        colnames(xy_new) <- c("x", "y"); rownames(xy_new) <- rn
+        as.data.frame(xy_new)
     }
 
-    # split to obain unique CRS
+    # split to obtain unique CRS
     grouped <- split(df, f = df[, c("from", "to")])
 
     unsplit(lapply(grouped, .transform), f = df[, c("from", "to")])

R/plot_trend.R

@@ -6,21 +6,38 @@
 #' lat (latitude), lon (longitude), slope and an additional user-defined column
 #' \code{column_name}.
 #' @param column_name name of the column to use for grouping the results.
+#' @param maptype maptype. The default is "toner-lite", see ?ggmap::get_stamenmap for other options 
+#' @param showmap set to FALSE to avoid plotting the map when running the function
 #'
-#' @return Two plots, side-by-side, the first showing the distribution of the
+#' @return Two plots, the first showing the distribution of the
 #' trend over a map, based on the slope of the linear model that describes the
 #' trend. The second plot shows a boxplot of the slope grouped based on the
-#' column Region. Region and slope can be user-defined.
+#' column \code{column_name} and slope can be user-defined 
+#' (notice that in the plot the very extreme slope values are not displayed to avoid skewed visualisations). 
 #'
 #' @export
 #'
 #' @examples
 #' \dontrun{
-#'   plot_trend(df, Region)
+#'   # some fake data around London 
+#'   df <- data.frame(lat = 51.5+runif(40,-0.3,0.3), 
+#'                    lon = 0+runif(40, -0.3,0.3),
+#'                    slope = rnorm(40, c(rep(-0.4,20),rep(0.4,20))), 
+#'                    g = factor(c(rep("a",20), rep("b",20))))
+#'   theplots <- plot_trend(df, "g", maptype = "terrain-background")
+#'   theplots$A # map
+#'   theplots$B + labs(subtitle = "Use ggplot usual commands to modify the plots") # boxplots
+#'   
 #' }
 
-plot_trend <- function(df, column_name) {
-
+plot_trend <- function(df, column_name, maptype = "toner-lite", showmap = TRUE) {
+  names(df) <- tolower(names(df))
+  if(!all(c("lat", "lon") %in% names(df))){
+    if(all(c("latitude", "longitude") %in% names(df))) {
+      df$lat <- df$latitude; df$lon <- df$longitude
+    }
+    else stop("need lat lon columns to proceed")
+  } 
   df$trend <- NA
   df$trend[df$slope >= 0]  <- "Positive"
   df$trend[df$slope < 0]  <- "Negative"
@@ -30,11 +47,30 @@ plot_trend <- function(df, column_name) {
 
   # Plot red dot for decreasing trend, blue dot for increasing trend
   tolerance <- (max(df$lon, na.rm = TRUE) - min(df$lon, na.rm = TRUE)) / 10
-  m <- ggmap::get_map(location = c(min(df$lon, na.rm = TRUE) - 2 * tolerance,
-                            min(df$lat, na.rm = TRUE) - 2 * tolerance,
-                            max(df$lon, na.rm = TRUE) + tolerance,
-                            max(df$lat, na.rm = TRUE)) + tolerance,
-               maptype = "toner-lite")
+  # m <- ggmap::get_map(location = c(min(df$lon, na.rm = TRUE) - 2 * tolerance,
+  #                           min(df$lat, na.rm = TRUE) - 2 * tolerance,
+  #                           max(df$lon, na.rm = TRUE) + tolerance,
+  #                           max(df$lat, na.rm = TRUE)) + tolerance,
+  #              maptype = "toner-lite")
+  ## ggmap only work with stamen, updating so one can choose maptype 
+
+  ## from ggmap 
+    location <- c(left = min(df$lon, na.rm = TRUE) - 2 * tolerance,
+                    bottom = min(df$lat, na.rm = TRUE) - 2 * tolerance,
+                    right= max(df$lon, na.rm = TRUE) + tolerance,
+                    top = max(df$lat, na.rm = TRUE)) + tolerance
+    lon_range <- location[c("left","right")]
+    lat_range <- location[c("bottom","top")]
+
+      # compute zoom
+    lonlength <- diff(lon_range)
+    latlength <- diff(lat_range)
+    zoomlon <- ceiling( log2( 360*2 / lonlength) )
+    zoomlat <- ceiling( log2( 180*2 / latlength) )
+    zoom <- max(zoomlon, zoomlat)
+
+  m <- ggmap::get_stamenmap(bbox = location,
+               maptype = maptype, zoom = zoom)
 
   # Plot map
   plot1 <- ggmap::ggmap(m, alpha = 0.5) +
@@ -46,7 +82,7 @@ plot_trend <- function(df, column_name) {
     ggplot2::theme(legend.position = "top") +
     ggplot2::ggtitle("A")
 
-  # Boxplot by NUTS1 region
+  # Boxplot by grouping variable (e.g. NUTS1 region)
   plot2 <- ggplot2::ggplot(df,
                            ggplot2::aes(x = eval(parse(text = column_name)),
                                         y = slope,
@@ -64,7 +100,7 @@ plot_trend <- function(df, column_name) {
                                                                           0,
                                                                           0))) +
     ggplot2::ggtitle("B")
-
+  if(showmap) print(plot1) 
   return(list("A" = plot1, "B" = plot2))
 
 }

R/rnrfa-package.R

@@ -9,7 +9,7 @@
 #' @docType package
 #' @title UK National River Flow Archive data from R
 #'
-#' @import rgdal
+# @import rgdal
 #' @importFrom curl has_internet
 #' @importFrom ggmap ggmap get_map
 #' @importFrom ggplot2 ggplot geom_point aes coord_flip scale_color_manual theme
@@ -20,7 +20,7 @@
 #' @importFrom jsonlite fromJSON
 #' @importFrom lubridate year
 #' @importFrom parallel parLapply
-#' @importFrom sp coordinates proj4string CRS spTransform
+#' @importFrom sf st_sfc st_point st_multipoint st_crs st_transform st_coordinates
 #' @importFrom stats aggregate glm quantile
 #' @importFrom tibble as_tibble
 #' @importFrom utils str