R/Region.R

In DSsim: Distance Sampling Simulations

#' @title  Class "Region" 
#' 
#' @description Class \code{"Region"} is an S4 class containing descriptions of the 
#' study area. The polygons describing the region are found in the 
#' coords slot and any gaps are described as polygons in the gaps slot.
#'
#' @name Region-class
#' @title S4 Class "Region"
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{make.region(region.name = "region.name", shapefile = region.shapefile)} 
#' @slot region.name Object of class \code{"character"}; giving the 
#'  name of the region.
#' @slot strata.name Object of class \code{"character"}; character 
#'  vector giving the names of the strata.
#' @slot units Object of class \code{"character"}; character describing
#'  the coordinate units ("km" or "m")
#' @slot area Object of class \code{"numeric"}; the area of the survey
#'  region
#' @slot box Object of class \code{"numeric"}; 4 values giving the x and
#'  y ranges of the region
#' @slot coords Object of class \code{"list"}; this list contains an
#'  element for each strata. Each of these list elements contains a list of 
#'  polygons defining the region.
#' @slot gaps Object of class \code{"list"};this list contains an
#'  element for each strata. Each of these list elements contains a list of 
#'  gaps in the region
#' @section Methods:
#' \describe{
#'  \item{\code{get.area}}{\code{signature(obj = "Region")}: retrieves the area 
#'  element }
#'  \item{\code{plot}}{\code{signature(x = "Region", y = "missing")}: plots
#'  the survey region defined by the object.}
#' }
#' @keywords classes
#' @seealso \code{\link{make.region}}
setClass(Class = "Region", 
         representation(region.name = "character", 
                        strata.name = "character",
                        units = "character", 
                        area = "numeric", 
                        box = "numeric", 
                        coords = "list", 
                        gaps = "list")
)

setMethod(
  f="initialize",
  signature="Region",
  definition=function(.Object, region.name = character(0), strata.name = character(0), units, area, shapefile = NULL, coords, gaps, check.LinkID){
    #Input pre-processing
    boundbox <- numeric(0)
    if(!is.null(shapefile)){
      # if no coordinates have been supplied then it uses the shapefile
      # put in to compensate for funny shapefiles generated by Distance
      shapefile <- check.shapefile(shapefile)
      # Now can extract coordinates
      polygons <- coords.from.shapefile(shapefile)
      coords <- polygons$coords
      gaps <- polygons$gaps  
      # Check the LinkID order
      if(check.LinkID){
        new.order <- check.LinkID.order(shapefile)
        if(!is.null(new.order)){
          #If there is a new order do the re-ordering
          warning("The LinkID values were not in sequential order in the shapefile attribute table. DSsim is reordering the strata to match that which Distance for Windows uses. This is important if your survey shapefiles were created in Distance. If you would not like them to be re-ordered please make the check.LinkID FALSE. See ?make.region for more information.", immediate. = TRUE, call. = FALSE)
          new.coords <- new.gaps <- list()
          for(i in seq(along = new.order)){
            index <- which(new.order == i)
            new.coords[[i]] <- coords[[index]]
            new.gaps[[i]] <- gaps[[index]]
          }
          coords <- new.coords
          gaps <- new.gaps
        }
      }
    }else if(length(coords) == 0 & is.null(shapefile)){
      #complains if neither the coordinates or the shapefile are supplied
      stop("You must provide either coordinates or a shapefile.", call. = FALSE)
    }
    #Gets the minimum bounding box
    boundbox <- get.bound.box(coords)
    #calculates the strata areas
    if(length(area) == 0){
      area <- calc.area(coords, gaps)
    }
    #Set slots
    .Object@region.name <- region.name
    .Object@strata.name <- strata.name
    .Object@units       <- units
    .Object@area        <- area
    .Object@box         <- boundbox
    .Object@coords      <- coords
    .Object@gaps        <- gaps
    #Check object is valid
    valid <- try(validObject(.Object), silent = TRUE)
    if(class(valid) == "try-error"){
      stop(attr(valid, "condition")$message, call. = FALSE)
    }
    # return object
    return(.Object) 
  }
)
setValidity("Region",
  function(object){
    if(length(object@strata.name) > 0){
      strata.name <- object@strata.name
      # Check that none are Total
      if(any(strata.name == "Total")){
        return("'Total' is not an accepted strata name, please ammend it.")
      }
      #check that they are all unique
      for(i in seq(along = strata.name)){
        temp <- strata.name[-i]
        for(j in seq(along = temp)){
          if(strata.name[i] == temp[j]){
            return("Stratum names must be unique")
          }
        }
      }
    }
    if(length(which(object@area <= 0)) > 0){
      return("All areas must be greater than 0")
    }    
    if(length(object@coords) != length(object@gaps)){
      return("The lengths of the coords and gaps lists differ, these must be the same and equal to the number of strata.")
    }
    if(length(object@coords) > 1 & length(object@coords) != length(object@strata.name)){
      return("Number of stratum names differs to number of strata in the shapefile.")
    }
    return(TRUE)
  }
)

