Merge pull request #19 from harphub/develop

Merge develop

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: harpCore
 Title: Core functions and methods for the harp ecosystem
-Version: 0.2.1
+Version: 0.2.2
 Authors@R: 
     person("Andrew", "Singleton", , "[email protected]", role = c("aut", "cre"))
 Description: harp is a collection of packages for reading, analysing and 
@@ -18,6 +18,7 @@ Imports:
     cli,
     dplyr (>= 1.1),
     glue,
+    lifecycle,
     magrittr,
     matrixStats,
     meteogrid,

NAMESPACE

@@ -61,6 +61,8 @@ S3method(geo_points,harp_grid_df)
 S3method(geo_regrid,geofield)
 S3method(geo_regrid,harp_geolist)
 S3method(geo_regrid,harp_grid_df)
+S3method(geo_reproject,data.frame)
+S3method(geo_reproject,harp_df)
 S3method(geo_subgrid,geodomain)
 S3method(geo_subgrid,geofield)
 S3method(geo_subgrid,harp_geolist)
@@ -206,6 +208,7 @@ export(geo_opts_xsection)
 export(geo_opts_zoom)
 export(geo_points)
 export(geo_regrid)
+export(geo_reproject)
 export(geo_subgrid)
 export(geo_transform)
 export(geo_upscale)

R/geo_transform.R

@@ -1075,6 +1075,168 @@ geo_transform.harp_geolist <- geo_transform_all
 #' @export
 geo_transform.harp_grid_df <- geo_transform_all
 
+
+########################
+# Reproject point data #
+########################
+
+## Helper function to get the projection
+get_projection <- function(x, caller = rlang::caller_env()) {
+
+  if (is_geolist(x) || meteogrid::is.geofield(x)) {
+    return(get_domain(x)[["projection"]])
+  }
+
+  if (meteogrid::is.geodomain(x)) {
+    return(x[["projection"]])
+  }
+
+  if (is.list(x) && !is.null(names(x)) && names(x)[1] == "proj") {
+    return(x)
+  }
+
+  if (is.character(x) && length(x) == 1 && substr(x, 1, 6) == "+proj=") {
+    return(meteogrid::proj4.str2list(x))
+  }
+
+  classes <- cli::cli_inform("")
+  types   <- "a proj string or a meteogrid proj list"
+  cli::cli_abort(c(
+    "Unknown type for {.arg proj}",
+    "i" = "{.arg proj} must be a {.cls geofield}, {.cls geodomain}, {.cls geolist},",
+    "i" = "{types}."
+    ),
+    call = caller
+  )
+}
+
+## Helper function to fix column names
+fix_col_names <- function(col_names, x_col, y_col) {
+  new_names <- gsub("projected_", "", col_names)
+  if (
+    length(which(new_names == x_col)) > 1 ||
+    length(which(new_names == y_col)) > 1
+  ) {
+    return(col_names)
+  }
+  new_names
+}
+
+#' Reproject from or to lat-lon coordinates
+#'
+#'
+#'
+#' @param x A data frame
+#' @param proj The projection. Can be a `geodomain`, a `geofield`, a `geolist`,
+#'   a projection string or a meteogrid projection list. When `inverse = FALSE`,
+#'   this is the projection to which locations in lat-lon coordinates are
+#'   reprojected, and when `inverse = TRUE`, this the projection from which
+#'   locations in projected coordinates are reprojected to lat-lon coordinates.
+#' @param x_col,y_col The names of the columns containing the x and y
+#'   coordinates to be reprojected. For `inverse = FALSE`, these should be
+#'   longitude and latitude in decimal degrees. For `inverse = TRUE`, these
+#'   should be eastings and northings in metres.
+#' @param crop When `proj` is a `geodomain`, `geofield` or `geolist`, set to
+#'   `TRUE` to crop the reprojected locations to only those locations within the
+#'   domain.
+#' @param inverse Set to `TRUE` to reprojected from projected coordinates to
+#'   lat-lon coordinates. The default is `FALSE`
+#'
+#' @return The input data frame with new columns for the reprojected
+#'   coordinates. The projection is added as an attribute. If the data are
+#'   cropped, the domain is also added.
+#' @export
+#'
+#' @examples
+#' geo_reproject(station_list, det_grid_df$fcst)
+#'
+#' # Crop to domain
+#' geo_reproject(station_list, det_grid_df$fcst, crop = TRUE)
+#'
+#' # inverse projection
+#' projected <- geo_reproject(station_list, det_grid_df$fcst)
+#' geo_reproject(
+#'   projected, det_grid_df$fcst, x_col = x, y_col = y, inverse = TRUE
+#' )
+#'
+geo_reproject <- function(
+  x, proj, x_col = "lon", y_col = "lat", crop = FALSE, inverse = FALSE
+) {
+  UseMethod("geo_reproject")
+}
+
+#' @export
+geo_reproject.data.frame <- function(
+    x, proj, x_col = "lon", y_col = "lat", crop = FALSE, inverse = FALSE
+) {
+
+  x_col <- rlang::ensym(x_col)
+  y_col <- rlang::ensym(y_col)
+
+  proj_list <- get_projection(proj)
+
+  x <- dplyr::mutate(
+    x,
+    projected = meteogrid::project(!!x_col, !!y_col, proj_list, inverse)
+  )
+
+  x_col_out <- "projected_x"
+  y_col_out <- "projected_y"
+  if (inverse) {
+    x[["projected"]] <- dplyr::rename(
+      x[["projected"]],
+      lon = dplyr::all_of("x"),
+      lat = dplyr::all_of("y")
+    )
+    x_col_out <- "projected_lat"
+    y_col_out <- "projected_lon"
+  }
+
+  x <- tidyr::unnest(x, dplyr::all_of("projected"), names_sep = "_")
+
+  if (!inverse) {
+    attr(x, "projection") <- meteogrid::proj4.list2str(proj_list)
+  }
+
+  if (!crop) {
+    colnames(x) <- fix_col_names(colnames(x), x_col, y_col)
+    return(x)
+  }
+
+  if (
+    !is_geolist(proj) &&
+    !meteogrid::is.geodomain(proj) &&
+    !meteogrid::is.geofield(proj)
+  ) {
+    cli::cli_warn(c(
+      "For {.arg crop = TRUE}, proj must be a",
+      "{.cls geofield}, {.cls geodomain}, or {.cls geolist}"
+    ))
+
+    colnames(x) <- fix_col_names(colnames(x), x_col, y_col)
+    return(x)
+  }
+
+  dom_ext <- meteogrid::DomainExtent(get_domain(proj))
+
+  x <- dplyr::filter(
+    x,
+    dplyr::between(.data[[x_col_out]], dom_ext[["x0"]], dom_ext[["x1"]]),
+    dplyr::between(.data[[y_col_out]], dom_ext[["y0"]], dom_ext[["y1"]])
+  )
+
+  colnames(x) <- fix_col_names(colnames(x), x_col, y_col)
+  attr(x, "domain") <- get_domain(proj)
+  x
+}
+
+#' @export
+geo_reproject.harp_df <- function(
+    x, proj, x_col = "lon", y_col = "lat", crop = FALSE, inverse = FALSE
+) {
+  as_harp_df(NextMethod())
+}
+
 #################
 # Make a domain #
 #################
