.packageName <- "GRASS"
# Copyright 2001 by Roger S. Bivand
#
# contour.G is a wrapper 
contour.G <- function (x, layer = NULL, xlab = "", ylab = "",
    reverse = NULL, add = FALSE, ...) 
{
    G <- x
    if (class(G) != "grassmeta") 
        stop("Data not a grass object")
    if (!add) {
        plot(G$xlim, G$ylim, xlim = G$xlim, ylim = G$ylim, asp = 1, 
            xlab = xlab, ylab = ylab, type = "n")
    }
    if (!is.null(layer)) {
        if (length(layer) != G$Ncells) 
            stop("GRASS object metadata do not match layer length")
        if (is.null(reverse)) 
            reverse <- reverse(G)
        contour(x = G$xseq, y = G$yseq, z = t(matrix(layer[reverse], 
            nrow = G$Nrow, ncol = G$Ncol, byrow = TRUE)), add = TRUE, 
            ...)
    }
}

# Copyright 2001-4 by Roger S. Bivand
#
# east() is an access function to return the eastings coordinates
# of raster cell centres from a grassmeta object
#

east <- function(object) UseMethod("east")

east.default <- function(object) stop("no default method for east")

east.grassmeta <- function(object)
{
    if (class(object) != "grassmeta") stop("No GRASS metadata object")
    if(is.loaded("eastG", PACKAGE="GRASS")) {
	east <- .Call("eastG", object, PACKAGE="GRASS")
    } else {
        east <- as.numeric(c(matrix(object$xseq, length(object$xseq), 
	    length(object$ryseq))))
    }
    invisible(east)
}


# Copyright 1999-2003 by Roger S. Bivand
#
# gmeta is a function that returns GRASS LOCATION metadata in a
# grassmeta object.
#
# CHANGED 000329 from g.region -g to g.region -p to get numbers of
# rows and columns directly
# CHANGED 000614 to replace dataframe xy by list xy
#
# CHANGED 000628 to introduce support for compiled "gmeta"
#
# CHANGED 000706 to move east, north, obsno and reverse to access functions
#
# CHANGED 000714 to add interp argument
#
# CHANGED 030109 to cope with interpreted problems when location
# is lat/lon

gmeta <- function(interp=FALSE) {
    if (length(Sys.getenv("GISBASE")) == 0) {
	stop("No GRASS environment detected - start GRASS before entering R")
    }
    if(is.loaded("gmeta", PACKAGE="GRASS") && (interp == FALSE)) {
	G <- .Call("gmeta", PACKAGE="GRASS")
    } else {
	G <- vector(mode="list")
	G$LOCATION <- system("g.gisenv LOCATION_NAME", intern=TRUE)
	G$MAPSET <- system("g.gisenv MAPSET", intern=TRUE)
	META <- system("g.region -p", intern=TRUE)
	ML <- length(META)
	if (ML > 8) G$proj <- paste(META[1:(ML-8)], sep=" ", collapse="; ")
	if (length(grep("Latitude-Longitude", G$proj)) == 0) {
	    META.i <- function(META, i) {
		as.numeric(unlist(strsplit(META[ML-i], ":"))[2])
	    }
	    G$n <- META.i(META, 7)
	    G$s <- META.i(META, 6)
	    G$w <- META.i(META, 5)
	    G$e <- META.i(META, 4)
	    G$nsres <- META.i(META, 3)
	    G$ewres <- META.i(META, 2)
	} else {
	    dmstodd <- function(str) {
		res <- vector(mode="list")
		res$dms <- unlist(strsplit(str, ":"))
		res$n <- length(res$dms)
		res <- getNSWE(res)
		ndms <- as.numeric(res$dms)
		dd <- 0
		for (i in 1:res$n) dd <- dd + ndms[i]*(60^(-(i-1)))
		if (!is.na(res$res) && (res$res == "S" || res$res == "W"))
		    dd <- -dd
		dd
	    }
	    getNSWE <- function(dms) {
		dum <- as.character(NA)
		if(length(grep("[NSWE]", dms$dms[dms$n])) > 0) {
		    dum <- substr(dms$dms[dms$n], nchar(dms$dms[dms$n]),
			nchar(dms$dms[dms$n]))
		    dms$dms[dms$n] <- substr(dms$dms[dms$n], 1,
			nchar(dms$dms[dms$n])-1)
		}
		dms$res <- dum
		dms
	    }
	    tmp <- unlist(strsplit(META[ML-7], " "))
	    G$n <- dmstodd(tmp[length(tmp)])
	    tmp <- unlist(strsplit(META[ML-6], " "))
	    G$s <- dmstodd(tmp[length(tmp)])
	    tmp <- unlist(strsplit(META[ML-5], " "))
	    G$w <- dmstodd(tmp[length(tmp)])
	    tmp <- unlist(strsplit(META[ML-4], " "))
	    G$e <- dmstodd(tmp[length(tmp)])
	    if (G$e < G$w) G$e <- G$e + 360
	    tmp <- unlist(strsplit(META[ML-3], " "))
	    G$nsres <- dmstodd(tmp[length(tmp)])
	    tmp <- unlist(strsplit(META[ML-2], " "))
	    G$ewres <- dmstodd(tmp[length(tmp)])
	}
	G$Nrow <- as.numeric(unlist(strsplit(META[ML-1], ":"))[2])
	G$Ncol <- as.numeric(unlist(strsplit(META[ML], ":"))[2])
	G$Ncells <- G$Nrow * G$Ncol
	G$xlim <- c(G$w, G$e)
	G$ylim <- c(G$s, G$n)
	G$xseq <- seq(from=G$w + (G$ewres/2), to=G$e - (G$ewres/2), by=G$ewres)
	G$yseq <- seq(from=G$s + (G$nsres/2), to=G$n - (G$nsres/2), by=G$nsres)
	G$ryseq <- rev(G$yseq)
	class(G) <- "grassmeta"
    }
    invisible(G)
}

