uasimg: Drone Image Utilities

Documented in uas_worldfile

#' Create world and projection files for UAS images
#'
#' Create world files and projection files for UAS images
#'
#' @param x A list of class 'uas_info'
#' @param flt Which flight(s) in x to process (character or numeric vector, default is all)
#' @param idx Which rows in x to process (default is all)
#' @param aux.xml Create an aux.xml file, logical. See details.
#' @param wld Create a world file, logical. See details.
#' @param wldext Extension for the world file, character. Ignored if \code{wld = FALSE}.
#' @param prj Create a prj file, logical. See details.
#' @param rotated Compute parameters that replicate the camera Yaw.
#' @param quiet Show messages.
#'
#' @details
#' This function creates \href{https://en.wikipedia.org/wiki/Sidecar_file}{sidecard}
#' \href{https://en.wikipedia.org/wiki/World_file}{world files} and/or
#' projection files that are required to view raw images in GIS software such as
#' ArcGIS and QGIS.
#'
#' The parameters this function uses to generate world file are taken exclusively from the image
#' metadata, including the altitude above the launch point, the focal length, and so on. They should be considered
#' \strong{approximate at best}. For a more accurate image placement, process the images with photogrammetry software.
#'
#'
#' To generate world files, you must tell it to compute footprints when generating the flight metadata object. In other words, when you
#' run \code{\link{uas_info}}, pass \code{fp = TRUE}. If you've already generated a metadata object for the flight, you
#' may also need to pass \code{update_cache = TRUE}. Computing footprints requires an altitude for each image, which not
#' all drones support (especially multispectral cameras). For details see \code{\link{uas_info}}.
#'
#' This function has been tested with JPG files from several DJI cameras. It has not yet been adapted for TIF files from
#' multispectral cameras, and may not work with those formats (contact the package author if you want to try).
#'
#'  \code{aux.xml} files are
#'  \href{http://desktop.arcgis.com/en/arcmap/latest/manage-data/raster-and-images/auxiliary-files.htm}{ESRI auxillary files}
#' for raster layers.  If your objective is to open the images in ArcGIS or QGIS, then generating
#' the \code{aux.xml} files (the default setting) should be all you need. To be read by GIS software, they should
#' have the same name as the image file with 'aux.xml' added on. aux.xml files are capable of storing a lot of info
#' about a raster layer, including statistics, the rotation info, coordinate
#' reference system, and other stuff. The files generated by this function will only contain projection
#' and rotation info only. Note also \strong{any existing aux.xml
#' files will be overwritten}. aux.xml files are the only sidecar file that ArcGIS
#' software seems to support for reading the projection info.
#'
#' World files are small text files with extensions \code{jpw} and \code{tfw} for
#' JPG and TIF files respectively. To be read by GIS software, they must
#' have the same basename as the image and be saved in the same directory. Both
#' ArcGIS and QGIS read world files, however they are not needed if the same
#' info is available in an aux.xml file.
#'
#' \code{prj} files contain just the Coordinate Reference System info. They do not
#' seem to be recognized for rasters by ArcGIS, however QGIS picks them up.
#'
#' @return A list of vectors of filenames generated.
#'
#' @seealso \code{\link{uas_info}}
#'
#' @importFrom sf st_crs st_coordinates
#' @import crayon
#' @export

