R/lxy.exp.mov.R

#' Create an animation from a LoCoH-xy object
#'
#' @param lxy A \link{LoCoH-xy} object
#' @param id The id value(s) to be on the plot
#' @param all.ids.at.once Display all the individual ids simultaneously. If False, an animation will be created for each id. T/F.
#' @param all.ids.col.unique Whether to use unique colors for each individual when plotting multiple individuals simultaneously, T/F
#' @param all.ids.col A named list of color values; the element names must match the name(s) of the ids
#' @param all.ids.legend Where to place a legend showing the color of each id (for an animation showing the movement of multiple individuals simultaneously): 'bottomright', 'bottom', 'bottomleft', 'left', 'topleft', 'top', 'topright', 'right', or 'center'. May also be \code{NULL}, in which case the legend will not be displayed. Ignored if \code{ids.legend.bln=FALSE} or \code{all.ids.col.unique = FALSE}.
#' @param all.ids.legend.cex The character expansion factor for the id legend. See parameter \code{all.ids.legend} above.
#' @param dt.start The starting date-time. An object of class POSIXt or one that can be coerced to POSIXt. If NULL, the earliest date-time in the series will be used.
#' @param dt.end The ending date-time. An object of class POSIXt or one that can be coerced to POSIXt. If NULL, the last date-time in the series will be used.
#' @param frame.method How each frame should be defined temporally, "time" - each frame represents a fixed amount of time, "location" each frame is a point in the series
#' @param frame.rtd The real-time duration of each frame (in seconds). If "auto" (default), the lowest median sampling frequency will be used
#' @param xlim A two-element numeric vector for the range of the x-axis (in map units)
#' @param ylim A two-element numeric vector for the range of the y-axis (in map units)
#' @param dt.label Add a label for the date of the frame
#' @param dt.label.x The x-coordinate for the date label (ignored if \code{dt.label=FALSE})
#' @param dt.label.y The y-coordinate for the date label (ignored if \code{dt.label=FALSE})
#' @param dt.label.col A color value/name for the date label (ignored if \code{dt.label=FALSE})
#' @param dt.label.bg The background color for the date label. Set to \code{NA} for a transparent background. (ignored if \code{dt.label=FALSE})
#' @param title A title for the map
#' @param title.show Whether to show the title. T/F.
#' @param axes.show Whether to show the axes (ticks, labels and titles). Can be over-written by \code{axes.ticks} and \code{axes.titles}. T/F.
#' @param axes.ticks Whether to show the tick marks and labels on the axes. T/F.
#' @param axes.titles Whether to show axes titles. T/F.
#' @param mar.map Margin settings for the map, see \code{\link{par}}
#' @param mgp.map Locations of the axes elements for the map, see \code{\link{par}}
#' @param col.xys.active Color of the active point. Ignored if \code{col.by.hour.of.day = TRUE}.
#' @param col.xys.background Color of the non-active points. To hide non-active points, set this to \code{NA}.
#' @param cex.xys.active The character expansion factor of the active point
#' @param cex.xys.background The character expansion factor of non-active points
#' @param tz.local The name of the time zone of the study area
#' @param tz.local.check Check whether tz.local is a valid timezone name (not implemented) T/F.
#' @param col.by.hour.of.day Whether to color the active point by the hour of day (i.e., dark colors at night, orange for day time locations). T/F
#' @param col.hod A vector of 24 color values used to symbolize the color of the active point by the hour of day. Ignored if \code{col.by.hour.of.day = FALSE}.
#' @param screen.test Create up to three sample frame(s) on the screen (instead of PNG files)
#' @param width The width of each frame in pixels (if \code{screen.test=FALSE}) or inches (if \code{screen.test=TRUE}).
#' @param height The height of each frame in pixels (if \code{screen.test=FALSE}) or inches (if \code{screen.test=TRUE}).
#' @param max.frames The maximum number of frames to produce.
#' @param png.pointsize The pointsize (in pixels) for the PNG image, equivalent to the height or width of a character in pixels (increase to make labels appear larger)
#' @param tmp.dir A directory where temporary PNG files for each frame will be created, character.
#' @param tmp.files.delete Delete the temporary PNG files when done, T/F
#' @param prompt.continue Whether to present a summary of the encoding settings and get user confirmation before continuing, T/F
#' @param fn.mov The path and filename of output file (either *.mov or *.mp4). If NULL a filename will be automatically generated
#' @param fn.mov.dir The directory where the animation will be saved (ignored if a value for \code{fn.mov} is passed)
#' @param fn.mov.exists What to do if the animation file already exists: "auto.increment", "overwrite", "stop", or "ask". 
#' @param duration The desired duration of the animation (in seconds)
#' @param fps A numeric value for frames per second 
#' @param skip Output every nth frame. To include every frame set skip=1. Integer.
#' @param ffmpeg The name of the ffmpeg executable. If NULL (the default), a default file name will be used. See notes.
#' @param create.mov Whether to actually create the mov file. Set to FALSE preview a few frames without actually encoding them.
#' @param show.cmd Whether to display the ffmpeg command line that stitches the images into an animation
#' @param fmt Video format: \code{'mov'} (Quicktime animation codec) or \code{'mp4'} (h.264)
#' @param info.only Only return info 
#' @param shp.csv The path and filename of a csv file that contains information about shapefiles, including layer names, file, and symbology.
#' @param layers The name(s) of layers in shp.csv to display in the background. Will be displayed using the symbology in shp.csv. Character vector or comma delimited string.
#' @param tiff.fn The path and name of a GeoTIFF file (e.g., satellite image) that will be displayed in the background. See notes.
#' @param tiff.bands A vector of threee integers corresponding to the bands of the GeoTIFF image that will be mapped to the red, green and blue color guns respectively.
#' @param tiff.col A vector of color values for plotting single-band images in the background. Ignored if using three bands.
#' @param tiff.pct Whether or to convert the GeoTIFF to an indexed 256 color RGB image, which may speed up drawing. T/F.
#' @param tiff.buff A numeric buffer distance in map units that the range of the plot will be expanded so the points are not right on the edge of the GeoTIFF.
#' @param tiff.fill.plot Whether to fill the entire plot area with the GeoTIFF. T/F.
#' @param bg2png Save the plot background elements as a static raster image (to improve speed), ignored if \code{screen.test=TRUE}
#' @param crop.layers.to.extent Whether to crop the shapefile layers to the view extent (may speed up drawing time)
#' @param date.bar The height of the lower section of the plot to devote to the time bar, in inches. To hide the time bar completely, set \code{date.bar=0}.
#' @param date.bar.bins The target number of bins (tick marks + 1) on the time bar (integer)
#' @param col.db A single color value for the date bar axes / tick labels, character
#' @param cex.axis.db Character expansion factor for the labels on the date bar axis.
#' @param beep Beep when done, T/F
#' @param report.time Show the time taken when done, T/F
#' @param status Show progress bar and status messages 
#'
#' @note This function creates an animation from a LoCoH-xy object. There are three general steps in the workflow: 
#' 1) design the frame layout, 2) choose values for the speed and duration of the animation, and 3) create the 
#' animation. 
#'
#' One will normally want to run the function a few times without actually encoding to tweak the frame design (e.g., where 
#' the date label and legend appear). To see what a frame in the output will approximately look like, 
#' set \code{screen.test=TRUE}. See Appendix VI of the T-LoCoH tutorial for further details.
#' 
#' If you have locations for all hours of the day, you can color the active point by the hour-of-day by setting
#' \code{col.by.hour.of.day=TRUE}. This can help you visualize daily patterns in movement, assuming that the 
#' time-stamps of each point are in local time. Locations at night will appear dark, sunrise and sunset reddish, and mid-day 
#' points blue (these colors can be adjusted by passing a different value for \code{col.hod}. 
#'
#' To create the animation, two and only two of the following parameters must be passed: \code{duration}, \code{fps}, and \code{skip}. 
#' The third parameter will be computed based on the other two. To include every frame, pass \code{fps}, set \code{skip=1}, and leave \code{duration} out.
#'
#' Larger values for \code{fps} will result in the animation running 'faster'. Values between 10 and 20 often work well; 
#' beyond 30 fps the eye can't keep up with the motion
#' Note if you pass values for \code{fps} and \code{duration}, an appropriate value for \code{skip} will be computed but 
#' the final duration of the animation may not be exactly equal to \code{duration} because only interger values of \code{skip} are allowed.
#'
#' \code{tmp.dir="."} (or another folder), \code{create.mov=FALSE}, and \code{tmp.files.delete=FALSE}. This will generate a few sample frames as PNG files 
#' but not delete them so you can inspect them using an image viewer. Once these look good, create the full animation by setting \code{create.mov=TRUE} and 
#' \code{max.frames=NULL}. 
#'
#' If frame.method is 'auto', the script will use time-based frames when multiple individuals are being animated simultaneously, and location-based frames otherwise.
#'
#' \code{duration} (the duration of the animation in seconds) should not be confused with \code{frame.rtd} which is the real-time duration 
#' of a single frame in seconds (e.g., \code{frame.rtd=3600} means each frame will represent 1 hour).
#'
#' If \code{date.bar} is too small or too large, you might get a 'margins too large' error. Try values around 1, or 
#' hide the date bar completely by setting \code{date.bar=0}.
#'
#' The output animation is encoded in QuickTime format or mp4. The Quicktime file is encoded using Quicktime's 'animation' codec, a lossless format that 'scrubs' 
#' well (i.e., you can drag the scroll bar 
#' to view frame by frame). The mp4 settings use the h.264 compression algorithm with a constant rate factor of 20 (good quality). In both cases, 
#' encoding requires installing the open source encoding program ffmpeg. ffmpeg is a command line program that 
#' can be downloaded from \url{http://ffmpeg.org/download.html} (select the package for your OS). 
#' After downloading, the ffmpeg executable file should be saved in the current working directory, or a directory on the PATH environment variable (e.g., on windows c:\\windows).
#' You may also pass the path and name to ffmpeg using the \code{ffmpeg} argument.
#'
#' If ffmpeg is not available, you can still use this function to generate the individual frames and then use another utility (e.g., ImageMagick, Quicktime Pro) 
#' to combine the frames into a video file. For best results use a 'lossless' compression method in the encoding program.
#' To create the individual frames only for encoding with another utility, set \code{tmp.dir="."} (the working directory) and \code{tmp.files.delete=FALSE}. 
#' Passing \code{show.cmd=TRUE} will display the ffmpeg command line that will generate animation, which can be useful if you want to do some additional
#' editing to the frames prior to encoding.
#'
#' If \code{fn.mov.exists = "auto.increment"}, a two-digit number will be appended to the filename to avoid overwriting an existing file
#'
#' @return A list with information about each file created. Each element of the list is another list with two 
#' elements: \code{fn} (the full filename) and \code{dim} (a two-element numeric vector with the frame width and height). If no animation file(s) were
#' created, returns NULL. 
#'
#' @export
#' @import sp