get.GRASSChkGISRC <- function() {
      get("ChkGISRC", env = .GRASS.meta)
}

GRASS.connect <- function(manual=FALSE, gisdbase=NULL, loc=NULL, mapset=NULL) {
      assign("ChkGISRC", FALSE, env = .GRASS.meta)
      if (manual) {
	      if (any(c(is.null(gisdbase), is.null(loc), is.null(mapset))))
                    stop("missing specification values")
	      z <- .Call("R_G__set_init", as.integer(1), PACKAGE="GRASS")
	      if (z != 1) stop("Init not set")
	      if (!file.exists(gisdbase)) stop(paste(gisdbase, "not found"))
              set.GISDBASE(gisdbase)
              set.LOCATION(loc)
              set.MAPSET(mapset)
              if (class(try(gmeta())) != "grassmeta")
		    stop("cannot retrieve window metadata from specificed location")
              assign("ChkGISRC", TRUE, env = .GRASS.meta)
              assign("gisrc", "faked", env = .GRASS.meta)
              cat("GRASS environment variables inserted manually\n")
            return()
      }
      if (Sys.getenv("GISRC") == "") return()
      gisrc <- .Call("R_G_get_gisrc_file", PACKAGE="GRASS")
      if (class(gisrc) == "gisrc") {  
              assign("ChkGISRC", TRUE, env = .GRASS.meta)
              assign("gisrc", gisrc, env = .GRASS.meta)
              cat("GRASS environment variables in:", gisrc, "\n")
      } else {
              cat("No GRASS environment found\n")
              assign("ChkGISRC", FALSE, env = .GRASS.meta)
              cat(paste("If GRASS.connect() fails in this way",
              "and you are running under CygWin,\nplease set the",
              "CygWin root file system prefix using:",
              "set.cygwinstring()", "\nand re-run GRASS.connect()\n"))
      }
}

make.maas.location <- function() {
	if (get.GRASSChkGISRC())
		stop ("Please run examples and checks outside GRASS")
	data(utm.maas)
	assign("maas.loc", FALSE, env = .GRASS.meta)
	rootdir <- tempdir()
	z <- .Call("R_G__set_init", as.integer(1), PACKAGE="GRASS")
	if (z != 1) stop("Init not set")
	set.GISDBASE(rootdir)
	set.LOCATION("maas")
	set.MAPSET("RGRASS")
	z <- .Call("R_G_make_maas", maas.metadata, PACKAGE="GRASS")
	if (z != 0) stop("error creating location")
	if (class(try(gmeta())) == "grassmeta")
		assign("maas.loc", TRUE, env = .GRASS.meta)
	z <- get("maas.loc", env = .GRASS.meta)
	z
}

set.cygwinstring <- function(cygwin) {
      if (!is.character(cygwin)) stop("Character string required")
      newcygwin <- .Call("R_G_set_cygwinstring", as.character(cygwin)[1],
	PACKAGE="GRASS")
      newcygwin
}

get.cygwinstring <- function() {
      cygwin <- .Call("R_G_get_cygwinstring", PACKAGE="GRASS")
      cygwin
}


set.LOCATION <- function(loc) {
      if (!is.character(loc)) stop("Character string required")
      newloc <- .Call("R_G_set_locstring", as.character(loc)[1], 
	PACKAGE="GRASS")
      newloc
}

get.LOCATION <- function() {
      loc <- .Call("R_G_get_location", PACKAGE="GRASS")
      loc
}


set.GISDBASE <- function(gisdbase) {
      if (!is.character(gisdbase)) stop("Character string required")
      newgisdbase <- .Call("R_G_set_gisdbasestring", as.character(gisdbase)[1],
	PACKAGE="GRASS")
      newgisdbase
}

get.GISDBASE <- function() {
      gisdbase <- .Call("R_G_get_gisdbase", PACKAGE="GRASS")
      gisdbase
}

set.MAPSET <- function(mapset) {
      if (!is.character(mapset)) stop("Character string required")
      newmapset <- .Call("R_G_set_mapset", as.character(mapset)[1], 
	PACKAGE="GRASS")
      newmapset
}

get.MAPSET <- function() {
      mapset <- .Call("R_G_get_mapset", PACKAGE="GRASS")
      mapset
}

# Copyright 1999-2000 by Roger S. Bivand
#
interp.new.G <- function(G, x, y, z, extrap = FALSE,
		duplicate = "error", dupfun = NULL, reverse=NULL) {
    require(akima)
    temp <- interp.new(x, y, z, xo=G$xseq, yo=G$yseq, 
	extrap = extrap, duplicate = duplicate, dupfun = dupfun)
    if(is.null(reverse)) reverse <- reverse(G)
    return(as.vector(temp$z)[reverse])
}
# GRASS adaptation Copyright 1999-2000 by Roger S. Bivand

