R/age_moto.R

In vein: Vehicular Emissions Inventories

#' Returns amount of vehicles at each age
#'
#' @family age
#'
#' @description \code{\link{age_moto}} returns amount of vehicles at each age
#'
#' @param x Numeric; numerical vector of vehicles with length equal to lines features of road network
#' @param name Character; of vehicle assigned to columns of dataframe
#' @param a Numeric; parameter of survival equation
#' @param b Numeric; parameter of survival equation
#' @param agemin Integer; age of newest vehicles for that category
#' @param agemax Integer; age of oldest vehicles for that category
#' @param k Numeric; multiplication factor. If its length is > 1, it must match the length of x
#' @param bystreet Logical; when TRUE it is expecting that 'a' and 'b' are numeric vectors with length equal to x
#' @param net SpatialLinesDataFrame or Spatial Feature of "LINESTRING"
#' @param verbose Logical;  message with average age and total numer of vehicles
#' @param namerows Any vector to be change row.names. For instance, name of
#' regions or streets.
#' @param time Character to be the time units as denominator, eg "1/h"
#' @return dataframe of age distrubution of vehicles
#' @importFrom sf st_sf st_as_sf
#' @note
#' The functions age* produce distribution of the circulating fleet by age of use.
#' The order of using these functions is:
#'
#' 1. If you know the distribution of the vehicles by age of use , use:  \code{\link{my_age}}
#' 2. If you know the sales of vehicles, or the registry of new vehicles,
#' use \code{\link{age}} to apply a survival function.
#' 3. If you know the theoretical shape of the circulating fleet and you can use
#' \code{\link{age_ldv}}, \code{\link{age_hdv}} or \code{\link{age_moto}}. For instance,
#' you dont know the sales or registry of vehicles, but somehow you know
#' the shape of this curve.
#' 4. You can use/merge/transform/adapt any of these functions.
#' @export
#' @examples \dontrun{
#' data(net)
#' MOTO_E25_500 <- age_moto(x = net$ldv, name = "M_E25_500", k = 0.4)
#' plot(MOTO_E25_500)
#' MOTO_E25_500 <- age_moto(x = net$ldv, name = "M_E25_500", k = 0.4, net = net)
#' plot(MOTO_E25_500)
#' }
age_moto <- function (x,
                      name = "age",
                      a = 0.2,
                      b = 17,
                      agemin = 1,
                      agemax = 50,
                      k = 1,
                      bystreet = FALSE,
                      net,
                      verbose = FALSE,
                      namerows,
                      time){
  # check na
  x[is.na(x)] <- 0

  x <- as.numeric(x)

  # length k
  if(length(k)  > 1 & length(k) < length(x)){
    stop("length of 'k' must be 1 ore equal to length of 'x'")
  }

  # check agemax
  if(agemax < 1) stop("Agemax should be bigger than 1")


  # bystreet = TRUE
  if (bystreet == TRUE){
    if(length(x) != length(a)){
      stop("Lengths of veh and age must be the same")
    }
    d <- suca <- list()
    for (i in seq_along(x)) {
      suca[[i]] <- function (t) {1/(1 + exp(a[i]*(t+b[i])))+1/(1 + exp(a[i]*(t-b[i])))}
      anos <- seq(agemin,agemax)
      d[[i]] <- (-1)*diff(suca[[i]](anos))
      d[[i]][length(d[[i]])+1] <- d[[i]][length(d[[i]])]
      d[[i]] <- d[[i]] + (1 - sum(d[[i]]))/length(d[[i]])
      d[[i]] <- d[[i]]*x[i]
    }
    df <- as.data.frame(matrix(0,ncol=length(anos), nrow=1))
    for (i in seq_along(x)) {
      df[i,] <- d[[i]]
    }

    df <- as.data.frame(cbind(as.data.frame(matrix(0,ncol=agemin-1,
                                                   nrow=length(x))), df))

    names(df) <- paste(name,seq(1,agemax),sep="_")

    df <- df*k

    if(verbose){
      message(paste("Average age of",name, "is",
                    round(sum(seq(1,agemax)*base::colSums(df)/sum(df)), 2),
                    sep=" "))
      message(paste("Number of",name, "is",
                    round(sum(df, na.rm = T), 3),
                    " [veh]",
                    sep=" "))
      cat("\n")
    }

    if(!missing(namerows)) {
      if(length(namerows) != nrow(df)) stop("length of namerows must be the length of number of rows of veh")
      row.names(df) <- namerows
    }



    # replace NA and NaN
    df[is.na(df)] <- 0

    if(!missing(net)){
      netsf <- sf::st_as_sf(net)
      if(!missing(time)){
        dfsf <- sf::st_sf(Vehicles(df*k, time = time), geometry = sf::st_geometry(netsf))
      } else {
        dfsf <- sf::st_sf(Vehicles(df*k), geometry = sf::st_geometry(netsf))
      }
      return(dfsf)
    } else {
      if(!missing(time)){
        return(Vehicles(df*k, time = time))
      } else {
        return(Vehicles(df*k))
      }
    }

    #bystreet = FALSE
  } else {
    suca <- function (t) {1/(1 + exp(a*(t+b)))+1/(1 + exp(a*(t-b)))}
    anos <- seq(agemin,agemax)
    d <- (-1)*diff(suca(anos))
    d[length(d)+1] <- d[length(d)]
    d <- d + (1 - sum(d))/length(d)
    df <- as.data.frame(as.matrix(x) %*%matrix(d,ncol=length(anos), nrow=1))

    df <- as.data.frame(cbind(as.data.frame(matrix(0,ncol=agemin-1,
                                                   nrow=length(x))), df))

    names(df) <- paste(name,seq(1,agemax),sep="_")

    df <- df*k

    if(verbose){
      message(paste("Average age of",name, "is",
                    round(sum(seq(1,agemax)*base::colSums(df)/sum(df)), 2),
                    sep=" "))
      message(paste("Number of",name, "is",
                    round(sum(df, na.rm = T)/1000, 2),
                    "* 10^3 veh", sep=" "))
      cat("\n")
    }

    if(!missing(namerows)) {
      if(length(namerows) != nrow(df)) stop("length of namerows must be the length of number of rows of veh")
      row.names(df) <- namerows
    }

    # replace NA and NaN
    df[is.na(df)] <- 0


    if(!missing(net)){
      netsf <- sf::st_as_sf(net)
      if(!missing(time)){
        dfsf <- sf::st_sf(Vehicles(df, time = time), geometry = sf::st_geometry(netsf))
      } else {
        dfsf <- sf::st_sf(Vehicles(df), geometry = sf::st_geometry(netsf))
      }
      return(dfsf)
    } else {
      if(!missing(time)){
        return(Vehicles(df, time = time))
      } else {
        return(Vehicles(df))
      }
    }
  }
}