@@ -1145,22 +1307,21 @@ tranverse_mercator <- function(ref_lon, ref_lat, tilt) {
 #'
 #' This function is used to define a new domain with a regular grid. At a
 #' minimum, the projection, geographic location of the centre of the domain and
-#' number of horizontal resolution of the grid points must be provided.
+#' number and horizontal resolution of the grid points must be provided.
 #'
-#' @param proj The projection of the domain. See details.
-#' @param centre_lon The longitude of the centre of the domain in decimal
-#'   degrees.
-#' @param centre_lat The latitude of the centre of the domain in decimal
-#'   degrees.
+#' @param proj The projection of the domain. This can be the name of the
+#'   projection or a projection string.
+#' @param centre_lon,centre_lat The longitude and latitude of the centre of the
+#'   domain in decimal degrees.
 #' @param nxny The number of grid points in the x and y directions. Should be a
 #'   vector of length 2 with the number of grid points in the x direction first.
 #'   If only 1 value is given it is assumed to be the same in both directions.
 #' @param dxdy The horizontal resolution of the grid in the x and y directions.
 #'   For lat-lon projections this should be in decimal degrees, otherwise should
 #'   be in metres. Should be a vector of length 2 with the resolution in the x
 #'   direction first.
-#' @param ref_lon The reference longitude of the projection if needed.
-#' @param ref_lat The reference latitude of the projection if needed.
+#' @param ref_lon,ref_lat The reference longitude and latitude of the projection
+#'   if relevant to the projection. Ignored if `proj` is a projection string.
 #' @param exey If defining a grid with an extension zone, a vector length 2 with
 #'   the number of grid points in the x and y directions of the extension zone.
 #'   If only 1 value is given it is assumed to be the same in both directions.
@@ -1178,6 +1339,13 @@ tranverse_mercator <- function(ref_lon, ref_lat, tilt) {
 #'
 #' dd <- define_domain(0, 0, c(360, 180), 1, "latlong") # Whole earth
 #' plot(dd)
+#'
+#' # Pass the projection as a proj string
+#' dd <- define_domain(
+#'   10, 60, 1000, 2500,
+#'   proj = "+proj=lcc +lon_0=15 +lat_0=63.3 +lat_1=63.3 +lat_2=63.3 +R=6371000"
+#' )
+#' plot(dd)
 define_domain <- function(
     centre_lon,
     centre_lat,
@@ -1196,22 +1364,31 @@ define_domain <- function(
     ...
 ) {
 
-  proj <- match.arg(proj)
-  proj <- switch(
-    proj,
-    "lambert"       = "lcc",
-    "mercator"      = "merc",
-    "lalo"          = ,
-    "longlat"       = "latlong",
-    "omerc"         = ,
-    "somerc"        = "tmerc",
-    "rot_longlat"   = ,
-    "rot_latlong"   = ,
-    "RotLatLon"     = "ob_tran",
-    "stereo"        = ,
-    "stereogrpahic" = "stere",
-    proj
-  )
+  is_proj_str <- FALSE
+  projTry <- try(match.arg(proj), silent = TRUE)
+  if (inherits(projTry, "try-error")) {
+    if (length(proj) == 1 && substring(proj, 1, 6) == "+proj=") {
+      is_proj_str <- TRUE
+    }
+  }
+  if (!is_proj_str) {
+    proj <- match.arg(proj)
+    proj <- switch(
+      proj,
+      "lambert"       = "lcc",
+      "mercator"      = "merc",
+      "lalo"          = ,
+      "longlat"       = "latlong",
+      "omerc"         = ,
+      "somerc"        = "tmerc",
+      "rot_longlat"   = ,
+      "rot_latlong"   = ,
+      "RotLatLon"     = "ob_tran",
+      "stereo"        = ,
+      "stereogrpahic" = "stere",
+      proj
+    )
+  }
 
   if (length(nxny) == 1) nxny <- rep(nxny, 2)
   if (length(dxdy) == 1) dxdy <- rep(dxdy, 2)
@@ -1234,25 +1411,29 @@ define_domain <- function(
     exey <- correct_length(round.POSIXt(exey))
   }
 
-  projection <- switch(
-    proj,
-    "lcc" = list(
-      proj = proj, lon_0 = ref_lon,
-      lat_0 = ref_lat, lat_1 = ref_lat, lat_2 = ref_lat
-    ),
-    "merc" = list(
-      proj = proj, lon_0 = ref_lon
-    ),
-    "tmerc" = tranverse_mercator(ref_lon, ref_lat, tilt),
-    "latlong" = list(proj = proj),
-    "ob_tran" = list(
-      proj = proj, o_proj = "latlong",
-      o_lat_p = -ref_lat, o_lon_p = 0, lon_0 = ref_lon
-    ),
-    "stere" = list(proj = proj, lon_0 = ref_lon, lat_0 = ref_lat)
-  )
+  if (is_proj_str) {
+    projection <- meteogrid::proj4.str2list(proj)
+  } else {
+    projection <- switch(
+      proj,
+      "lcc" = list(
+        proj = proj, lon_0 = ref_lon,
+        lat_0 = ref_lat, lat_1 = ref_lat, lat_2 = ref_lat
+      ),
+      "merc" = list(
+        proj = proj, lon_0 = ref_lon
+      ),
+      "tmerc" = tranverse_mercator(ref_lon, ref_lat, tilt),
+      "latlong" = list(proj = proj),
+      "ob_tran" = list(
+        proj = proj, o_proj = "latlong",
+        o_lat_p = -ref_lat, o_lon_p = 0, lon_0 = ref_lon
+      ),
+      "stere" = list(proj = proj, lon_0 = ref_lon, lat_0 = ref_lat)
+    )
 
-  projection <- c(projection, list(R = R, ...))
+    projection <- c(projection, list(R = R, ...))
+  }
 
   res <- list(
     projection = projection, nx = nxny[1], ny = nxny[2],