#
# kde2d and bandwidth.nrd copyright 1994-9 W.N.Venables & B.D.Ripley
#
kde2d.G <- function (G, x, y, h, reverse=NULL, Z=NULL) 
{
    bandwidth.nrd <- function (x) 
    {
        r <- quantile(x, c(0.25, 0.75))
        h <- (r[2] - r[1])/1.34
        4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
    }
    nx <- length(x)
    if (length(y) != nx) 
        stop("Data vectors must be the same length")
    if (missing(h)) 
        h <- c(bandwidth.nrd(x), bandwidth.nrd(y))
    h <- h/4
    ax <- outer(G$xseq, x, "-")/h[1]
    ay <- outer(G$yseq, y, "-")/h[2]
    z <- matrix(dnorm(ax), G$Ncol, nx) %*%
            t(matrix(dnorm(ay), G$Nrow, nx))/(nx * h[1] * h[2])
    if (!is.null(Z)) {
        if (length(Z) != nx) 
            stop("Data vectors must be the same length")
        z1 <- matrix(dnorm(ax), G$Ncol, nx) %*% diag(Z) %*%
            t(matrix(dnorm(ay), G$Nrow, nx))/(nx * h[1] * h[2])
        z <- z1 / z
    }
    if (is.null(reverse)) reverse <- reverse(G)
    return(as.vector(z)[reverse])
}


# Copyright 2001 by Roger S. Bivand
#

krige.G <- function(point.obj, at, var.mod.obj, G, mask=NULL) 
{
    require(sgeostat)
    require(spatial)
    if (!inherits(point.obj, "point")) 
        stop("point.obj must be of class, \"point\".\n")
    if (!inherits(var.mod.obj, "variogram.model")) 
        stop("var.mod.obj must be of class, \"variogram.model\".\n")
    if (class(G) != "grassmeta") 
        stop("G not a grass object")
    at.val <- point.obj[[match(at, names(point.obj))]]
    nugget <- var.mod.obj$parameters[1]
    sill <- var.mod.obj$parameters[2]
    range <- var.mod.obj$parameters[3]
    names(nugget) <- NULL
    names(range) <- NULL
    names(sill) <- NULL
    se <- sqrt(nugget + sill)
    alpha <- nugget / (nugget + sill)
    d <- range
    model <- attr(var.mod.obj, "type")
    kr.model <- NULL
    if (model == "exponential")
        kr.model <- surf.gls(0, expcov, x=point.obj$x, y=point.obj$y,
            z=at.val, d=d, alpha=alpha, se=se)
    else if (model == "gaussian")
        kr.model <- surf.gls(0, gaucov, x=point.obj$x, y=point.obj$y,
            z=at.val, d=d, alpha=alpha, se=se)
    else if (model == "spherical")
        kr.model <- surf.gls(0, sphercov, x=point.obj$x, y=point.obj$y,
            z=at.val, d=d, alpha=alpha, se=se, D=2)
    else stop("unsupported variogram model")
    if (!is.null(mask)) {
        if (length(mask) != G$Ncells)
            stop ("mask length does not equal grid size")
        s <- cbind(east(G)[!is.na(mask)], 
	    north(G)[!is.na(mask)])
    }
    else s <- cbind(east(G), north(G))
    zhat <- prmat2.G(kr.model, s)
    sehat <- semat2.G(kr.model, s)
    if (!is.null(mask)) {
        res1 <- as.numeric(rep(NA,length(mask)))
        res1[!is.na(mask)] <- zhat
        res2 <- as.numeric(rep(NA,length(mask)))
        res2[!is.na(mask)] <- sehat
        return(list(kr=kr.model, zhat=res1, sehat=res2))
    }
    return(list(kr=kr.model, zhat=zhat, sehat=sehat))
}
prmat2.G <- function (obj, s) 
{
    predval <- function(obj, xp, yp) {
        npt <- length(xp)
        .C("VR_krpred", z = double(npt), as.double(xp), as.double(yp), 
            as.double(obj$x), as.double(obj$y), as.integer(npt), 
            as.integer(length(obj$x)), as.double(obj$yy), PACKAGE="spatial")$z
    }
    require(spatial)
    if (!inherits(obj, "trgls")) 
        stop("object not from kriging")
    .C("VR_frset", as.double(obj$rx[1]), as.double(obj$rx[2]), 
        as.double(obj$ry[1]), as.double(obj$ry[2]), PACKAGE="spatial")
    alph <- obj$alph
    if (length(alph) <= 1) {
        mm <- 1.5 * sqrt((obj$rx[2] - obj$rx[1])^2 + (obj$ry[2] - 
            obj$ry[1])^2)
        alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
    }
    .C("VR_alset", as.double(alph), as.integer(length(alph)), PACKAGE="spatial")
    z <-  predict.trls(obj, s[,1], s[,2]) + predval(obj, s[,1], s[,2])
    invisible(z)
}
semat2.G <- function (obj, s, se) 
{
    seval <- function(obj, xp, yp) {
        npt <- length(xp)
        np <- obj$np
        npar <- ((np + 1) * (np + 2))/2
        .C("VR_prvar", z = double(npt), as.double(xp), as.double(yp), 
            as.integer(npt), as.double(obj$x), as.double(obj$y), 
            as.double(obj$l), as.double(obj$r), as.integer(length(obj$x)), 
            as.integer(np), as.integer(npar), as.double(obj$l1f), 
	    PACKAGE="spatial")$z
    }
    require(spatial)
    if (!inherits(obj, "trgls")) 
        stop("object not from kriging")
    .C("VR_frset", as.double(obj$rx[1]), as.double(obj$rx[2]), 
        as.double(obj$ry[1]), as.double(obj$ry[2]), PACKAGE="spatial")
    alph <- obj$alph
    if (length(alph) <= 1) {
        mm <- 1.5 * sqrt((obj$rx[2] - obj$rx[1])^2 + (obj$ry[2] - 
            obj$ry[1])^2)
        alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
    }
    .C("VR_alset", as.double(alph), as.integer(length(alph)), PACKAGE="spatial")
    np <- obj$np
    npar <- ((np + 1) * (np + 2))/2
    if (missing(se)) 
        se <- sqrt(sum(obj$W^2)/(length(obj$x) - npar))
    z <- se * sqrt(seval(obj, s[,1], s[,2]))
    invisible(z)
}
# Copyright 2003 by Roger S. Bivand
#
list.grass <- function(type="cell") {
	types <- c("cell", "site_lists", "dig_plus")
	if (!(type %in% types)) stop(paste(type, "unknown value"))
	mpsts <- get.mapsets()
	res <- vector(mode="list", length=length(mpsts))
	names(res) <- mpsts
	for (i in 1:length(mpsts))
		if (file.exists(paste(get.GISDBASE(), get.LOCATION(), 
			mpsts[i], type, sep="/"))) 
			res[[i]] <- list.files(paste(get.GISDBASE(), 
				get.LOCATION(), mpsts[i], type, sep="/"))
	class(res) <- "glist"
	attr(res, "type") <- type
	res
}