lxy.exp.mov <- function(lxy, id=NULL, all.ids.at.once=TRUE, all.ids.col.unique=all.ids.at.once, all.ids.col=NULL, 
                        all.ids.legend=c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center")[5],
                        all.ids.legend.cex=0.8, dt.start=NULL, dt.end=NULL, frame.method=c("auto", "time", "location")[1], frame.rtd="auto", xlim=NULL, ylim=NULL, 
                        dt.label=TRUE, dt.label.x=NULL, dt.label.y=NULL, dt.label.col="black", dt.label.bg="gray90",
                        title=NULL, title.show=TRUE,
                        axes.show=TRUE, axes.ticks=axes.show, axes.titles=FALSE, 
                        mar.map=c(0.7 + (if (axes.ticks) 0.9 else 0) + (if (axes.titles) 1.3 else 0),
                                  0.5 + (if (axes.ticks) 0.9 else 0) + (if (axes.titles) 1.3 else 0), 
                                  if (title.show) 2.1 else 0.5, 
                                  0.5), 
                        mgp.map=c(0.4 + if (axes.ticks) 1.2 else 0, 0.4, 0),
                        col.xys.active="red", col.xys.background="gray80", 
                        cex.xys.active=1.5, cex.xys.background=0.5, 
                        tz.local=NULL, tz.local.check=TRUE, 
                        col.by.hour.of.day=FALSE, col.hod=colorRampPalette(colors()[c(24,30,553,121,26,121,553,30,24)])(24),
                        width=if (screen.test) 7 else 608, height=NULL, max.frames=NULL, png.pointsize=16+(width-480)/80, 
                        screen.test=FALSE, tmp.dir=NULL, tmp.files.delete=TRUE, prompt.continue=TRUE,
                        fn.mov=NULL, fn.mov.dir=getwd(), fn.mov.exists=c("auto.increment", "overwrite", "stop", "ask")[1], 
                        duration=NULL, fps=NULL, skip=NULL, ffmpeg=NULL, create.mov=TRUE, show.cmd=FALSE, fmt=c("mov","mp4")[1], info.only=TRUE, 
                        shp.csv=NULL, layers=NULL, tiff.fn=NULL, tiff.bands=c(3,2,1), tiff.col=gray(0:255/255), tiff.pct=FALSE, tiff.buff=0, tiff.fill.plot=TRUE, 
                        bg2png=!is.null(layers), crop.layers.to.extent=TRUE, 
                        date.bar=0.85, date.bar.bins=12, col.db="darkblue", cex.axis.db=0.7, 
                        beep=FALSE, report.time=TRUE, status=TRUE) {

    
    
    ## info.only (not implemented), will *not* make the mov, only return a list of the filename, width, and height, 

    ## col.hour.of.day is NULL or a 24 color values for each hour of the day
    ## tz.local can be NULL or a named time zone. If not null, dt will be converted to local time

    ## This will take an lxy object and export it as mov file
    if (!inherits(lxy, "locoh.lxy")) stop("lxy should be of class \"locoh.lxy\"")
    if (bg2png) if(!requireNamespace("png", quietly=TRUE)) stop("Package png required for bg2png")
    if (!fmt %in% c("mp4","mov")) stop("fmt must be 'mov' or 'mp4'")

    if (is.null(lxy[["pts"]][["dt"]])) stop("Can't animate without date values")
    if (is.null(id)) {
        id <- levels(lxy[["pts"]][["id"]])
    } else {
        if (FALSE %in% (id %in% levels(lxy[["pts"]][["id"]]))) stop("id value(s) not found")
    }
    if (length(id)==1) all.ids.at.once <- FALSE
    
    if (!is.null(ylim) && !is.null(height)) stop("You can only pass ylim or height, not both")
    if (!is.null(duration) && !is.null(fps) && !is.null(skip)) stop(cw("You can pass at most two of the following parameters: duration, fps, and skip"), final.cr=F)
    if (!screen.test && is.null(duration) && is.null(fps) && is.null(skip)) stop("You must only pass two of the following parameters: duration, fps, and skip")    
    if (length(col.hod) != 0 && length(col.hod) != 24) stop("col.hod must be a vector with 24 color values or NULL")
    ## if (!is.null(tz.local) && tz.local.check) {if (!tz.valid(tz.local)) stop("Unknown value for tz")}
    if (!is.null(dt.end) && !is(dt.end, "POSIXt")) stop("dt.end must be of class POSIXt")
    if (!is.null(dt.start) && !is(dt.start, "POSIXt")) stop("dt.start must be of class POSIXt")
    if (frame.method=="auto") frame.method <- if (length(id) > 1 && all.ids.at.once) "time" else "location"
    
    #if (!is.null(fn.mov) && length(fn.mov) != length(hs)) stop("The length of fn.out must equal the number of runs in lhs")

    if (is.null(width)) stop("Please supply a value for width")
    if (screen.test) {
        if (is.null(max.frames)) max.frames <- 3
        if (max.frames>10) stop("Can't create more than 10 frames for a screen test")
        if (width > 18) {
            if (status) cat("   Width too big for screen, resetting to 7 \n")
            width <- 7
        }
        round.up.to.nearest <- 1
    } else {
        if (width < 100) {
            ans <- readline(prompt = paste("Do you really want to create an animation that is only", width, "pixels wide? y/n "))
            if (ans != "y") return(invisible(NULL))
        }
        ## Round width up to nearest multiple of 16 to help encoder
        round.up.to.nearest <- 16
        if (width %% 16 != 0) {
            width <- ceiling(width/round.up.to.nearest) * round.up.to.nearest
            if (status) cat("   Increasing width to ", width, " (nearest multiple of 16) for encoding \n", sep="")
        }
        
    }
    
    if (date.bar > 0) {
        ## if (date.bar > 0.5) stop("Date bar should be 0..0.5")
        ## date.bar.cm <- paste(date.bar * 2.54, " cm", sep="")
    }
    
    ## Create a list of the indices 
    if (all.ids.at.once) {
        idx.ids.lst <- list()
        idx.ids.lst[[paste(id, collapse=".", sep="")]] <- list(id=id, idx=which(lxy[["pts"]][["id"]] %in% id))
    
        ## Define a vector of colors that is in the same sequence as levels(lxy[["id"]]) (if needed)
        if (all.ids.col.unique) {
            if (col.by.hour.of.day) stop(cw("If you are displaying each individual with a unique color, you can't also color the points by hour of day."))
            if (is.null(all.ids.col)) {
                ids.cols <- rainbow(length(id), end=5/6)
            } else {
                if (!is.list(all.ids.col) || length(all.ids.col) != length(id)) stop("all.ids.col must be a named list with a color value for each id")
                if (FALSE %in% (names(all.ids.col) %in% levels(lxy[["pts"]][["id"]]))) stop(cw("The names of the elements in all.ids.col must match the names of the individuals in lxy", final.cr=FALSE))
                ids.cols <- rep(NA, nlevels(lxy[["pts"]][["id"]]))
                for (id.name in names(all.ids.col)) ids.cols[which(levels(lxy[["pts"]][["id"]])==id.name)] <- all.ids.col[[id.name]]
            }
            lxy.id.int <- as.numeric(lxy[["pts"]][["id"]])
        } 
    } else {
        idx.ids.lst <- lapply(id, function(x) list(id=x, idx=which(lxy[["pts"]][["id"]] == x)))
        names(idx.ids.lst) <- id
        all.ids.col.unique <- FALSE
    }
    
    ## Error check legend settings
    ids.legend.bln <- all.ids.at.once && all.ids.col.unique && !is.null(all.ids.legend)
    if (ids.legend.bln) {
        if (!all.ids.legend %in% c("bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right", "center")) {
            stop("Invalid value for all.ids.legend")    
        }
        legend.str <- levels(lxy[["pts"]][["id"]])
     }

    ## Create presets for the margins of the map and date bar
    mar.db <- c(2, 5, 0, 2)
    mgp.db <- c(0, 0.4, 0)
    axes.lbl <- colnames(coordinates(lxy[["pts"]]))
    tick.show <- if (axes.ticks) "s" else "n"
    
    ## Get the temp folder
    if (!screen.test) {
        if (is.null(tmp.dir)) {
            tmp.dir <- tempdir()
        } else {
            if (!file.exists(tmp.dir)) {
                dir.made <- dir.create(tmp.dir)
                if (!dir.made) stop(paste("Couldn't make output folder ", tmp.dir, sep=""))            
            }
        }
    }
    

    ## Expand xlim and ylim by tiff.buff 
    if (is.null(tiff.fn)) {
        range.expand <- rep(0,2)
    } else {
        if (!requireNamespace("rgdal", quietly=TRUE)) stop("package rgdal required to display a tiff in the background, please install")
        if (!file.exists(tiff.fn)) stop(paste(tiff.fn, "not found"))
        range.expand <- c(-tiff.buff, tiff.buff)
        tiff.sgdf <- NULL
        if (length(tiff.bands) > 3) stop("tiff.bands can not be longer than 3")
    }

   
    ## Prepare GIS features
    gis.features.full.extent <- shp.layers(layers, shp.csv=shp.csv)
    
    ## Convert date stamps to local time if needed
    if (!is.null(tz.local) && attr(lxy[["pts"]][["dt"]], "tzone") != tz.local) {
        lxy[["pts"]][["dt"]] <- as.POSIXct(format(lxy[["pts"]][["dt"]], tz=tz.local), tz=tz.local)   
    }

    res <- list()
    
    for (idx.set in 1:length(idx.ids.lst)) {
       
        idx.this.loop <- idx.ids.lst[[idx.set]][["idx"]]
        
        ## Define the first and last dates for this group
        dt.start.use <- if (is.null(dt.start)) min(lxy[["pts"]][["dt"]][idx.this.loop]) else dt.start
        dt.end.use <- if (is.null(dt.end)) max(lxy[["pts"]][["dt"]][idx.this.loop]) else dt.end
    
        if (frame.method=="time") {
            ## Get the frame time-window (real time duration, not animation duration)
            if (frame.rtd == "auto") {
                frame.rtd.use <- min(lxy[["rw.params"]][lxy[["rw.params"]][["id"]] %in% idx.ids.lst[[idx.set]][["id"]], "time.step.median"])
            } else {
                frame.rtd.use <- frame.rtd
            }
    
            ## Next define the start time of each frame. dt.frames is a vector of POSIXct objects in local time
            dt.frames <- seq(from=dt.start.use, to=dt.end.use, by=frame.rtd.use)
            
            ## We need to add one more date at the end of dt.frames, which we will delete later, to prevent 
            ## points beyond the range of dates from appearing on the final frame
            ## THIS WILL SCREW-UP THE TIME ZONE ATTRIBUTE, BETTER TO MANUALLY DELETE POINTS PAST DT.END.USE BELOW
            ## if (dt.end.use > dt.frames[length(dt.frames)]) dt.frames <- 
    
            ## Get the frame where each point falls, then make that into a list
            lxyidx.dtframeidx <- findInterval(lxy[["pts"]][["dt"]][idx.this.loop], dt.frames)        
            dt.frames.lxyidx.lst <- lapply(1:length(dt.frames), function(x) which(lxyidx.dtframeidx==x))
            
            ## Delete any points from the very last list element that were sampled after dt.end.use
            dt.frames.lxyidx.lst[[length(dt.frames)]] <- dt.frames.lxyidx.lst[[length(dt.frames)]][lxy[["pts"]][["dt"]][dt.frames.lxyidx.lst[[length(dt.frames)]]] <= dt.end.use]
            
        } else if (frame.method=="location") {
            dt.frames.lxyidx.lst <- as.list(idx.this.loop)
            dt.frames <- lxy[["pts"]][["dt"]][idx.this.loop]
        } else {
            stop("Unknown value for 'frame.method'")
        }
        
        ## Set the spatial range of the entire plot
        xys <- data.frame(coordinates(lxy[["pts"]]))[idx.this.loop, , drop=FALSE]
        xys.range.lower.upper <- as.data.frame(sapply(xys, range))
        height.use <- height

        ## Define xlim and ylim
        if (is.null(xlim)) xlim <- xys.range.lower.upper[,1] + range.expand
        if (is.null(ylim)) {
            ylim <- xys.range.lower.upper[,2] + range.expand
            
            ## Compute frame height based on aspect ratio of points
            if (!is.null(height.use)) {
                ylim.diff.based.on.width.height.ratio <-  diff(xlim) * height.use / width 
                if (ylim.diff.based.on.width.height.ratio > diff(ylim)) {
                    ylim <- mean(ylim) + c(-ylim.diff.based.on.width.height.ratio / 2, ylim.diff.based.on.width.height.ratio / 2)
                } else {
                    ## The passed height is too small to fit all the points, set it to NULL so it will be recculated later
                    ## (Alternative is to reduce ylim)
                    height.use <- NULL
                }
            }
        }

        ## Prepare the tiff
        if (!is.null(tiff.fn)) {
            if (tiff.fill.plot) {
                half.plot.size <- c(-0.5, 0.5) * max(diff(xlim), diff(ylim))
                rx.tiff <- half.plot.size + mean(xlim)
                ry.tiff <- half.plot.size + mean(ylim)
            } else {
                rx.tiff <- xlim
                ry.tiff <- ylim
            }
            
            tiff.sgdf <- readpartgdal(tiff.fn, xlim=rx.tiff, ylim=ry.tiff, band=tiff.bands, silent=TRUE)
            if (is.null(tiff.sgdf)) {
                tiff.fn <- NULL
            } else {
                if (tiff.pct) {
                    if (length(tiff.sgdf@data)!=3) {
                        cat("   Incorrect number of bands, can't convert image to indexed RGB\n")
                        tiff.pct <- FALSE
                    } else if (bg2png) {
                        cat("   Can't convert background image to a RGB when bg2png=TRUE\n")
                        tiff.pct <- FALSE
                    }
                }
                if (tiff.pct) {
                    tiff.sgdf.cols <- rgdal::SGDF2PCT(tiff.sgdf, adjust.bands=FALSE)
                    tiff.sgdf$idx <- tiff.sgdf.cols$idx
                }
            }
            
        }

        
        ## Determine if we're not going to crop gis.layers to the plot extent
        if (crop.layers.to.extent && length(gis.features.full.extent) > 0 ) {
            gis.layers.ready <- FALSE
        } else {
            gis.features <- gis.features.full.extent
            gis.layers.ready <- TRUE
        }
        
        bg.png <- NULL
        fn.background.png <- tempfile(fileext=".png")
        
        ## Get device resolution and height of default-sized characters in inches
        dpi <- if (screen.test)  1 else 72
        no.device.open <- dev.cur() == 1 
        line.height.in <- par("csi")
        if (no.device.open) dev.off()
        
        ## Compute height of the device if needed
        if (is.null(height.use)) {
            top.plus.bottom.margins.pixels <- dpi * sum(mar.map[c(1,3)]) * line.height.in
            left.plus.right.margins.pixels <- dpi * sum(mar.map[c(2,4)]) * line.height.in
            needed.height.of.plot.area.pixels <- (width - left.plus.right.margins.pixels) * diff(ylim) / diff(xlim)
            if (date.bar > 0) needed.height.of.plot.area.pixels <- needed.height.of.plot.area.pixels + (dpi * date.bar)
            height.use <- top.plus.bottom.margins.pixels + needed.height.of.plot.area.pixels
            
            ## If plotting to PNG, round height up to the nearest multiple of 16 (to help the video encoder)
            if (!screen.test) height.use <- ceiling(round(height.use)/round.up.to.nearest) * round.up.to.nearest
        } 
        
        ## Get the height of the date bar as a proportion of the total height
        date.bar.pth <- date.bar / (height.use / dpi)
        
        ## Find a section of the plot area where there are no xy points so we can place the date label there
        ## Basically the approach here is to divide up the plot into nxn grids, then scan through them from top to bottom and find the first 
        ## grid no points in it. If there are no grids without points, the one with with the fewest number of points in it is selected
        ## and that becomes the location of the date label
        if (dt.label) {
            if (is.null(dt.label.x) && is.null(dt.label.y)) {
                divide.axes.into.n <- 5
                label.area.num.pts <- matrix(9999, ncol=divide.axes.into.n, nrow=divide.axes.into.n)
                label.areas.xs <- seq(from=xlim[1], to=xlim[2], by=diff(xlim) / divide.axes.into.n)
                label.areas.ys <- seq(from=ylim[1], to=ylim[2], by=diff(ylim) / divide.axes.into.n)
                label.area.use.idx <- NULL
                for (j in divide.axes.into.n:1) {
                    for (i in 1:divide.axes.into.n) {
                        if (is.null(label.area.use.idx)) {
                            pts.in.this.label.area <- which(xys[,1] >= label.areas.xs[i] & xys[,1] <= label.areas.xs[i+1] & xys[,2] >= label.areas.ys[j] & xys[,2] <= label.areas.ys[j+1])
                            label.area.num.pts[i,j] <- length(pts.in.this.label.area) 
                            if (length(pts.in.this.label.area)==0) {
                                ## First one that's empty is good enough
                                label.area.use.idx <- c(i,j)
                            }
                        }
                    }
                }
                if (is.null(label.area.use.idx)) {
                    label.area.use.idx <- as.numeric(which(label.area.num.pts == min(label.area.num.pts), arr.ind = TRUE)[1,])
                }
                dtlabel.pt <- data.frame(x=mean(label.areas.xs[c(label.area.use.idx[1], label.area.use.idx[1]+1)]), 
                                         y=mean(label.areas.ys[c(label.area.use.idx[2], label.area.use.idx[2]+1)]))
            } else {
                dtlabel.pt <- data.frame(x=dt.label.x, y=dt.label.y)
            }
        }

        ## Compute duration, skip and fps
        duration.use <- NULL
        skip.use <- NULL
        fps.use <- NULL
        
        if (is.null(duration)) {
            if (is.null(fps)) fps.use <- 10 else fps.use <- fps
            if (is.null(skip)) skip.use <- 1 else skip.use <- skip
            frames.to.use <- seq(from=1, to=length(dt.frames.lxyidx.lst), by=skip.use)
            if (!is.null(max.frames) && length(frames.to.use) > max.frames) frames.to.use <- frames.to.use[1:max.frames]
            duration.use <- length(frames.to.use) / fps.use
        } else {
            ## Duration is provided. Either fps or skip is NULL and needs to be computed
            duration.use <- duration
            if (is.null(fps)) {
                if (is.null(skip)) skip.use <- 1 else skip.use <- skip
                frames.to.use <- seq(from=1, to=length(dt.frames.lxyidx.lst), by=skip.use)
                fps.use <- length(frames.to.use) / duration.use
                if (fps.use > 60) stop("Required frame rate higher than maximum value of 60. Please increase the duration or skip parameter")
                ## fps can be a decimal value, so no need to recompute duration
            } else {
                fps.use <- fps
                ## Frame rate and duration are both provided by the user, need to compute skip and then recompute duration
                skip.use <- floor(length(dt.frames.lxyidx.lst) / (fps.use * duration.use))
                if (skip.use < 1) stop("there are not enough locations to produce enough frames for that duration and frame rate")
                frames.to.use <- seq(from=1, to=length(dt.frames.lxyidx.lst), by=skip.use)
                duration.use <- length(frames.to.use) / fps.use
            }
        }
        
        # If needed, reduce the number of frames, but we won't compute the frame rate
        if (!is.null(max.frames) && length(frames.to.use) > max.frames) frames.to.use <- frames.to.use[1:max.frames]
        if (title.show) {
            if (is.null(title)) {
                str.title <- paste(idx.ids.lst[[idx.set]][["id"]], collapse=", ", sep="")
            } else {
                str.title <- title
            }
        } else {
            str.title <- NULL
        }

        if (screen.test) {
            if (.Platform$OS.type == "windows") {
          		  windows(width=width, height=height.use, record=T)
       	    } else {
                dev.new(width=width, height=height.use)
            }	
            opar <- par(mar=mar.map, mgp=mgp.map)
            if (date.bar > 0) layout(matrix(1:2, ncol=1), heights=c(1 - date.bar.pth, date.bar.pth))
            if (bg2png) bg2png <- FALSE
        } else {
            ## Construct a filename for the mov
            fn.mov.base.exp <- expression(paste(paste(idx.ids.lst[[idx.set]][["id"]], collapse=".", sep=""), 
                                                ".", format(min(lxy[["pts"]][["dt"]][idx.ids.lst[[idx.set]][["idx"]]]), format="%Y-%m-%d"), 
                                                ".", format(max(lxy[["pts"]][["dt"]][idx.ids.lst[[idx.set]][["idx"]]]), format="%Y.%m.%d"), 
                                                ".n", length(frames.to.use), ".", width, "x", height.use, sep=""))
                                          
            if (is.null(fn.mov)) {
                fn.mov.base <- eval(fn.mov.base.exp)
                fn.mov.full <- file.path(fn.mov.dir, paste(fn.mov.base, ".", fmt, sep=""))
            } else {
                ## User provided a value for fn.mov. Check if it is a full file name include extension
                if (substr(fn.mov, nchar(fn.mov) - 3, nchar(fn.mov)) != paste(".", fmt, sep="")) {
                    fn.mov.full <- file.path(paste(fn.mov, ".", fmt, sep=""))
                } else { 
                    fn.mov.full <- fn.mov
                }
            }
            
            if (file.exists(fn.mov.full)) {
                if (fn.mov.exists == "stop") {
                    stop(paste(fn.mov.full, "exists"))
                } else if (fn.mov.exists == "overwrite") {
                    cat("  ", fn.mov.full, "exists, overwriting. \n")
                    if (!file.remove(fn.mov.full)) stop(paste("Unable to overwrite", fn.mov.full))
                } else if (fn.mov.exists == "auto.increment") {
                    i <- 0
                    while (file.exists(fn.mov.full)) {
                        i <- i + 1
                        if (i > 99) stop("Have tried 99 unique file names for the animation file, but all exist. Bailing")
                        
                        ## Look for a .nn pattern in the file name, then either update the 
                        ## auto-increment digits in the file name, or append new ones
                        auto.incr.dot <- substr(fn.mov.full, nchar(fn.mov.full) - 6, nchar(fn.mov.full) - 6)
                        auto.incr.digits <- substr(fn.mov.full, nchar(fn.mov.full) - 5, nchar(fn.mov.full) - 4)
                        if (auto.incr.dot=="." && length(grep("[^0-9]", auto.incr.digits, value = TRUE)) == 0) {
                            fn.mov.full <- paste(substr(fn.mov.full, 1, nchar(fn.mov.full) - 6), sprintf("%02d", i), ".", fmt, sep="")    
                        } else {
                            fn.mov.full <- paste(substr(fn.mov.full, 1, nchar(fn.mov.full) - 4), ".", sprintf("%02d", i), ".", fmt, sep="")    
                        }
                    }
                } else if (fn.mov.exists == "ask") {
                    ans <- readline(prompt = paste(fn.mov.full, " exists. Overwrite? y/n ", sep=""))
                    if (ans == "y") {
                        if (!file.remove(fn.mov.full)) stop(paste("Unable to overwrite", fn.mov.full))
                    } else {
                        return(invisible(NULL))
                    }
                } else {
                    stop("Unknown value for fn.mov.exists")
                }
                
            }

            fn.png <- paste(tempfile(pattern="img", tmpdir=tmp.dir), "%04d.png", sep="")
            
            ## Find ffmpeg.exe and create cmd line
            if (create.mov || show.cmd) {
                if (is.null(ffmpeg)) {
                    if (.Platform$OS.type == "windows") {
                        ffmpeg <- "ffmpeg.exe"
                    } else {
                        ffmpeg <- "ffmpeg"
                    }
                }
              
                ffmpeg.exec <- findonpath(ffmpeg)
                if (is.null(ffmpeg.exec)) {
                    cat(cw("Cant find ffmpeg. Please make sure this file is downloaded and saved either in the working directory or a directory on the PATH environment variable (e.g., c:/windows)", final.cr=T))
                    return(invisible(NULL))
                }
                if (fmt=="mov") {
                    cmd <- paste(ffmpeg.exec, " -r ", fps.use, " -i \"", fn.png, "\" -s ", width, "x", height.use, " -g 300 -pix_fmt rgb555be -vcodec qtrle -an \"", fn.mov.full, "\"", sep="")
                } else {
                    cmd <- paste(ffmpeg.exec, " -r ", fps.use, " -i \"", fn.png, "\" -s ", width, "x", height.use, " -c:v libx264 -crf 20 -pix_fmt yuv420p -preset slow -tune animation -an -f mp4 \"", fn.mov.full, "\"", sep="")
                }
            }
            
            if (prompt.continue) {
                cat("Ready to generate animation \n", sep="")
                cat("  Num frames=", length(frames.to.use), ". Duration=", round(duration.use, 1), "secs. fps=", fps.use, ". Skip=", skip.use, ". Frame size: ", 
                    width, "x", height.use, "\n", sep="")
                cat("  Temp folder: ", tmp.dir, "\n  Delete temp files: ", tmp.files.delete, "\n  Record background as PNG: ", bg2png, "\n", sep="")
                if (!bg2png && !is.null(tiff.fn)) cat(cw("  Tip: Creating frames with a background image may be faster if you set bg2png=TRUE", indent=4, exdent=4, final.cr=T))
                if (col.by.hour.of.day) cat("  Color points by the hour-of-day: TRUE\n")
                cat("  Format: ", fmt, "\n  File: ", if (create.mov) fn.mov.full else "<skipped>", "\n", sep="")
                ans <- readline(prompt = "Continue? y/n ")
                if (ans != "y") return(invisible(NULL))
            }
            
            start.time <- Sys.time()
            png(filename=fn.png, width=width, height=height.use, pointsize=png.pointsize)
            ##if (record.plot) dev.control(displaylist="enable")
            opar <- par(mar=mar.map , mgp=mgp.map)
            if (date.bar > 0) layout(matrix(1:2, ncol=1), heights=c(1-date.bar.pth, date.bar.pth))
        }
        on.exit(par(opar))
        
        ## Set up date bar parameters
        if (date.bar > 0) {
            db.axis <- dateticks(dt.frames, max.ticks=date.bar.bins + 1)
            db.axis.int.gmt <- as.numeric(db.axis[["tick.all"]])
            dt.frames.int.gmt <- as.numeric(dt.frames)
            db.xlim.int.gmt <- c(min(db.axis.int.gmt[1], dt.frames.int.gmt[1]), max(db.axis.int.gmt[length(db.axis.int.gmt)], dt.frames.int.gmt[length(dt.frames.int.gmt)]))
        }
                    
        ## Set up the colors for the active point if coloring locations by the hour of day
        if (col.by.hour.of.day) {
            col.xys.active.hod <- factor(col.hod[as.POSIXlt(dt.frames)[["hour"]] + 1])
        } 
        
        ## Make the frames
        cat("  Creating frames...\n")
        if (status) pb <- txtProgressBar(min = 0, max = length(frames.to.use), style = 3)
        make.frames.start <- Sys.time()
        
        for (frmidx in 1:length(frames.to.use)) {
            
            if (status) setTxtProgressBar(pb, frmidx)
            i <- frames.to.use[frmidx]
            
            ## If bg2png=T and this is the first pass, create a PNG file of just the plot area containing the map background, then load it to memory
            if (bg2png && is.null(bg.png)) {
            
                ## Switch to the top figure and apply map margins 
                par(mfg=c(1,1), mar=mar.map)
                
                ## Get the active figure area, width and height, in inches
                par.fin <- par("fin")
                
                ## Subtract the current (map) margins (in inches)
                par.mai <- par("mai")
                par.fin[1] <- par.fin[1] - sum(par.mai[c(2,4)])
                par.fin[2] <- par.fin[2] - sum(par.mai[c(1,3)])
                                
                ## Calculate the size of the figure area in pixels
                plot.region.pixels <- round(par.fin * dpi)
                                
                ## Switch off the current png device, will reopen it shortly
                dev.off()
                
                ## Create a new png device for just the figure area
                png(filename=fn.background.png, width=plot.region.pixels[1], height=plot.region.pixels[2], pointsize=png.pointsize)
                
                ## Create a blank plot with no margins and no axes
                par(mar=c(0,0,0,0))
                plot(NULL, xlim=xlim, ylim=ylim, axes=F, asp=1)

                ## Display the tiff
                if (!is.null(tiff.fn)) {
                    if (length(tiff.sgdf@data)==3) {
                        image(tiff.sgdf, red=1, green=2, blue=3, add=TRUE)
                    } else {
                        if (length(tiff.col)==0) stop("Something is wrong. tiff.col should not be an empty vector")
                        image(tiff.sgdf, col=tiff.col, add=TRUE)
                    }
                }

                ## Crop the GIS layers on first pass if needed
                if (!gis.layers.ready) {
                    gis.features <- shp.layers.crop(gis.features.full.extent, matrix(par("usr"), byrow=F, ncol=2))
                    gis.layers.ready <- TRUE
                }
                
                ## Put down the polygon layers the GIS layers
                for (featname in names(gis.features)[sapply(gis.features, function(x) x[["type"]]=="polygon")]) {
                    with(gis.features[[featname]], plot(sdf, lty=lty, pch=pch, cex=cex, col=col, border=border, lwd=lwd, add=TRUE))
                }
                
                ## Put down the background points
                if (!is.na(col.xys.background)) points(xys, col=col.xys.background, pch=20, cex=cex.xys.background)
                
                ## Put down the GIS point and line layers
                for (featname in names(gis.features)[sapply(gis.features, function(x) x[["type"]] %in% c("point", "line"))]) {
                    with(gis.features[[featname]], plot(sdf, lty=lty, pch=pch, cex=cex, col=col, border=border, lwd=lwd, add=TRUE))
                }

                ## Plot the legend
                if (ids.legend.bln) legend(all.ids.legend, legend=legend.str, col=ids.cols, pch=20, bg="white", title="id", cex=all.ids.legend.cex) 
                
                ## We're done with all the background elements. Close the device and read the png file back into memory as a 3-dimension matrix
                dev.off()
                bg.png <- png::readPNG(fn.background.png)
                
                ## Reopen the original png device with the original dimensions
                png(filename=fn.png, width=width, height=height.use, pointsize=png.pointsize)
                if (date.bar > 0) layout(matrix(1:2, ncol=1), heights=c(1-date.bar.pth, date.bar.pth))
                
            }
            
            ## Create a blank plot
            par(mar=mar.map, mgp=mgp.map)
            plot(NULL, xlim=xlim, ylim=ylim, xlab=if (axes.titles) axes.lbl[1] else "", ylab=if (axes.titles) axes.lbl[2] else "", xaxt=tick.show, yaxt=tick.show, cex.axis=0.8, main=str.title, asp=1)
            usr.map <- par("usr")
            
            if (bg2png) {
                ## Display the raster image with the static background elements
                rasterImage(bg.png, usr.map[1], usr.map[3], usr.map[2], usr.map[4], interpolate=FALSE)
                box()
            } else {
            
                ## Display the tiff
                if (!is.null(tiff.fn)) {
                    if (length(tiff.sgdf@data)==3) {
                        image(tiff.sgdf, red=1, green=2, blue=3, add=TRUE)
                    } else {
                        if (length(tiff.col)==0) stop("Something is wrong. tiff.col should not be an empty vector")
                        image(tiff.sgdf, col=tiff.col, add=TRUE)
                    }
                }
                
                ## Crop the GIS layers on first pass if needed
                if (!gis.layers.ready) {
                    gis.features <- shp.layers.crop(gis.features.full.extent, matrix(par("usr"), byrow=F, ncol=2))
                    gis.layers.ready <- TRUE
                }
                
                ## Put down the polygon layers the GIS layers
                for (featname in names(gis.features)[sapply(gis.features, function(x) x[["type"]]=="polygon")]) {
                    with(gis.features[[featname]], plot(sdf, lty=lty, pch=pch, cex=cex, col=col, border=border, lwd=lwd, add=TRUE))
                }
                
                ## Put down the background points
                if (!is.na(col.xys.background)) points(xys, col=col.xys.background, pch=20, cex=cex.xys.background)
                
                ## Put down the points and lines layers
                for (featname in names(gis.features)[sapply(gis.features, function(x) x[["type"]] %in% c("point", "line"))]) {
                    with(gis.features[[featname]], plot(sdf, lty=lty, pch=pch, cex=cex, col=col, border=border, lwd=lwd, add=TRUE))
                }

                ## Plot the legend
                if (ids.legend.bln) legend(all.ids.legend, legend=legend.str, col=ids.cols, pch=20, bg="white", title="id") 
            
                ## We'll come back and plot the active point later
            }
            
            ## Set up the date bar
            if (date.bar > 0) {
                par(mar=mar.db, mgp=mgp.db)
                plot(NULL, xlim=db.xlim.int.gmt, ylim=c(0,1), axes=F, xlab="", ylab="")
                usr.db <- par("usr")
                axis(side=1, at=db.axis.int.gmt, labels=format(db.axis[["tick.all"]], format=db.axis[["format.str"]]), cex.axis=cex.axis.db, col.axis=col.db, col=col.db, col.ticks=col.db)
            }
            
            ## Time to plot the active point(s). Set the focus to the top frame and apply map margins
            par(mfg=c(1,1), mar=mar.map, mgp=mgp.map, usr=usr.map)
            
            ## Plot the active point
            if (col.by.hour.of.day) {
                points(coordinates(lxy[["pts"]])[dt.frames.lxyidx.lst[[i]], ,drop=FALSE], col=if (all.ids.col.unique) ids.cols[lxy.id.int[dt.frames.lxyidx.lst[[i]]]] else as.character(col.xys.active.hod[dt.frames.lxyidx.lst[[i]]]), pch=20, cex=cex.xys.active)
            } else {
                points(coordinates(lxy[["pts"]])[dt.frames.lxyidx.lst[[i]], ,drop=FALSE], col=if (all.ids.col.unique) ids.cols[lxy.id.int[dt.frames.lxyidx.lst[[i]]]] else col.xys.active, pch=20, cex=cex.xys.active)
            }
            
            ## Plot the date-time label
            if (dt.label) {
                dt.text <- format(dt.frames[i], format="%b. %d, %Y. %H:%M")
                if (!is.na(dt.label.bg)) {
                    dt.buff <- strwidth("m", units = "user", font=2) * 0.2
                    dt.halfwidth <- dt.buff + strwidth(dt.text, units = "user", font=2) / 2 
                    dt.halfheight <- dt.buff + strheight(dt.text, units = "user", font=2)  / 2
                    dt.boxX <- dtlabel.pt[1,1] + dt.halfwidth * c(-1,1,1,-1,-1)
                    dt.boxY <- dtlabel.pt[1,2] + dt.halfheight * c(1,1,-1,-1,1)
                    polygon(dt.boxX, dt.boxY, border=NA, col=dt.label.bg)
                }
                text(dtlabel.pt, labels=dt.text, font=2, col=dt.label.col)
            }
            
            ## Put a line on the date bar
            if (date.bar > 0) {
                ## Set the focus to the bottom frame
                par(mfg=c(2,1), mar=mar.db, mgp=mgp.db, usr=usr.db)
                
                #For some funny reason it won't plot the date bar line without the following:
                box(lty=0)
                
                points(x=rep(dt.frames.int.gmt[i], 2), y=c(0,0.5), type="l", col=col.xys.active, lwd=3)
            }
            
        }
        if (status) close(pb)
        
        if (screen.test) {
            return(invisible(NULL))
        } else {
            dev.off()
            if (report.time) {
                time.taken.frames = difftime(Sys.time(), make.frames.start, units="auto")
                time.taken.frames.secs = as.numeric(difftime(Sys.time(), make.frames.start, units="secs"))
                cat("   Time to make frames:", round(time.taken.frames, 1), units(time.taken.frames), " (", round(time.taken.frames.secs / length(frames.to.use), 1), " seconds/frame) \n", sep = " ") 
            }    

            if (create.mov) {         
                if (.Platform$OS.type == "windows") {    
                    shell.ffmpeg <- system(cmd, wait=TRUE, invisible=TRUE, intern=FALSE, minimized=FALSE)
                } else {
                    shell.ffmpeg <- system(cmd, wait=TRUE, intern=FALSE)
                }
                if (shell.ffmpeg == 0) {
                    cat("  Created ", fn.mov.full, " (", round(file.info(fn.mov.full)[["size"]] / 2^20, 3), "Mb)\n", sep="")    
                } else {
                    cat("Creating the animation file appears to have failed \n")
                    return(NULL)
                }
                res[[names(idx.ids.lst)[idx.set]]] <- list(fn=fn.mov.full, dim=c(width, height.use))
            }
            if (tmp.files.delete) file.remove(sprintf(fn.png, 1:length(frames.to.use)))                      
        }
        
    }

    if (report.time && !screen.test) {
        time.taken.total = difftime(Sys.time(), start.time, units="auto")
        cat("   Total time:", round(time.taken.total, 1), units(time.taken.total), "\n", sep = " ")    
    } 
    
    if (show.cmd) {
        cat("\nffmpeg command to make video file:\n")
        cat(cmd, "\n")
    }   

    if (beep) {
        flush.console()
        for (i in 1:3) {
            alarm()
            Sys.sleep(0.8)
        }
    }

    return(invisible(res))
}

Try the tlocoh package in your browser

Any scripts or data that you put into this service are public.

tlocoh documentation built on May 2, 2019, 5:27 p.m.