zoom <- function(x, y, data, t=1,..., same = FALSE, mark = "+") {

# Zoom in on either a two-dimensional map (matrix) with p longitude and q 
# latitude space dimensions or a three-dimensional array (field) with 
# p longitude and q latitude as the first two space dimensions and n as the
# third time dimension, by selecting two corners that define a rectangular
# region of the space and a subset of time slices 
#
# Description:
#
#      Returns a subset of a two-dimensional map (matrix) or a
#      subset of a three-dimensional array (field) and re-plot data
#      of the first time slice of a rectangular region selected by
#      clicking on an image.
#
# Usage:
#
#      zoom(x, y, data, t=1, ..., same = FALSE, mark = "+")
#
# Arguments:
#
#      x: Vector of longitudes 
#
#      y: Vector of latitudes
#
#      data: two-dimensional map (matrix) or three-dimensional array of data.
#            Space dimensions of `data' must be `length(x)' by `length(y)'
#            respectively. 
#            Three-dimensional array must have the following dimention order:
#            first dimension for longitude, second dimension for latitude,
#            and third dimension for time
#
#      t: Subset of times slices to be extracted. Default is 1 (i.e. extract
#         only the first time slice
#
#    ...: Additional arguments passed to `image'
#
#   same: Logical: if `FALSE' (the default), zoomed images are re-plotted
#         in a new window; if `TRUE', each image is plotted in the same
#         window without cleaning the longitude and latitude coordinates
#         previously plotted  
#
#   mark: A vector of numbers or strings with which selected cells are
#         labelled. If shorter than the number of grid cells, values are
#         recycled; if `NULL', selections are not marked
#
# Details:
#
#      Clicking on two grid cells in an image plot created with grid
#      lines at locations `x' and `y' causes the image of the first time slice
#      to be re-plotted for the rectangular region thus defined.
#
# Output:
#
#      A list with three components is returned invisibly:
#
#         x: Horizontal (lon) coordinates of grid lines used to plot the final
#            zoomed image
#
#         y: Vertical (lat) coordinates of grid lines used to plot the final
#            zoomed image
#
#         data: Matrix or three-dimensional array containing the selected data
#
#         t: Extracted times slices 
#
# See Also:
#
#      `identify', `image'
#
# Author:
#
#      Chris Ferro <c.a.t.ferro@reading.ac.uk> 27 October 2005
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
#
# Examples:
#
#      # Zoom in three-dimensional array 
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 100)
#      data <- array(rnorm(prod(dim)), dim)
#      output <- zoom(x, y,data,t=c(1:3))
#      output <- zoom(x, y,data,t=c(2,4,6))
#      output <- zoom(x, y,data,t=c(10,20,30,40))
#
#      # Zoom in matrix
#      matr <- data[,,1]
#      output <- zoom(x, y,matr)


if (length(dim(data))==2) {
  xnew <- x
  ynew <- y 
  image(x, y, data)
  if (min(x)<0){map("world",add=TRUE)}
  else{map("world2",add=TRUE)}
  cat("Left-click to select two corners of zoom region; right-click to finish.\n")
  repeat {
    grid <- expand.grid(y = ynew, x = xnew)[, 2:1]
    if(is.null(mark)) lab <- 1:nrow(grid)
    else lab <- rep(mark, length.out = nrow(grid))
    pts <- identify(grid, labels = lab, n = 2, plot = !is.null(mark), offset = 0)
    if(length(pts) < 2) break
    loc <- grid[pts, ]
    xnew <- x[x >= min(loc$x) & x <= max(loc$x)]
    ynew <- y[y >= min(loc$y) & y <= max(loc$y)]
    datanew <- data[match(xnew, x), match(ynew, y), drop = FALSE]
    if(!same) x11()
    par(mfg = par("mfg"))
    image(xnew, ynew, datanew, ...)
    if (min(x)<0){map("world",add=TRUE)}
    else{map("world2",add=TRUE)}    
  }
}  

if(length(dim(data))==3) {
  z <- data[,,1]
  xnew <- x
  ynew <- y
  tnew <- t
  image(x, y, data[,,1])
  if (min(x)<0){map("world",add=TRUE)}
  else{map("world2",add=TRUE)}
  cat("Left-click to select two corners of zoom region; right-click to finish.\n")
  repeat {
    grid <- expand.grid(y = ynew, x = xnew)[, 2:1]
    if(is.null(mark)) lab <- 1:nrow(grid)
    else lab <- rep(mark, length.out = nrow(grid))
    pts <- identify(grid, labels = lab, n = 2, plot = !is.null(mark), offset = 0)
    if(length(pts) < 2) break
    loc <- grid[pts, ]
    xnew <- x[x >= min(loc$x) & x <= max(loc$x)]
    ynew <- y[y >= min(loc$y) & y <= max(loc$y)]
    znew <- z[match(xnew, x), match(ynew, y), drop = FALSE]
    datanew <- data[match(xnew, x), match(ynew, y), t, drop = FALSE]
    if(!same) x11()
    par(mfg = par("mfg"))
    image(xnew, ynew, znew, ...)
    if (min(x)<0){map("world",add=TRUE)}
    else{map("world2",add=TRUE)}
  } 
} 

invisible(list(x = xnew, y = ynew, data = datanew, t = t))
}