R/canvec.qplot.R

In rcanvec: Access and Plot CanVec and CanVec+ Data for Rapid Basemap Creation in Canada

#' Quickly Plot Canvec Data
#' 
#' Quickly plot CanVec data with options to change the plotting style of each. 
#' If data does not exist in the cache it will be downloaded. Simplest usage
#' uses \code{searchbbox()} to find an appropriate bounding box (e.g.
#'  \code{canvec.qplot(bbox=searchbbox("Wolfville NS"))}). Be careful plotting
#'  feature-intensive layers (e.g. "road", "building") over large areas (e.g.
#'  \code{searchbbox("toronto, on")}). This will happily run but plotting the map
#'  may take up to 20 minutes!
#' 
#' @param ntsid One or more NTS References as generated by \code{nts()}
#' @param bbox A bounding box describing the desired extent. If no \code{ntsid} is
#' provided, nts(bbox=bbox) will be invoked to find the appropriate NTS Reference(s)
#' @param layers A list of layers as defined in \code{canvec_layers$id} in the order in
#' which they should be plotted
#' @param options A list object containing the options for each layer in the form
#' \code{options$layerid <- list(col="lightblue")}
#' @param cachedir Pass a specific cache directory in which files have been extracted.
#'                  Default value is that returned by \code{canvec.cachedir()}
#' @param data A list object that contains the loaded Spatial* data to be plotted.
#' This should always be an object that was returned by \code{canvec.qplot()}
#' @param plotdata TRUE if data should be plotted, FALSE if data should just be loaded.
#' @param atscale One of \code{nts.SCALE50K} (CanVec data) or \code{nts.SCALE250K} (CanVec+ data)
#' @param stoponlargerequest Stop if a large (greater than 4 tiles) area is requested. Defaults to
#' \code{TRUE}.
#' @param epsg The epsg code in which to plot the data, or \code{NULL} for automatic. Use
#' \code{epsg=3857} to layer on Open Street Map tiles, or \code{epsg=269XX} (where \code{XX}
#' is the UTM Zone) for a UTM projection. Defaults to no projection, although \code{sp::plot}
#' adjusts the aspect such that the default does not appear distorted.
#' @param ... A list of graphical parameters passed to the inital call to \code{plot()}. Use
#' \code{add=TRUE} to layer on an existing plot.
#' @return A list object that contains the Spatial* data that was plotted
#' 
#' @examples
#' \donttest{
#'  #simplest use using searchbbox() from {prettymapr}
#'  library(prettymapr)
#'  wolfville <- searchbbox("Wolfville NS", source="google")
#'  canvec.qplot(bbox=wolfville)
#'  canvec.qplot(bbox=wolfville, layers=c("waterbody", "forest"))
#' 
#'  #can also plot by NTS sheet and use bbox= or xlim, ylim to zoom.
#'  canvec.qplot(nts("21h1"), layers=c("waterbody", "forest", "contour", "river", "road"))
#'  canvec.qplot(bbox=makebbox(45.1, -64.35, 45.05, -64.4), 
#'         layers=c("waterbody", "contour", "river", "building", "building_poly", "road"))
#'
#'  #method returns plot data argument so data does not need to be loaded each time. 
#'  #this will not work when changing nts sheets.
#'  plotdata <- canvec.qplot(nts("21h1"), layers=c("waterbody", "forest", "contour", "river"))
#'  plotdata <- canvec.qplot(bbox=makebbox(45.1, -64.35, 45.05, -64.4), 
#'                         layers=c("waterbody", "contour", "river"),
#'                         data=plotdata)
#'  
#'  #use with prettymapr::addscalebar() and prettymapr::addnortharrow()
#'  library(prettymapr)
#'  wolfville <- searchbbox("Wolfville NS", source="google")
#'  canvec.qplot(bbox=wolfville)
#'  addscalebar()
#'  addnortharrow()
#'  
#'  #or use with prettymapr::prettymap() to set margins and add north arrow/
#'  #scalebar
#'  prettymap(canvec.qplot(bbox=wolfville))
#'  
#'  }
#' 
#' @export
#' 
canvec.qplot <- function(ntsid=NULL, bbox=NULL, 
                          layers=c("waterbody", "forest", "contour", "river", "road"), 
                          options=NULL, data=NULL, cachedir=NULL, plotdata=TRUE, atscale=nts.SCALE50K, 
                          stoponlargerequest=TRUE, epsg=NULL, ...) {
  
  if(!is.null(epsg)) {
    crs <- sp::CRS(paste0("+init=epsg:", epsg))
  } else {
    crs <- NULL
  }
  
  if(!is.null(bbox)) {
    if(is.null(ntsid)) {
      ntsid <- nts(bbox=bbox, atscale = atscale)
      if(class(ntsid)=="list" && length(ntsid)>4 && stoponlargerequest) 
        stop("Current requrest requires loading of ", length(ntsid),
             " mapsheets. This may download a large amount of data and/or take a",
             " very long time to load. Rerun with stoponlargerequest=FALSE to continue.")
    }
    if(!is.null(crs)) {
      #transform bbox
      coords <- sp::coordinates(
        sp::spTransform(
          sp::SpatialPoints(sp::coordinates(t(bbox)), sp::CRS("+init=epsg:4326")), crs))
      bbox <- t(coords)
    }
    xlim <- bbox[1,]
    ylim <- bbox[2,]
  } else if(is.null(ntsid)) {
    stop("No arguments specified for data to plot")
  }
  
  if(is.null(cachedir)) {
    cachedir <- canvec.cachedir()
  }
  
  if(class(ntsid) != "list") {
    ntsid <- list(ntsid)
    #check to make sure arguments passed are indeed NTS refs
    for(singleid in ntsid) {
      if((class(singleid) != "character") || 
         (length(singleid)<2) || 
         (length(singleid)>3)) stop("One or more arguments specified was not an NTS reference: ",
                                    singleid, ". Did you mean bbox=... or nts(...) instead?")
    }
  }
  
  if(is.null(data)) {
    #download
    canvec.download(ntsid, cachedir=cachedir)
    data <- list()
  }
  
  if(is.null(options)) {
    options <- list()
  }
  
  #load data
  for(layerid in layers) {
    if(is.null(data[[layerid]]))
      data[[layerid]] <- canvec.load(ntsid, layerid)
  }
  
  if(plotdata) {
    #plot corners of mapsheet extents 
    if(class(ntsid)=="list") {
      if(length(ntsid)==0) stop("Cannot plot background for zero mapsheets")
      bbox1 <- nts.bbox(ntsid[[1]])
      for(singleid in ntsid) {
        bbox2 <- nts.bbox(singleid)
        bbox1 <- matrix(c(min(bbox1[1,1], bbox2[1,1]), min(bbox1[2,1], bbox2[2,1]),
                          max(bbox1[1,2], bbox2[1,2]), max(bbox1[2,2], bbox2[2,2])), ncol=2, byrow=FALSE)
      }
    } else {
      bbox1 = nts.bbox(ntsid)
    }
    
    if(!is.null(crs)) {
      coords <- sp::coordinates(t(bbox1))
      spoints <- sp::SpatialPoints(coords, proj4string = sp::CRS("+proj=longlat +ellps=GRS80 +no_defs"))
      spoints <- sp::spTransform(spoints, crs)
      coords <- sp::coordinates(spoints)
      bbox1 <- t(coords)
    } else {
      coords <- sp::coordinates(t(bbox1))
      spoints <- sp::SpatialPoints(coords, proj4string = sp::CRS("+proj=longlat +ellps=GRS80 +no_defs"))
    }
    
    plotargs <- list(...)
    if(is.null(bbox)) {
      if(is.null(plotargs$xlim))
        xlim <- bbox1[1,]
      if(is.null(plotargs$ylim))
        ylim <- bbox1[2,]
    }
    
    sp::plot(spoints, pch=".", xlim=xlim, ylim=ylim, ...)
    
    #plot data
    for(layerid in layers) {
      layeropts <- options[[layerid]]
      if(is.null(layeropts))
        layeropts <- canvec.defaultoptions(layerid)
      canvec.plot(data[[layerid]], layeropts, crs=crs, add=TRUE)
    }
    
  }

  invisible(data)
}