uas_worldfile <- function(x, flt = NULL, idx = NULL, aux.xml = TRUE, wld = FALSE, wldext = "auto", prj = FALSE,
                             rotated = TRUE, quiet = FALSE) {

    if (!inherits(x, "uas_info")) stop("x should be of class \"uas_info\"")
    reslt <- list()

    wldext_lst <- list("jpg" ="jpw", "tif" = "tfw", "bil" = "blw")

    if (is.null(flt)) {
      flt_idx_use <- 1:length(x)
    } else {
      if (is.numeric(flt)) {
        if (max(flt) > length(x)) stop("flt should not be larger than the number of flights saved in the uas image collection object")
        flt_idx_use <- flt
      } else if (is.character(flt)) {
        if (FALSE %in% (flt %in% names(x))) stop("flight name not found in the uas image collection object")
        flt_idx_use <- which(names(x) %in% flt)
      } else {
        stop("Invalid value for flt")
      }
    }

    for (flt_idx in flt_idx_use) {

      if (identical(x[[flt_idx]]$pts, NA)) {
        warning(paste0("Centroids not found for ", names(x)[flt_idx], ". Can not create world file."))
        next
      }

      ## Make sure there's footprint info saved
      if (is.null(x[[flt_idx]]$fp)) {
        stop("Image footprints are required to create worldfiles, but were not found in the image info object. Try re-running `uas_info()` with fp=TRUE")
      }

      if (identical(x[[flt_idx]]$fp, NA)) {
        warning(paste0("Can not generate world file(s). ", names(x)[flt_idx], " does not have footprints."))
        next
      }

      if (!"gsd" %in% names(x[[flt_idx]]$pts)) {
        warning(paste0("Can not generate world file(s). ", names(x)[flt_idx], " does not have gsd."))
        next
      }

      files_gen <- NULL

      ## Get the CRS which will be used to generate the prj files.
      img_wkt <- (x[[flt_idx]]$pts %>% st_crs())$proj4string

      ## Loop through the rows
      if (is.null(idx)) {
        idx_use <- 1:nrow(x[[flt_idx]]$pts)
      } else {
        idx_use <- idx
      }
      for (i in idx_use) {

        ## Get the input file name (minus path)
        img_fn_in <- x[[flt_idx]]$pts[i, "file_name", drop = TRUE]
        img_ext <- tolower(substr(img_fn_in, nchar(img_fn_in)-2, nchar(img_fn_in)))
        img_fnfull <- x[[flt_idx]]$pts[i, "img_fn", drop = TRUE]

        ## Extract the GSD in map units (meter)
        img_gsd_m <- x[[flt_idx]]$pts[i, "gsd", drop = TRUE] / 100

        ## Extract the footprint width, height, and center
        img_fp_width <- x[[flt_idx]]$fp[i, "fp_width", drop = TRUE]
        img_fp_height <- x[[flt_idx]]$fp[i, "fp_height", drop = TRUE]
        #img_ctr <- coordinates(x[[flt_idx]]$pts[i,])
        img_ctr <- x[[flt_idx]]$pts %>% slice(i) %>% sf::st_coordinates()

        if (rotated) {
          ## Extract the yaw
          img_yaw <- x[[flt_idx]]$pts[i, "yaw", drop = TRUE]

          ## Convert compass angle to radians of rotation on the Cartesian plane
          theta <- (180 - img_yaw) *  pi / 180

          ## Define the rotation matrix
          (rot_mat <- matrix(data=c(cos(theta), -sin(theta), sin(theta),
                                    cos(theta)),
                             nrow=2, byrow=TRUE))
        } else {
          theta <- 0
        }

        ## Find the coordinates of the lower right corner relative to the ctr
        lr_unrotated_xy <- matrix(data=c((img_fp_width / 2) - (img_gsd_m / 2),
                                          (- img_fp_height / 2) + (img_gsd_m / 2)),
                                   byrow=TRUE, ncol=2)

        if (rotated) {
          ## Compute rotated coordinates of the lower right
          lr_rotated_xy <- t(rot_mat %*% t(lr_unrotated_xy))

          ## Compute the absolute coordinates of the LR
          lr_rot_ctr_x <- img_ctr[1] + lr_rotated_xy[1]
          lr_rot_ctr_y <- img_ctr[2] + lr_rotated_xy[2]

        } else {
          ## Compute the absolute coordinates of the LR
          lr_rot_ctr_x <- img_ctr[1] + lr_unrotated_xy[1]
          lr_rot_ctr_y <- img_ctr[2] + lr_unrotated_xy[2]
        }

        wld_params <- c(- img_gsd_m * cos(theta),
                        - img_gsd_m * sin(theta),
                        - img_gsd_m * sin(theta),
                        img_gsd_m * cos(theta),
                        lr_rot_ctr_x,
                        lr_rot_ctr_y)

        # Line 1: A: x-component of the pixel width (x-scale)
        # Line 2: D: y-component of the pixel width (y-skew)
        # Line 3: B: x-component of the pixel height (x-skew)
        # Line 4: E: y-component of the pixel height (y-scale), typically negative
        # Line 5: C: x-coordinate of the center of the original
        #            image's upper left pixel transformed to the map
        # Line 6: F: y-coordinate of the center of the original image's
        #            upper left pixel transformed to the map

        if (aux.xml) {
          ## Write the aux.xml
          xml_fn <- paste0(img_fnfull, ".aux.xml")

          ## Note the contents of the GeoTransform tag are the same as the
          ## six coefficients from the world file, but in a different order
          ## c(5, 1, 3, 6, 2, 4)
          writeLines(paste0("<PAMDataset>\n<SRS>", img_wkt, "</SRS>\n",
                            "<GeoTransform>", paste(wld_params[c(5, 1, 3, 6, 2, 4)],
                                                    collapse=", "), "</GeoTransform>\n",
                            "</PAMDataset>"),
                     xml_fn)
          files_gen <- c(files_gen, xml_fn)
        }

        if (wld) {
          # Write the world file
          if (wldext == "auto") {
            if (img_ext %in% names(wldext_lst)) {
              wld_fn_ext <- wldext_lst[[img_ext]]
            } else {
              warning(paste0("World file extension for ", img_ext, " files not known. Resorting to '.wld'."))
              wld_fn_ext <- "wld"
            }
          } else {
            wld_fn_ext <- wldext
          }

          wld_fn <- paste0(substr(img_fnfull, 0, nchar(img_fnfull) - 3),
                           wld_fn_ext)
          writeLines(as.character(wld_params), wld_fn)
          files_gen <- c(files_gen, wld_fn)
        }

        if (prj) {
          ## Write the prj file
          prj_fn <- paste0(substr(img_fnfull, 0, nchar(img_fnfull) - 3), "prj")
          writeLines(img_wkt, prj_fn)
          files_gen <- c(files_gen, prj_fn)
        }

        #wld_params[c(5, 1, 3, 6, 2, 4)]

      } # for 1 in 1:nrow

      ## Append this list of world files to the result
      reslt[[names(x)[flt_idx]]] <- files_gen

    } #for (flt_idx in 1:length(x))

    if (!quiet) message(crayon::green("Done."))

    invisible(reslt)
}