R/plots_auxiliar.R

In roliveros-ramos/kali: Tools for habitat and niche modeling

# create grid axes (lat, lon)
createGridAxes = function(lat,lon, dx, dy=dx, center=FALSE) {
  # Create a rectangular grid given lat, lon and dxy.
  # No correction by Earth curvature
  if(dx <= 0 || dy <= 0) stop("dx and dy must be positive.")
  
  if(isTRUE(center)) {
    lat[which.min(lat)] = lat[which.min(lat)] + 0.5*dy
    lat[which.max(lat)] = lat[which.max(lat)] - 0.5*dy
    lon[which.min(lon)] = lon[which.min(lon)] + 0.5*dx
    lon[which.max(lon)] = lon[which.max(lon)] - 0.5*dx
  }
  
  nx = 1 + round(diff(lon)/dx, 0)
  ny = 1 + round(diff(lat)/dy, 0)

  lats.rho = seq(from=min(lat),to=max(lat), length=ny)
  lons.rho = seq(from=min(lon),to=max(lon), length=nx)
  
  lats.psi = seq(from=min(lat) - 0.5*dy, to=max(lat) + 0.5*dy, length=ny+1)
  lons.psi = seq(from=min(lon) - 0.5*dx, to=max(lon) + 0.5*dx, length=nx+1)
  
  rho = list(lat=lats.rho, lon=lons.rho)
  psi = list(lat=lats.psi, lon=lons.psi)
  
  nlat = length(rho$lat)
  nlon = length(rho$lon)
  
  LAT = matrix(rho$lat, ncol=nlat, nrow=nlon, byrow=TRUE)
  LON = matrix(rho$lon, ncol=nlat, nrow=nlon)
  
  area = (111*dy)*(111*cos(LAT*pi/180)*dx)
  
  output = list(lon=lons.rho, lat=lats.rho, rho=rho, psi=psi, 
                LON=LON, LAT=LAT, area=area)
  
  return(output)
}

# create time axis
createTimeAxis = function(start, end, frequency=12, center=FALSE, shift=TRUE,
                          day=1) {
  
  times  = seq(from=start[1] + (start[2]-1)/frequency, 
               to= end[1] + end[2]/frequency, by=1/frequency)
  if(center) {
    out = list(bounds = times, center = times[-length(times)] + 0.5/frequency)
    out$year = floor(out$center)
    out$month = round(frequency*(out$center-out$year) + 0.5, 0)
    out$season = out$month
    seasons = matrix(if(shift) 1:12 else c(12, 1:11), ncol=3, byrow=TRUE)
    out$season[out$season %in% seasons[1,]] = "summer"
    out$season[out$season %in% seasons[2,]] = "fall"
    out$season[out$season %in% seasons[3,]] = "winter"
    out$season[out$season %in% seasons[4,]] = "spring"
    out$dates = as.Date(paste(out$year, out$month, 1, sep="-"))
    out$start = out$year + 
      (lubridate::yday(out$dates)-1)/(365+lubridate::leap_year(out$year))
    times = out
  } 
  return(times)
}

# rename map.axes an overwrite map's map.axes?
map.axes2 = function(sides=c(1,2), cex.axis=0.75, line=-0.4) {
  
  .axis.map(sides[1], "lon", las=1, cex.axis=cex.axis, line=line, tick=FALSE)
  .axis.map(sides[2], "lat", las=1, cex.axis=cex.axis, line=line, tick=FALSE)
  axis(1, labels=FALSE)
  axis(2, labels=FALSE)
  axis(3, labels=FALSE)
  axis(4, labels=FALSE)
  box()
  
  return(invisible(NULL))
}


# this function change log coordinates from -180,+180 to 0-360
# change name
rotateMap =  function(x) {
  out = x$lon>0
  out[is.na(out)] = FALSE
  x$lon[out] = x$lon[out] - 360
  return(x)
}


# change the center of a map database
map2 =  function(database, center, ...) {
  # change the center of the map (e.g. to show complete Pacific Ocean)
  # From stackoverflow (check reference and how to cite)
  Obj = maps::map(database, ..., plot=FALSE)
  coord = cbind(Obj[[1]],Obj[[2]])
  # split up the coordinates
  id = rle(!is.na(coord[,1]))
  id = matrix(c(1,cumsum(id$lengths)),ncol=2,byrow=TRUE)
  polygons = apply(id, 1, function(i) {coord[i[1]:i[2],]} )
  # split up polygons that differ too much
  polygons = lapply(polygons, function(x) {
    x[,1] = x[,1] + center
    x[,1] = ifelse(x[,1]>180,x[,1]-360,x[,1])
    if(sum(diff(x[,1])>300,na.rm=TRUE) >0){
      id = x[,1] < 0
      x = rbind(x[id,],c(NA,NA),x[!id,])
    }
    
    return(x)
  })
  
  # reconstruct the object
  polygons = do.call(rbind,polygons)
  Obj[[1]] = c(polygons[,1], c(NA, NA), polygons[,1] + 360) 
  Obj[[2]] = c(polygons[,2], c(NA, NA), polygons[,2])
  
  maps::map(Obj,...)
}
 
# create method fill??
fill.map = function(x, mask=1, FUN="max") {
  
  # fill a map with the max (min, mean, median) value of the map
  
  FUN = match.fun(FUN)
  maximo = FUN(x, na.rm=TRUE)
  x[is.na(x)] = maximo
  x = x*mask
  
  return(x)
}

# create method convexHull, this will be convexHull.array
convexHull = function(x, mask) {
  
  dim.old = dim(x)
  out = apply(x, 3, convexHull.map, mask=mask)  
  dim(out) = dim.old
  
  return(out)
}


convexHull.map = function(x, mask, excludeZeros=TRUE) {
  
  if(all(is.na(x))) {
    return(x)
  } else {
    ind = if(excludeZeros) (!!x & !is.na(x)) else !is.na(x) 
    xc = col(x)[ind]
    xr = row(x)[ind]
    coord = cbind(xr, xc)
    xvar = x[coord]
    
    x.int = interp(x=xr, y=xc, z=xvar, 
                   xo=seq(len=nrow(x)), yo=seq(len=ncol(x)),
                   duplicate="strip", linear=TRUE,
                   extrap=FALSE)$z * mask
    x.int[mask==0] = NA
    
    return(x.int)    
  }
}

# change resolution of a map
regridMap = function(old, new, normalize=TRUE, ...) {
  
  # generalize to any kind of input and output data
  
  new$mask[new$mask==0] = NA # correct mask
  
  ind = !apply(old, 1, function(x) any(is.na(x)))
  old = old[ind, ] # remove NAs
  
  newp   = interpp(old$lon, old$lat, old$z, 
                   xo=as.numeric(new$lon), 
                   yo=as.numeric(new$lat),
                   ...)$z
  
  newmap = matrix(newp, 
                  ncol=ncol(new$mask), 
                  nrow=nrow(new$mask))
  
  newmap[is.na(newmap)] = 0
  newmap = newmap*new$mask
  if(normalize) newmap = newmap/sum(newmap, na.rm=TRUE)
  
  return(newmap)
}