Transformations

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: crstools
 Title: Tools to Work with Projections
-Version: 1.0.0.9002
+Version: 1.0.0.9004
 Authors@R: c(
     person("Andrea", "Manica", , "am315@cam.ac.uk", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-1895-450X")),
@@ -21,11 +21,14 @@ Description: Choose the appropriate map projection for a given application,
   Pozzi et al. (2026) <doi:10.31223/X5K17P> (pre-print).
 License: GPL (>= 3)
 Encoding: UTF-8
+URL: https://github.com/EvolEcolGroup/crstools,
+    https://evolecolgroup.github.io/crstools/
+BugReports: https://github.com/EvolEcolGroup/crstools/issues
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
 Depends:
     R (>= 3.5.0)
 Imports:
+    fields,
     sf,
     ggplot2,
     terra,
@@ -41,7 +44,10 @@ Suggests:
     rnaturalearthdata,
     tidyterra,
     svglite,
-    spelling
+    spelling,
+    scales
 Config/testthat/edition: 3
 VignetteBuilder: knitr
 Language: en-GB
+Config/roxygen2/version: 8.1.0
+RoxygenNote: 7.3.3

NAMESPACE

@@ -5,5 +5,7 @@ export(extract_coords)
 export(find_gcp_coords)
 export(geom_tissot)
 export(georeference_img)
+export(get_gcp_residuals)
+export(get_pts_coords)
 export(plot_gcp)
 export(suggest_crs)

NEWS.md

@@ -1,5 +1,10 @@
 # crstools dev
 * implement plotting function for GCPs.
+* Allow to choose transformation method when georeferencing rasters.
+* Add `get_points_coords()` function to georeference individual points with 
+  a chosen transformation method.
+* Add `get_gcp_residuals()` function` to get the residuals of polynomials to
+  assess the quality of the transformation.
 
 # crstools 0.1.0
 * Initial CRAN submission.

R/choose_gcp.R

@@ -27,12 +27,13 @@
 #'   corresponding geographic coordinates.
 #' @export
 #'
-#' @examplesIf rlang::is_interactive() # Get the path to an example image
-#'   included in the package and choose GCPs img_path <-
-#'   system.file("extdata/europe_map.jpeg", package = "crstools") # this will
-#'   open a new window where you can choose some points gcp_europe <-
-#'   choose_gcp(img_path) # after the first set of points is chosen, we can add
-#'   more points gcp_europe <- choose_gcp(img_path, gcp = gcp_europe)
+#' @examplesIf rlang::is_interactive() 
+#'   # Get the path to an example image included in the package and choose GCPs
+#'   img_path <- system.file("extdata/europe_map.jpeg", package = "crstools") 
+#'   # this will open a new window where you can choose some points
+#'   gcp_europe <- choose_gcp(img_path)
+#'   # after the first set of points is chosen, we can add more points 
+#'   gcp_europe <- choose_gcp(img_path, gcp = gcp_europe)
 choose_gcp <- function(image_obj, gcp = NULL, col = "red") {
   # check if image is a file path or an array
   if (is.character(image_obj)) {

R/find_gcp_coords.R

@@ -1,16 +1,16 @@
 #' Find the coordinates (longitude and latitude) of the ground control points
 #' (GCPs) in a given image.
 #'
-#'
-#' @details NOTE: There are two conventions on how to define pixel coordinates. In this
+#' @details NOTE: There are two conventions on how to define pixel coordinates.
+#' In this
 #' function (and more generally throughout `crstools`), the origin is defined as
 #' the bottom left corner of the image, with x increasing to the right and y
 #' increasing upwards. This is consistent with the convention used in many image
 #' processing libraries. However, some libraries (like OpenCV and GDAL) define
 #' the origin at the top left corner, with y increasing downwards. Be sure to
 #' check which convention your image processing library uses when working with
 #' pixel coordinates.
-
+#'
 #' @param gcp A data frame containing the GCPs with columns `id`, `x`, `y`,
 #'   `longitude`, and `latitude`.
 #' @param sf_obj the reference map, as an sf object (already cut to the extend

R/georeference_img.R

@@ -17,36 +17,60 @@
 #'   - `id`: An identifier for each GCP (numeric).
 #'   - `x`: The x-coordinate of the GCP (in pixel space).
 #'   - `y`: The y-coordinate of the GCP (in pixel space).
-#'   - `lon`: The longitude of the GCP (georeferenced).
-#'   - `lat`: The latitude of the GCP (georeferenced).
+#'   - `longitude`: The longitude of the GCP (georeferenced).
+#'   - `latitude`: The latitude of the GCP (georeferenced).
 #'
 #' @param output_path A character string representing the file path to the input
 #'   image. (`_warp.tif`) will be appended to it.
 #'
+#' @param transform_method A character string specifying the transformation
+#'   method to be used for warping the image. Options are "poly_1" (first order
+#'   polynomial), "poly_2" (second order polynomial), "poly_3" (third order
+#'   polynomial), "tps" (thin plate spline), or "auto" (the default, allowing
+#'   GDAL to choose of a polynomial of the appropriate order based on the
+#'   number of available GCP. Polynomials are best for standard maps directly
+#'   captured from a publication (a first or second order polynomial is often
+#'   sufficient), tps allows for scanning artefacts, but it is badly affected
+#'   by any incorrect GCP.
+#'
 #' @return A character string representing the path to the newly created warped
 #'   TIFF image file (`_warp.tif`). This file contains the georeferenced image.
 #'
 #' @export
 #'
 #' @examplesIf rlang::is_interactive()
 #' # get the path to an example image included in the package
-#' img_path <- system.file("extdata/europe_map.jpeg", package = "crstools")
+#' img_path <- system.file("extdata/europe_map.jpeg",
+#'   package = "crstools"
+#' )
 #' # load a set of GCPs (or we could create them using the choose_gcp()
 #' # and find_gcp() functions)
 #' gcp_df <- readRDS(system.file(
 #'   "extdata/europe_gcp_georef.RDS",
 #'   package = "crstools"
 #' ))
-#' #' # Assuming you have a set of GCPs in gcp_df and an image file "image.jpg"
+#' # Assuming you have a set of GCPs in gcp_df and an image file "image.jpg"
 #' warped_img <- georeference_img(
 #'   image_obj = img_path, gcp = gcp_df,
 #'   output_path = tempfile(
-#'     patter = "georef_img_",
-#'     tmpdir = tempdir(),
+#'     pattern = "georef_img_", tmpdir = tempdir(),
 #'     fileext = ".tif"
 #'   )
 #' )
-georeference_img <- function(image_obj, gcp, output_path = NULL) {
+georeference_img <- function(image_obj, gcp, output_path = NULL,
+                             transform_method = c(
+                               "auto", "poly_1", "poly_2",
+                               "poly_3", "tps"
+                             )) {
+  transform_method <- match.arg(transform_method)
+  # now convert transform method into the appropriate GDAL option
+  gdal_transform_option <- switch(transform_method,
+    "poly_1" = c("-order", "1"),
+    "poly_2" = c("-order", "2"),
+    "poly_3" = c("-order", "3"),
+    "tps" = "-tps",
+    "auto" = NULL
+  )
   # check if gcp is a dataframe with the right columns
   # nolint start
   if ((!is.data.frame(gcp)) ||
@@ -109,19 +133,25 @@ georeference_img <- function(image_obj, gcp, output_path = NULL) {
     options = c(as.vector(t(cbind("-gcp", gcp[, -1]))), "-of", "GTiff")
   )
 
+  warp_options <- c(
+    "-s_srs",
+    "EPSG:4326",
+    "-t_srs",
+    "EPSG:4326",
+    "-overwrite"
+  )
+
+  # if we have a transform method, add it to the option
+  if (!is.null(gdal_transform_option)) {
+    warp_options <- c(gdal_transform_option, warp_options)
+  }
+
   # Warp the image into a spatial reference system (EPSG:4326)
   sf::gdal_utils(
     "warp",
     source = map_tif,
     dest = map_warp_tif,
-    options = c(
-      "-tps",
-      "-s_srs",
-      "EPSG:4326",
-      "-t_srs",
-      "EPSG:4326",
-      "-overwrite"
-    )
+    options = warp_options
   )
 
   return(map_warp_tif)