print.glist <- function(x, ...) {
	mpsts <- names(x)
	if (any(!sapply(x, is.null))) {
		cat("----------------------------------------------\n")
		for (i in 1:length(mpsts))
			if (!is.null(x[[i]])) {
				cat(attr(x, "type"), 
					" files avalable in mapset ",
					mpsts[i], ":\n", sep="")
				print(x[[i]], quote=FALSE)
				cat("\n")
			}
		cat("----------------------------------------------\n")
	}
	invisible(x)
}

# Copyright 2003 by Roger S. Bivand

get.mapsets <- function() {
	res <- .Call("R_G_get_mapsets", PACKAGE="GRASS")
	res
}

refresh.mapsets <- function() {
	res <- .Call("R_G_refresh_mapsets", PACKAGE="GRASS")
	res
}

# Copyright 2000 by Roger S. Bivand
#
# north() is an access function to return the northings coordinates
# of raster cell centres from a grassmeta object
#
north <- function(object) UseMethod("north")

north.default <- function(object) stop("no default method for north")

north.grassmeta <- function(object)
{
    if (class(object) != "grassmeta") stop("No GRASS metadata object")
    if(is.loaded("northG", PACKAGE="GRASS")) {
	north <- .Call("northG", object, PACKAGE="GRASS")
    } else {
        north <- as.numeric(c(matrix(object$ryseq, length(object$xseq),
	    length(object$ryseq), byrow=TRUE)))
    }
    invisible(north)
}

# Copyright 2000 by Roger S. Bivand
#
# obsno() is an access function to return the numbers in sequence
# of raster cells. The function is used when NA cells need to be dropped
# for analysis but reinstated later for display or transfer back to GRASS.
#
obsno <- function(G)
{
    if (class(G) != "grassmeta") stop("No GRASS metadata object")
    if(is.loaded("obsnoG", PACKAGE="GRASS")) {
	obsno <- .Call("obsnoG", G, PACKAGE="GRASS")
    } else {
        obsno <- as.integer(1:G$Ncells)
    }
    invisible(obsno)
}

# Copyright 1999-2000 by Roger S. Bivand
#
#
# plot.grassmeta provides a simple interface between grass data
# objects and the image() function; category layers may be plotted
# by taking unclass() of the layer, and setting zlim to non-default values.
# If layer is not set, a blank base map is plotted, for instance for use
# with points().
#
plot.grassmeta <- function(x, layer=NULL, xlab="", ylab="",
    reverse=NULL, add=FALSE, breaks=NULL, ...) {
    G <- x
    if (class(G) != "grassmeta") stop("Data not a grass object")
    if (!add) {
        plot(G$xlim, G$ylim, xlim=G$xlim, ylim=G$ylim, asp=1, xlab = xlab,
            ylab = ylab, type = "n")
    }
    if (!is.null(layer)) {
	if (length(layer) != G$Ncells)
	    stop("GRASS object metadata do not match layer length")
	if (!is.null(breaks)) {
	    if (any(is.na(breaks) | is.nan(breaks))) 
		stop ("NAs found in breaks")
	    if (is.unsorted(breaks)) 
	        stop("`breaks' must be sorted non-decreasingly")
	    if (breaks[1] == -Inf) breaks[1] <- -(.Machine$double.xmax)
	    if (breaks[length(breaks)] == Inf) 
		breaks[length(breaks)] <- .Machine$double.xmax
	}
        if (is.null(reverse)) reverse <- reverse(G)
	if (is.null(breaks)) {
	    image(x=G$xseq, y=G$yseq, z=t(matrix(layer[reverse],
            	nrow=G$Nrow, ncol=G$Ncol, byrow=TRUE)), add=TRUE, 
		...)
	} else {
	    image(x=G$xseq, y=G$yseq, z=t(matrix(layer[reverse],
            	nrow=G$Nrow, ncol=G$Ncol, byrow=TRUE)), add=TRUE, 
		breaks=breaks, ...)
	}
    }
}

