R/derivate.R

In santiagoh719/ClimFunctions: Functions For Climate Data

derivate <- function(data_array, lat, lon, lat_lon = TRUE,
                     dx = TRUE, dy = TRUE, spherical_x = FALSE,
                     a = 6371000, rad = FALSE){

  checks <- makeAssertCollection()
  assert_array(data_array,mode = 'numeric',d = 2, add = checks)
  assert_numeric(lat, add = checks)
  assert_numeric(lon, add = checks)
  assert_numeric(a, len = 1, finite = TRUE,add = checks)

  assert_logical(dx, add = checks)
  assert_logical(dy, add = checks)

  assert_logical(spherical_x, add = checks)

  assert_logical(lat_lon, add = checks)
  assert_logical(rad, add = checks)
  reportAssertions(checks)

  if(lat_lon){
    if(dim(data_array)[1] != length(lat) | dim(data_array)[2] != length(lon)){
      stop('The dimetions of data_array do not match with the length of lat and lon, check if the data is order lat_lon or not')
    }
  }else{
    if(dim(data_array)[1] != length(lat) | dim(data_array)[2] != length(lon)){
      stop('The dimetions of data_array do not match with the length of lat and lon, check if the data is order lat_lon or not')
    }
    data_array <- aperm(a = data_array,perm = c(2,1))
  }

  if(!rad){
    lat <- lat / 180 * pi
    lon <- lon / 180 * pi
  }

  if(dx){

    delta_lon <- (lon[3] - lon[1]) * a
    lat_fact <- 1/cos(lat)

    delta_data_array <- cbind(NA,t(diff(t(data_array), lag = 2)),NA)
    if(spherical_x) {
      delta_data_array[,1] <- data_array[,2] -
        data_array[,ncol(data_array)]

      delta_data_array[,ncol(data_array)] <- data_array[,1] -
        data_array[,ncol(data_array) - 1]
    }

    dvar_dx <- (delta_data_array / delta_lon) * array(data = lat_fact, dim(delta_data_array))

  }else{
    dvar_dx <- NULL
  }

  if(dy){

    delta_lat <- (lat[3] - lat[1]) * a

    delta_data_array <- rbind(NA,diff(data_array, lag = 2),NA)
    dvar_dy <- (delta_data_array / delta_lat)

  }else{
    dvar_dy <- NULL
  }

  if(!lat_lon){
    if(dx){
      dvar_dx <- aperm(a= dvar_dx,perm = 2:1)
    }
    if(dy){
      dvar_dy <- aperm(a= dvar_dy,perm = 2:1)
    }
  }

  return(list(dx = dvar_dx,
              dy = dvar_dy))
}