identifyfield <- function(x, y, data, fun=NULL, ..., mark = "*") {

# Identify points by clicking on an image plot of the first time slice of
# a given three-dimensional data array (i.e. a field) with first two space 
# dimensions (first dimension longitude, and second dimension latitude) and 
# third time dimension, or an image plot of a two-dimensional
# data matrix (i.e. a map), with first dimension longitude, and second dimension
# latitude. Also apply a user specified function to each of 
# the selected points of the three-dimensional data array.
# If a two-dimensional data matrix (i.e. a map) is provided, there is no need
# to specify a function to apply at each selected point, and
# the function returns the data values for the points selected by clicking on 
# the image plot. 
#
# Description:
#
#      Applies a function to data associated with cells selected by
#      clicking on an image plot if data is a three-dimentional array.
#      If data is a matrix (i.e. a map), the function returns the data 
#      values for the points selected by clicking on the image plot. 
#      Also returns the location of the selected points.
#
# Usage:
#
#      identifyfield(x, y, data, fun, ..., mark = "*")
#
# Arguments:
#
#      x: Longitude vector to plot the image.
#
#      y: Latitude vector to plot the image.
#
#   data: A three-dimensional array (i.e. a field) or matrix (i.e. a map)
#         containing data associated with each cell of the image. The first
#         two dimensions must be `length(x)' and `length(y)' respectively.
#
#    fun: The function to be applied to the three-dimensional array. Its 
#         first argument must be the data vector associated with a cell.
#
#    ...: Additional arguments passed to `fun'.
#
#   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 an image plot created with grid lines at locations `x'
#      and `y' causes `fun' to be called repeatedly with each of the
#      vectors of `data' corresponding to selected grid cells. Returned
#      values from `fun' are stored and returned invisibly.
#
# Value:
#
#      A list with an element for each grid-cell selection is returned
#      invisibly. Each element is a list with two or three components:
#
#       loc: The (longitude,latiture)-coordinates of the selected grid cell.
#
#       pos: The positions of the (x,y)-coordinates in `x' and `y'.
#
#      data: The vector of `data' associated with that cell.
#
#     value: The object (if any) returned by `fun' for the three-dimensional
#            array data ata each selected point.
#
# See Also:
#
#      `apply', `identify', `image'
#
# Author:
#
#      Chris Ferro <c.a.t.ferro@reading.ac.uk> 27 October 2005
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk>
#
# Examples:
#
#      # Three-dimensional array example (field)
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 100)
#      data <- array(rnorm(prod(dim)), dim)      
#      fun <- function(x, ...) { x11(); plot(x, ...); mean(x) } # plot time
#             series and compute mean value of this time series
#      out <- identifyfield(x, y, data, fun, type = "l", mark = 1:prod(dim[1:2]))
#      fun <- function(x, index, ...) { x11(); plot(index, x, ...) }
#      out <- identifyfield(x, y, data, fun, index = rexp(100))
#
#      # Matrix example (map)
#      out <- identifyfield(x, y, data[,,1])

  if (length(dim(data))==3) {
  fun <- match.fun(fun)
  image(x,y,data[,,1])
    if (min(x)<0){map("world",add=TRUE)}
    else{map("world2",add=TRUE)}
  cat("Left-click to select grid cells; right-click to finish.\n")
  grid <- expand.grid(y = y, x = x)[, 2:1]
  if(is.null(mark)) lab <- 1:nrow(grid)
  else lab <- rep(mark, length.out = nrow(grid))
  pts <- identify(grid, labels = lab, plot = !is.null(mark), offset = 0)
  series <- NULL
  loc <- grid[pts, ]
  ix <- match(loc[, 1], x)
  iy <- match(loc[, 2], y)
  #fun <- match.fun(fun)
  for(i in 1:length(pts)) {
    series[[i]] <- list(loc = loc[i, ])
    series[[i]]$pos <- c(x = ix[i], y = iy[i])
    series[[i]]$data <- data[ix[i], iy[i], ]
    series[[i]]$value <- fun(series[[i]]$data, ...)
  }
  }
  if (length(dim(data))==2) {
    image(x,y,data)
    if (min(x)<0){map("world",add=TRUE)}
    else{map("world2",add=TRUE)}
  cat("Left-click to select grid cells; right-click to finish.\n")
  grid <- expand.grid(y = y, x = x)[, 2:1]
  if(is.null(mark)) lab <- 1:nrow(grid)
  else lab <- rep(mark, length.out = nrow(grid))
  pts <- identify(grid, labels = lab, plot = !is.null(mark), offset = 0)
  series <- NULL
  loc <- grid[pts, ]
  ix <- match(loc[, 1], x)
  iy <- match(loc[, 2], y)  
  for(i in 1:length(pts)) {
    series[[i]] <- list(loc = loc[i, ])
    series[[i]]$pos <- c(x = ix[i], y = iy[i])
    series[[i]]$data <- data[ix[i], iy[i]]    
  }
  }
  invisible(series)
}