legtext <- function(break.levels)
{
    x <- break.levels
    n <- length(x)
    cx <- as.character(x)
    legend <- character(length=(n-1))
    for (i in 1:length(legend)) legend[i] <- paste(x[i], "-", x[i+1], sep="")
    legend
}

leglabs <- function(x1, under="under", over="over", between="-") {
	x <- x1
	res <- character(length(x)-1)
	res[1] <- paste(under, x[2])
	for (i in 2:(length(x)-2)) res[i] <- paste(x[i], between, x[i+1])
	res[length(x)-1] <- paste(over, x[length(x)-1])
	res
}

#findInterval2 <- function (x2, vec, rightmost.closed = FALSE,
# all.inside = TRUE) {
#    x <- x2
#    nx <- length(x)
#    if (any(is.na(vec) | is.nan(vec))) stop ("NAs found in vec")
#    if (is.unsorted(vec)) 
#        stop("`vec' must be sorted non-decreasingly")
#    if (vec[1] == -Inf) vec[1] <- -(.Machine$double.xmax)
#    if (vec[length(vec)] == Inf) 
#	vec[length(vec)] <- .Machine$double.xmax
#    .C("find_interv_vec", xt = as.double(vec), n = length(vec), 
#        x = as.double(x), nx = nx, as.logical(rightmost.closed), 
#        as.logical(all.inside), index = integer(nx), DUP = FALSE,
#	PACKAGE = "base")$index
#}