# GENERIC METHODS DEFINITIONS --------------------------------------------

#' Returns the area of the region
#' 
#' @param object object of class \code{Region}
#' @return numeric value specifying the area of the region
#' @rdname get.area-methods 
setGeneric(name = "get.area", def = function(object){standardGeneric ("get.area")})

#' @rdname get.area-methods
setMethod(
  f="get.area",
  signature="Region",
  definition=function(object){
    return(object@area)
  }
)

#' Plot
#' 
#' Plots an S4 object of class 'Region'
#' 
#' @param x object of class Region
#' @param y not used
#' @param add logical indicating whether it should be added to 
#'  existing plot
#' @param plot.units allows for units to be converted between m
#'  and km
#' @param region.col fill colour for the region
#' @param gap.col fill colour for the gaps
#' @param main character plot title
#' @param ... other general plot parameters 
#' @rdname plot.Region-methods
#' @importFrom graphics polygon plot axTicks axis
#' @exportMethod plot
setMethod(
  f="plot",
  signature="Region",
  definition=function(x, y, add = FALSE, plot.units = character(0), region.col = NULL, gap.col = NULL, main = "", ...){
    # If main is not supplied then take it from the object
    if(main == ""){
      main <- x@region.name
    }
    plot.list <- function(list.coords, border = 1, fill.col = NULL){
      #lapply(list.coords, FUN = lines, type = type, col = col)
      lapply(list.coords, FUN = polygon, border = border, col = fill.col)
      invisible(list.coords)                          
    }
    #Set up plot
    if(length(plot.units) == 0){
      plot.units <- x@units
    }
    if(!add){      
      xlabel <- paste("X-coords (",plot.units[1],")", sep = "")
      ylabel <- paste("Y-coords (",plot.units[1],")", sep = "")
      plot(c(x@box[["xmin"]], x@box[["xmax"]]), c(x@box[["ymin"]], x@box[["ymax"]]), col = "white", xlab = xlabel, ylab = ylabel, main = main, yaxt = "n", xaxt = "n", ...)
      xticks <- axTicks(1)
      yticks <- axTicks(2)
      #Set up axes
      if(plot.units != x@units){
        #convert units
        if(x@units == "m" & plot.units == "km"){ 
          axis(1, at = xticks, labels = xticks/1000)
          axis(2, at = yticks, labels = yticks/1000)
        }else if(x@units == "km" & plot.units == "m"){
          axis(1, at = xticks, labels = xticks*1000)
          axis(2, at = yticks, labels = yticks*1000)
        }else{
          warning("The requested conversion of units is not currently supported, this option will be ignored.", call. = FALSE, immediate. = TRUE)
        }
      }else{
        #no unit conversion needed
        axis(1, at = xticks, labels = xticks)
        axis(2, at = yticks, labels = yticks)
      }
    }
    lapply(x@coords, FUN = plot.list, fill.col = region.col)
    lapply(x@gaps, FUN = plot.list, fill.col = gap.col)
    invisible(x)
  }    
)