# Copyright 1999-2003 by Roger S. Bivand
#
#
# rast.get moves one or more GRASS 5.0 raster files to a list, returning
# the filled object. Setting catlabels to TRUE imports category labels 
# instead of codes, and requires more memory.
#
rast.get <- function(G, rlist, catlabels=NULL, drop.unused.levels=FALSE, 
	make.ordered=TRUE, debug=FALSE, interp=FALSE) 
{
    if (class(G) != "grassmeta") stop("No GRASS metadata object")
    if (! is.character(rlist))
	stop("character vector of GRASS data base file names required")
    if (! is.null(catlabels)) {
	if (! is.logical(catlabels))
	    stop("catlabels should be logical vector")
	if (length(catlabels) != length(rlist))
	    stop("catlabels should be same length as rlist")
    } else catlabels <- rep(FALSE, length(rlist))
    
    if(is.loaded("rastget", PACKAGE="GRASS") && (interp == FALSE)) {
	data <- .Call("rastget", G=G, layers=rlist, flayers=catlabels,
		PACKAGE="GRASS")
    } else {

	G.list <- unlist(list.grass(type="cell"))
	res <- rlist %in% G.list
	if (! all(res)) {
		warning("The following GRASS data base files were not found:")
		print(rlist[res == FALSE])
		stop("transfer terminated")
	}
	data <- vector(mode="list", length=length(rlist))
	ndata <- character(length=length(rlist))
	for (i in 1:length(rlist)) {
	    FILE <- tempfile("GRtoR")
	    if (catlabels[i]) {
		rstats <- "r.stats -1ql fs=\":\" input="
		rstats <- paste(rstats, rlist[i], ",", sep="")
		rstats <- paste(rstats, " output=", FILE, sep="")
		system(rstats)
		x <- scan(FILE, what=(list(double(0), character(0))), 
		    sep=":", na.strings="*", quiet=TRUE)
		if (length(x[[1]]) != G$Ncells)
		    stop("Number of rows imported does not match metadata")
		x[[2]][is.na(x[[1]])] <- NA
		ndata[i] <- paste(rlist[i], ".f", sep="")
		data[[i]] <- factor(x[[2]], levels=unique(x[[2]]), ordered=TRUE)
		rm(x)
	    } else {
		rstats <- "r.stats -1q fs=\":\" input="
		rstats <- paste(rstats, rlist[i], ",", sep="")
		rstats <- paste(rstats, " output=", FILE, sep="")
		system(rstats)
		x <- scan(FILE, na.strings="*", quiet=TRUE)
		if (length(x) != G$Ncells)
		    stop("Number of rows imported does not match metadata")
		ndata[i] <- rlist[i]
		data[[i]] <- x
		rm(x)
	    }
	    if (!debug) unlink(FILE)
	}
	names(data) <- ndata
    }
    ndata <- names(data)
    names(data) <- make.names(names=ndata, unique=TRUE)
    if (!is.null(catlabels)) {
	for (i in 1:length(data)) {
	    if (catlabels[i]) {
		if (drop.unused.levels) data[[i]] <- data[[i]][, drop=TRUE]
		if (!make.ordered) class(data[[i]]) <- "factor"
	    }
	}
    }
    invisible(data)
}
# Copyright 1999-2001 by Roger S. Bivand
#
#
# rast.put moves a single numeric vector to GRASS, using the metadata
# retrieved by gmeta() from the GRASS data base.
#
rast.put <- function(G, lname="", layer, title="", cat=FALSE, DCELL=FALSE,
	breaks=NULL, col=NULL, nullcol=NULL, defcol=NULL, debug=FALSE,
	interp=FALSE, check=TRUE) 
    {
    if (class(G) != "grassmeta") stop("Data not a grass object")
    if (length(lname) != 1)
	stop("Single new GRASS data base file name required")
    if (length(layer) != G$Ncells)
	stop("GRASS object metadata do not match layer length")
    if (!(is.numeric(layer) || is.factor(layer)))
	stop("layer is neither numeric nor factor")
    if(is.loaded("rastput", PACKAGE="GRASS") && (interp == FALSE)) {
	if(!is.logical(check)) stop("check must be logical")
	if(is.null(nullcol)) nullcol <- "honeydew"
	if(is.null(defcol)) defcol <- "pale turquoise"
	nullcolor <- as.integer(col2rgb(nullcol[1]))
	defcolor <- as.integer(col2rgb(defcol[1]))
	if(cat || is.factor(layer)) {
	    if(!is.null(breaks)) warning("breaks ignored for factor layers")
	    if (is.null(col)) {
		col <- rev(grey(1:length(levels(layer))/length(levels(layer))))
	    }
	    else if(length(levels(layer)) != length(col))
		stop("number of colors must equal number of factor levels")
	    color <- as.integer(col2rgb(col))
	    layer.range <- range(na.omit(unclass(layer)))
	        x <- .Call("rastput", G=G, layer=as.integer(unclass(layer)),
		isfactor=TRUE, DCELL=FALSE, check=as.logical(check), 
		levels=levels(layer), output=lname, title=title, breaks=NULL,
		color=as.integer(color), nullcolor=as.integer(nullcolor),
		as.integer(defcolor), range=as.integer(layer.range),
		PACKAGE="GRASS")
	} else {
# check col/breaks Jos Agustn Garca Garca 12/12-03
	    if(is.null(breaks)) {
		breaks <- pretty(as.double(na.omit(layer)), n=20, min.n=10)
	    } else {
		if (!is.numeric(breaks))
		    stop("non-numeric breaks not accepted")
	    }
	    if (is.null(col)) {
		col <- rev(grey(1:(length(breaks)-1)/(length(breaks)-1)))
	    } else if(length(breaks) != (length(col)-1))
		stop("number of colors must equal one less than the number of breaks")
	    layer.levels <- character((length(breaks)-1))
	    for (i in 1:(length(breaks)-1)) {
		layer.levels[i] <- paste("(", signif(breaks[i]), ",",
		    signif(breaks[i+1]), "]", sep="")
	    }
	    col <- as.integer(col2rgb(col))
	    layer.range <- range(breaks)
	    x <- .Call("rastput", G=G, layer=as.double(layer), isfactor=FALSE, 
		DCELL=FALSE, check=as.logical(check), levels=layer.levels, 
		output=lname, title=title, breaks=as.double(breaks), 
		color=as.integer(col), nullcolor=as.integer(nullcolor), 
		defcolor=as.integer(defcolor),	range=as.double(layer.range),
		PACKAGE="GRASS")
	}
    } else {

	G.list <- list.grass(type="cell")
	res <- lname %in% G.list[[match(get.MAPSET(), get.mapsets())]]
	if (any(res))
		stop(paste(lname, ": GRASS raster file already exists ", 
			"in mapset: ", get.MAPSET(), sep=""))
	if (length(layer) != G$Ncells)
		stop("GRASS object metadata do not match layer length")
	FILE <- tempfile("RtoGR")
	outstr <- paste("north:   ", G$n, "\nsouth:   ", G$s, "\neast:    ", 
		G$e, "\nwest:    ", G$w, "\nrows:    ", G$Nrow, 
		"\ncols:    ", G$Ncol, "\n", sep="")
	cat(outstr, file=FILE)
	if (cat) write(t(matrix(as.integer(unclass(layer)), nrow=G$Nrow,
		ncol=G$Ncol, byrow=TRUE)), file=FILE, append=TRUE,
		ncolumns=G$Ncol)
	else write(t(matrix(as.double(layer), nrow=G$Nrow, ncol=G$Ncol, 
		byrow=TRUE)), file=FILE, append=TRUE, ncolumns=G$Ncol)
	if (cat) system(paste("r.in.ascii -i input=", FILE, " nv=NA output=", 
		lname, " title=\"", title, "\"", sep=""))
	else {
	    if (DCELL)
	        system(paste("r.in.ascii -d input=", FILE, " nv=NA output=", 
	  	    lname, " title=\"", title, "\"", sep=""))
	    else
	        system(paste("r.in.ascii -f input=", FILE, " nv=NA output=", 
	  	    lname, " title=\"", title, "\"", sep=""))
	}
	if (!debug) unlink(FILE)
    }
}
# Copyright 2000 by Roger S. Bivand
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
# reverse() is an access function to return the reversed order of
# raster cells for use with image()
#
reverse <- function(G)
{
    if (class(G) != "grassmeta") stop("No GRASS metadata object")
    if(is.loaded("reverseG", PACKAGE="GRASS")) {
	reversed <- .Call("reverseG", G, PACKAGE="GRASS")
    } else {
	eastG <- east(G)
	northG <- north(G)
        reversed <- as.integer(order(northG, eastG))
    }
    invisible(reversed)
}




# Copyright 1999-2003 by Roger S. Bivand
#
# sites.get moves one GRASS 5.0 sites file to a data frame, returning
# the filled object. 
#
sites.get <- function(G, slist = "", all.sites=FALSE, collapse.labels=TRUE, 
	debug=FALSE, interp=FALSE) {
	if (class(G) != "grassmeta") stop("No GRASS metadata object")
	if (! is.character(slist))
		stop("character GRASS data base file name required")
	if (!is.loaded("sitesget", PACKAGE="GRASS")) interp <- TRUE
	if (interp) {

		G.list <- unlist(list.grass(type="site_lists"))
		res <- slist %in% G.list
		if (! all(res)) stop(paste("transfer terminated: ",
			slist[res == FALSE], "GRASS data base file not found"))
		if (all.sites) allsites <- " -a"
		else allsites <- ""
		FILE <- tempfile("GRtoR")
		system(paste("s.out.ascii -d", allsites, " sites=", slist,
			" > ", FILE, sep=""))
		data <- read.table(FILE, na.strings="*")
# CHANGE 000329 RSB Only expect eastings and northings, not id as before
		nc2 <- ncol(data) - 2
		nlist <- character(0)
		for (i in 1:nc2) nlist <- c(nlist, paste("var", i, sep=""))
		names(data) <- c("east", "north", nlist)
		if (!debug) unlink(FILE)
	} else {
		res <- .Call("sitesget", G, slist, all.sites, PACKAGE="GRASS")
		nncols <- attr(res, "nncols")
		if (length(res) == 1) {
			data <- as.data.frame(res)
			names(data) <- c("east", "north")
		} else {
			data <- as.data.frame(res[[1]])
			nlist <- c("east", "north")
			if (length(res[[1]]) > 2) {
				for (i in 3:length(res[[1]])) nlist <- c(nlist, 
					paste("dim", i, sep=""))
			}
			data <- cbind(data, as.data.frame(res[[2]]))
			nlist <- c(nlist, "id")
			if (nncols[3] > 0) {
				data <- cbind(data, as.data.frame(res[[3]]))
				for (i in 1:length(res[[3]])) nlist <- c(nlist, 
					paste("num", i, sep=""))
			}
			if (nncols[4] > 0) {
# changes 2003/3/28 to accept old-style labels
				if (nncols[2] < 0 &&
				    collapse.labels) {
				    oldatt <- res[[4]][[1]]
				    for (i in 2:length(res[[4]]))
					oldatt <- paste(oldatt, res[[4]][[i]])
				    data <- cbind(data, as.data.frame(oldatt))
				    nlist <- c(nlist, paste("str1", sep=""))
				} else {
				    data <- cbind(data, as.data.frame(res[[4]]))
				    for (i in 1:length(res[[4]])) nlist <- 
					c(nlist, paste("str", i, sep=""))
				}
			}
			if (!is.null(attr(res, "labels"))) {
				labs <- unlist(strsplit(attr(res, 
					"labels"), " "))
				if (length(labs) == ncol(data)) {
					colnames(data) <- make.names(labs)
				} else {
					colnames(data) <- make.names(nlist)
				}
			} else {
				colnames(data) <- make.names(nlist)
			}
		}
	}
	attr(data, "nncols") <- nncols
	invisible(data)
}
# Copyright 1999-2003 by Roger S. Bivand
#
# sites.put moves a single variable site to GRASS, using the metadata
# prepared when layers were got from the GRASS data base.
#
sites.put <- function(G, lname="", east, north, var, 
	debug=FALSE) {
	warning("This function may be withdrawn, consider using sites.put2()")
	if (class(G) != "grassmeta") stop("Data not a grass object")
	if (length(lname) != 1)
		stop("Single new GRASS data base file name required")

	G.list <- list.grass(type="site_lists")
	res <- lname %in% G.list[[match(get.MAPSET(), get.mapsets())]]
	if (any(res))
		stop(paste(lname, ": GRASS sites file already exists ", 
		"in mapset: ", get.MAPSET(), sep=""))
	if (length(east) != length(north))
		stop("Different numbers of eastings and northings")
	if (length(east) != length(var))
		stop("Different numbers of coordinates and observations")
	inregion <- (east >= G$w & east <= G$e) & (north >= G$s & north <= G$n)
	if(all(!inregion)) stop("None of the site locations are inside the current GRASS region")
	if(any(!inregion)) warning("Some site locations are outside the current GRASS region")
	FILE <- tempfile("RtoGR")
	if (is.numeric(var))
		a <- data.frame(x=east, y=north, z=as.character(paste("%", var, sep="")))
	else
		a <- data.frame(x=east, y=north, z=as.character(paste("@", var, sep="")))
	write.table(a, row.names=FALSE, col.names=FALSE, quote=FALSE, file=FILE)
	system(paste("s.in.ascii input=", FILE, " sites=", lname, sep=""))
	if (!debug) unlink(FILE)
}

sites.put2 <- function(G, data, id=NULL, dims, lname="", all.sites=FALSE,
	check=TRUE) {
	if (class(G) != "grassmeta") stop("Data not a grass object")
	if (class(data) != "data.frame") stop("Data frame required")
	nas <- unlist(lapply(data, function(x) any(is.na(x))))
	if (any(nas)) stop("NAs cannot be moved to GRASS sites files")
	if (!is.loaded("sitesput", PACKAGE="GRASS")) stop(paste("sitesput",
		"compiled function not loaded, use old sitesput() function"))	

	if(!is.logical(check)) stop("check must be logical")
	n <- nrow(data)
	if (is.null(id)) {
		ids <- as.integer(1:n)
		idname <- "id"
		xid <- NULL
	} else {
		if (length(id) > 1) stop ("single id required")
		if (is.character(id)) {
			xid <- match(id, colnames(data))
			if (is.na(xid))
				stop ("id not found")
			ids <- data[,xid]
			idname <- id
		} else {
			if (is.integer(id) && id > 0 && id <= ncol(data)) {
				ids <- data[,id]
				idname <- colnames(data)[id]
				xid <- id
			} else {
				stop ("id not valid number")
			}
		}
	}
	if(!is.numeric(ids)) stop ("id not numeric")
	if(is.integer(ids)) cattype <- 0
	else cattype <- 2

	ndims <- length(dims)
	if (ndims < 2) stop("less than two dimensions")
	if (is.character(dims)) {
		xdims <- match(dims, names(data))
		if (any(is.na(xdims)))
			stop ("dims not found")
		dims.mat <- as.matrix(data[, dims])
		dimsnames <- dims
	} else if (is.integer(dims)) {
		if (any(dims < 1) || any (dims > ncol(data))) 
			stop ("dims not valid")
		dims.mat <- as.matrix(data[, dims])
		dimsnames <- names(data)[dims]
		xdims <- dims
	} else stop ("dims not valid")
	if (!is.numeric(dims.mat)) stop ("dims not numeric")
	if (is.integer(dims.mat)) dims.mat <- as.numeric(dims.mat)

	if (length(lname) != 1)
		stop("Single new GRASS data base file name required")
	if (!is.logical(all.sites)) stop("all.sites not logical")
	if (!is.logical(check)) stop("check not logical")

	if (is.null(xid)) drops <- xdims
	else drops <- c(xid, xdims)
	dnames <- colnames(data)
	data.attr <- as.data.frame(data[, -drops])
	da.ncol <- ncol(data.attr)
	dblnames <- NULL
	dbl.mat <- NULL
	strnames <- NULL
	str.mat <- NULL
	if (da.ncol == 0) {
		xnumeric <- NULL
		xother <- NULL
		ndbls <- 0
		nstrs <- 0
		dblnames <- NULL
		strnames <- NULL
		warning("No attributes transferred, only dimensions")
	} else {
		dnames <- dnames[-drops]

		xnumeric <- which(unlist(lapply(data.attr, is.numeric)))
		xother <- which(unlist(lapply(data.attr,
			function(x) !is.numeric(x))))
		ndbls <- length(xnumeric)
		
		if (length(xnumeric) > 0) {
			if (da.ncol == 1) {
				dbl.mat <- matrix(data.attr[,1], ncol=1, nrow=n)
			} else {
				dbl.mat <- as.matrix(data.attr[, xnumeric])
			}
			dblnames <- dnames[xnumeric]
			if (is.integer(dbl.mat)) dbl.mat <- as.numeric(dbl.mat)
		}

		nstrs <- length(xother)
	
		if (nstrs > 0) {
			if (da.ncol == 1) {
				str.mat <- matrix(as.character(data.attr[,1]), 
					ncol=1, nrow=n)
			} else {
				str.mat <- as.matrix(data.attr[, xother])
			}
			strnames <- dnames[xother]
		}
	}

	labs <- paste(c(dimsnames, idname, dblnames, strnames), collapse=" ")
	n.args <- as.integer(c(cattype, n, ndims, ndbls, nstrs))
	call <- deparse(match.call(), width=500)

	res <- list(G=G, lname=lname, n.args=n.args, all.sites=all.sites, 
		labs=labs, ids=ids, dims.mat=dims.mat, dbl.mat=dbl.mat, 
		str.mat=str.mat, call=call, check=check)
	xx <- .Call("sitesput", res, PACKAGE="GRASS")

	invisible(xx)
}
# Copyright 1999-2000 by Roger S. Bivand
#
# summary.grassmeta displays the metadata prepared when map layers are moved to R.
#
# CHANGED 000329 RSB added projection output
summary.grassmeta <- function(object, ...) {
	G <- object
	if (class(G) != "grassmeta") stop("Data not a grass object")
	cat("Data from GRASS 5.0 LOCATION ", G$LOCATION, " with ", G$Ncol,
	" columns and ", G$Nrow, " rows;\n", G$proj, "\nThe west-east range is: ",
	G$w, ", ", G$e, ",\nand the south-north: ",
	G$s, ", ", G$n,
	";\nWest-east cell sizes are ", G$ewres,
	" units,\nand south-north ", G$nsres,
	" units.\n", sep="")
}
# GRASS adaptation Copyright 1999-2001 by Roger S. Bivand
#

trmat.G <- function (G, obj, east=NULL, north=NULL) 
{
    require(spatial)
    if (!inherits(obj, "trls")) 
        stop("object not a fitted trend surface")
    if (class(G) != "grassmeta") 
        stop("Data not a grass object")
    if (is.null(east)) east <- east(G)
    if (is.null(north)) north <- north(G)
    z <- predict(obj, east, north) 
    invisible(z)
}


# Copyright 2000-3 by Roger S. Bivand. 
#

.GRASS.meta <- new.env(FALSE, globalenv())
assign("ChkGISRC", FALSE, env = .GRASS.meta)
assign("maas.loc", FALSE, env = .GRASS.meta)

.First.lib <- function(lib, pkg) {
	library.dynam("GRASS", pkg, lib)
	.Call("R_G_init", "RGRASS_INTERFACE", PACKAGE="GRASS")
	GRASS.connect()
}

