R/rast_link.R

In rgrass7: Deprecated Interface Between GRASS Geographical Information System and R

# Interpreted GRASS 7, 8 raster interface functions
# Copyright (c) 2022 Roger S. Bivand
#

read_RAST <- function(vname, cat=NULL, NODATA=NULL, 
    ignore.stderr=get.ignore.stderrOption(), return_format="SGDF", 
    close_OK=return_format=="SGDF", flags=NULL) {
        .Deprecated(new="read_RAST", package="rgrass", old="read_RAST",
            msg="Package rgrass7 transitioning to package rgrass for GRASS 7+.")
    if (!is.null(cat))
        if(length(vname) != length(cat)) 
	    stop("vname and cat not same length")
    if (get.suppressEchoCmdInFuncOption()) {
        inEchoCmd <- set.echoCmdOption(FALSE)
    }
    if (close_OK) {
        openedConns <- as.integer(row.names(showConnections()))
    }
    stopifnot(is.logical(ignore.stderr))

    if (!is.null(NODATA)) {
        if (any(!is.finite(NODATA)) || any(!is.numeric(NODATA)))
            stop("invalid NODATA value")
	if (NODATA != round(NODATA)) 
	    warning("NODATA rounded to integer")
	NODATA <- round(NODATA)
        if (length(NODATA) != length(vname)) 
            NODATA <- rep(NODATA[1], length.out=length(vname))
     }

    if (return_format == "SGDF") {
        if (!(requireNamespace("sp", quietly=TRUE)))
            stop("sp required for SGDF output")
        resa <- .read_rast_non_plugin_ng(vname=vname, cat=cat, NODATA=NODATA, 
            ignore.stderr=ignore.stderr)
        if (close_OK) { #closeAllConnections()
            openConns_now <- as.integer(row.names(showConnections()))
            toBeClosed <- openConns_now[!(openConns_now %in% openedConns)]
            for (bye in toBeClosed) close(bye)
        }
    } else {
        if (!(requireNamespace("terra", quietly=TRUE))) 
            stop("terra required for SpatRaster output")
        drv <- "RRASTER"
        fxt <- ".grd"
        ro <- FALSE
        o <- execGRASS("r.out.gdal", flags="l", intern=TRUE)
        oo <- grep("RRASTER", o)
        if (length(oo) == 0L) ro <- TRUE
        if (!ro) {
            RR <- o[oo]
            RRs <- strsplit(RR, " ")[[1]]
            if (length(grep("\\(rw", RRs)) == 0L) ro <- TRUE
        }
        if (ro) {
            drv <- "GTiff"
            fxt <- ".tif"
        }
        reslist <- vector(mode="list", length=length(vname))
        names(reslist) <- vname
        tmplist <- vector(mode="list", length=length(vname))
        names(tmplist) <- vname
        for (i in seq(along=vname)) {
# 130422 at rgdal 0.8-8 GDAL.close(DS)
# 061107 Dylan Beaudette NODATA
# 071009 Markus Neteler's idea to use range
            if (is.null(NODATA)) {
	        tx <- execGRASS("r.info", flags="r", map=vname[i], intern=TRUE, 
	            ignore.stderr=ignore.stderr)
	        tx <- gsub("=", ":", tx)
	        con <- textConnection(tx)
                res <- read.dcf(con)
	        close(con)
	        lres <- as.list(res)
	        names(lres) <- colnames(res)
	        lres <- lapply(lres, function(x)
                ifelse(x == "NULL", as.numeric(NA), as.numeric(x)))
	        if (!is.numeric(lres$min) || 
	           !is.finite(as.double(lres$min))) 
                     NODATAi <- as.integer(999)
	        else {
	            lres$min <- floor(as.double(lres$min))
	            NODATAi <- floor(lres$min) - 1
	        }
            } else NODATAi <- NODATA[i]
            tmplist[[i]] <- tempfile(fileext=fxt)
            if (is.null(flags)) flags <- c("overwrite", "c", "m")
            if (!is.null(cat) && cat[i]) flags <- c(flags, "t")
            execGRASS("r.out.gdal", input=vname[i], output=tmplist[[i]],
                format=drv, nodata=NODATAi, flags=flags,
                ignore.stderr=ignore.stderr)
            reslist[[i]] <- getMethod("rast", "character")(tmplist[[i]])
        }         
        resa <- getMethod("rast", "list")(reslist)
    }
    if (get.suppressEchoCmdInFuncOption()) tull <- set.echoCmdOption(inEchoCmd)
    
    resa
}


.read_rast_non_plugin_ng <- function(vname, cat, NODATA, ignore.stderr) {
    
    pid <- as.integer(round(runif(1, 1, 1000)))
        
    gLP <- getLocationProj()
    if (gLP == "XY location (unprojected)")
        p4 <- sp::CRS(as.character(NA))
    else
        p4 <- sp::CRS(gLP)

    reslist <- vector(mode="list", length=length(vname))
    names(reslist) <- vname

    for (i in seq(along=vname)) {

        glist <- execGRASS("r.info", flags="g", map=vname[i],
            intern=TRUE, ignore.stderr=ignore.stderr)
        whCELL <- glist[grep("datatype", glist)]
	to_int <- length(which(unlist(strsplit(
	    whCELL, "=")) == "CELL")) > 0
        Dcell <- length(which(unlist(strsplit(
	    whCELL, "=")) == "DCELL")) > 0

	gtmpfl1 <- dirname(execGRASS("g.tempfile",
	    pid=pid, intern=TRUE, ignore.stderr=ignore.stderr))
	rtmpfl1 <- ifelse(.Platform$OS.type == "windows" &&
            (Sys.getenv("OSTYPE") == "cygwin"), 
	    system(paste("cygpath -w", gtmpfl1, sep=" "), 
	    intern=TRUE), gtmpfl1)
	gtmpfl11 <- paste(gtmpfl1, vname[i], sep=.Platform$file.sep)
	rtmpfl11 <- paste(rtmpfl1, vname[i], sep=.Platform$file.sep)

# 130422 at rgdal 0.8-8 GDAL.close(DS)
# 061107 Dylan Beaudette NODATA
# 071009 Markus Neteler's idea to use range
    if (is.null(NODATA)) {
	tx <- execGRASS("r.info", flags="r", map=vname[i], intern=TRUE, 
	    ignore.stderr=ignore.stderr)
	tx <- gsub("=", ":", tx)
	con <- textConnection(tx)
        res <- read.dcf(con)
	close(con)
	lres <- as.list(res)
	names(lres) <- colnames(res)
	lres <- lapply(lres, function(x)
        ifelse(x == "NULL", as.numeric(NA), as.numeric(x)))
	if (!is.numeric(lres$min) || 
	    !is.finite(as.double(lres$min))) 
            NODATAi <- as.integer(999)
	else {
	    lres$min <- floor(as.double(lres$min))
	    NODATAi <- floor(lres$min) - 1
	}
    }

        rOutBinFlags <- "b"
# 120118 Rainer Krug
        if (to_int) rOutBinFlags <- c(rOutBinFlags, "i")
        else rOutBinFlags <- c(rOutBinFlags, "f")
        rOutBinBytes <- 4L
        if (Dcell) rOutBinBytes <- 8L
        tryCatch({
	    execGRASS("r.out.bin", flags=rOutBinFlags, 
                input=vname[i], output=gtmpfl11, bytes=rOutBinBytes,
                null=as.integer(NODATAi), ignore.stderr=ignore.stderr)
                        
                gdal_info <- bin_gdal_info_ng(rtmpfl11, to_int)
                        
                what <- ifelse(to_int, "integer", "double")
                n <- gdal_info[1] * gdal_info[2]
                size <- gdal_info[10]/8
                        
                reslist[[i]] <- readBinGridData(rtmpfl11, what=what,
                    n=n, size=size, endian=attr(gdal_info, "endian"),
                    nodata=attr(gdal_info, "nodata"))
            }, finally = {
                unlink(paste(rtmpfl1, list.files(rtmpfl1,
                    pattern=vname[i]), sep=.Platform$file.sep))
            }
        )
    }

    co <- unname(c((gdal_info[4] + (gdal_info[6]/2)),
        (gdal_info[5] + (gdal_info[7]/2))))
    grid <- sp::GridTopology(co, unname(c(gdal_info[6], gdal_info[7])),
        unname(c(gdal_info[2], gdal_info[1])))

    if (length(unique(sapply(reslist, length))) != 1)
        stop ("bands differ in length")

    df <- as.data.frame(reslist)

    resa <- sp::SpatialGridDataFrame(grid=grid, data=df, proj4string=p4)

    if (!is.null(cat)) {
	for (i in seq(along=cat)) {
	    if (cat[i] && is.integer(resa@data[[i]])) {

	        rSTATS <- execGRASS("r.stats", flags=c("l", "quiet"),
	            input=vname[i], intern=TRUE,
                    ignore.stderr=ignore.stderr)

		cats <- strsplit(rSTATS, " ")
		catnos <- sapply(cats, function(x) x[1])
		catlabs <- sapply(cats, 
		    function(x) paste(x[-1], collapse=" "))
		if (any(!is.na(match(catnos, "*")))) {
		    isNA <- which(catnos == "*")
		    catnos <- catnos[-isNA]
		    catlabs <- catlabs[-isNA]
		}
                if (length(catlabs) > length(unique(catlabs))) {
                    catlabs <- paste(catlabs, catnos, sep="_")
                    warning("non-unique category labels; category number appended")
                }
# https://files.nc.gov/ncdeq/Energy+Mineral+and+Land+Resources/Geological+Survey/1985_state_geologic_map_500000_scale.pdf (catnos vector polygon IDs)
		resa@data[[i]] <- factor(resa@data[[i]], 
		    levels=catnos, labels=catlabs)
            }
	} 
    }
    return(resa)
}




bin_gdal_info_ng <- function(fname, to_int) {
	if (!file.exists(fname)) stop(paste("no such file:", fname))
	if (!file.exists(paste(fname, "hdr", sep="."))) 
		stop(paste("no such file:", paste(fname, "hdr", sep=".")))
	if (!file.exists(paste(fname, "wld", sep="."))) 
		stop(paste("no such file:", paste(fname, "wld", sep=".")))
	con <- file(paste(fname, "hdr", sep="."), "r")
	l8 <- readLines(con, n=8)
	close(con)
	l8 <- read.dcf(con <- textConnection(gsub(" ", ":", l8)))
	close(con)
	lres <- as.list(l8)
	names(lres) <- colnames(l8)
	lres$nrows <- as.integer(lres$nrows)
	lres$ncols <- as.integer(lres$ncols)
	lres$nbands <- as.integer(lres$nbands)
	lres$nbits <- as.integer(lres$nbits)
	lres$skipbytes <- as.integer(lres$skipbytes)
	lres$nodata <- ifelse(to_int, as.integer(lres$nodata), 
		as.numeric(lres$nodata))
	lres$byteorder <- as.character(lres$byteorder)
	endian <- .Platform$endian
	if ((endian == "little" && lres$byteorder == "M") ||
		(endian == "big" && lres$byteorder == "I")) endian <- "swap"
	con <- file(paste(fname, "wld", sep="."), "r")
	l6 <- readLines(con, n=6)
	close(con)
	lres$ewres <- abs(as.numeric(l6[1]))
	lres$nsres <- abs(as.numeric(l6[4]))
	lres$n_cc <- as.numeric(l6[6])
	lres$w_cc <- as.numeric(l6[5])
	lres$s_cc <- lres$n_cc - lres$nsres * (lres$nrows-1)

    outres <- numeric(10)
    outres[1] <- lres$nrows
    outres[2] <- lres$ncols
    outres[3] <- lres$nbands
    outres[4] <- lres$w_cc - (lres$ewres/2)
    outres[5] <- lres$s_cc - (lres$nsres/2)
    outres[6] <- lres$ewres
    outres[7] <- lres$nsres
    outres[10] <- lres$nbits
    attr(outres, "endian") <- endian
    attr(outres, "nodata") <- lres$nodata
#	grid = GridTopology(c(lres$w_cc, lres$s_cc), 
#		c(lres$ewres, lres$nsres), c(lres$ncols,lres$nrows))
    outres
}

readBinGridData <- function(fname, what, n, size, endian, nodata) {
    map <- readBin(fname, what=what, n=n, size=size, signed=TRUE,
	endian=endian)
    is.na(map) <- map == nodata
    map
}

repair_cats <- function(x, layer, vname) {
    if (!(requireNamespace("terra", quietly=TRUE))) 
        stop("terra required for to repair terra object cats")
    if (!inherits(x, "SpatRaster")) stop("x not a SpatRaster object")
    xcats <- getMethod("cats", "SpatRaster")(x)
    nms <- getMethod("names", "SpatRaster")(x)
    if (length(xcats) > 1) {
        nms <- getMethod("names", "SpatRaster")(x)
        if (!(layer %in% nms)) stop("layer not found")
    }
    xcat <- xcats[[match(layer, nms)]]
    if (is.null(xcat)) warning("no cats in layer")
    zcats <- cats_internals(vname)
    res <- vector(mode="list", length=1)
    names(res) <- layer
    res[[1]] <- getMethod("classify", "SpatRaster")(x[[layer]], zcats$rcp_out)
    getMethod("set.cats", "SpatRaster")(res[[1]], layer=1, value=zcats$levs)
    getMethod("rast", "list")(res)
}

cats_internals <- function(vname) {
    cats <- execGRASS("r.category", map=vname, separator=":", intern=TRUE)
    cats1 <- strsplit(cats, ":")
    cats2 <- data.frame(ID=sapply(cats1, "[", 1L), 
        category=sapply(cats1, "[", 2))
    cats2$Icat <- as.integer(factor(cats2$category,
        levels=unique(cats2$category)))
    cats3 <- split(cats2, cats2$Icat)
    rcp_out <- do.call("rbind", lapply(cats3, function(x) {
        cbind(as.integer(x[[1]]), as.integer(x[[3]]))
    }))
    levs <- do.call("rbind", lapply(cats3, function(x) {
        ux <- unique(x[, 2:3])
        data.frame(cat=as.integer(ux[[2]]), lev=ux[[1]])
    }))
    list(rcp_out=rcp_out, levs=levs)
}

read_cat_colors <- function(vname) {
    zcats <- cats_internals(vname)
    cols <- execGRASS("r.colors.out", map=vname, intern=TRUE)
    cols12 <- strsplit(cols, " ")
    cat_col1 <- match(as.character(zcats$rcp_out[,1]), sapply(cols12, "[", 1))
    if (length(cat_col1) != (length(cols12) - 2L)) 
        warning("possible category mismatch")
    coltb <- sapply(cols12[zcats$rcp_out[order(zcats$rcp_out[,1]), 2]], "[", 2L)
    df <- as.data.frame(t(sapply(strsplit(coltb, ":"), as.integer)))
    pal <- rgb(df$V1, df$V2, df$V3, maxColorValue=255)
    list(coltab=df, pal=pal)
}


write_RAST <- function(x, vname, zcol = 1, NODATA=NULL, flags=NULL, 
    ignore.stderr = get.ignore.stderrOption(), overwrite=FALSE, verbose=TRUE) {
        .Deprecated(new="write_RAST", package="rgrass", old="write_RAST",
            msg="Package rgrass7 transitioning to package rgrass for GRASS 7+.")
    if (get.suppressEchoCmdInFuncOption()) {
        inEchoCmd <- set.echoCmdOption(FALSE)
    }
    stopifnot(is.logical(ignore.stderr))
    stopifnot(is.character(vname))
    if (!is.null(flags)) stopifnot(is.character(flags))

    if (inherits(x, "SpatialGridDataFrame")) {
        if (!(requireNamespace("sp", quietly=TRUE))) 
            stop("sp required for SGDF input")
        stopifnot(length(vname) == 1L)
        pid <- as.integer(round(runif(1, 1, 1000)))
        gtmpfl1 <- dirname(execGRASS("g.tempfile", pid=pid,
            intern=TRUE, ignore.stderr=ignore.stderr))
                
        rtmpfl1 <- ifelse(.Platform$OS.type == "windows" &&
            (Sys.getenv("OSTYPE") == "cygwin"), 
            system(paste("cygpath -w", gtmpfl1, sep=" "), intern=TRUE), 
            gtmpfl1)
                
        fid <- paste("X", pid, sep="")
        gtmpfl11 <- paste(gtmpfl1, fid, sep=.Platform$file.sep)
        rtmpfl11 <- paste(rtmpfl1, fid, sep=.Platform$file.sep)
        if (!is.numeric(x@data[[zcol]])) 
            stop("only numeric columns may be exported")
        if (overwrite && !("overwrite" %in% flags))
            flags <- c(flags, "overwrite")
        res <- writeBinGrid_ng(x, rtmpfl11, attr = zcol, na.value = NODATA)
                        
        flags <- c(res$flag, flags)
                        
        execGRASS("r.in.bin", flags=flags, input=gtmpfl11, output=vname,
            bytes=as.integer(res$bytes), north=as.numeric(res$north), 
            south=as.numeric(res$south), east=as.numeric(res$east), 
            west=as.numeric(res$west), rows=as.integer(res$rows), 
            cols=as.integer(res$cols), anull=as.numeric(res$anull), 
            ignore.stderr=ignore.stderr)
        if (verbose) cat("SpatialGridDataFrame read into GRASS using r.in.bin\n")

    } else if (inherits(x, "SpatRaster")) {
        if (!(requireNamespace("terra", quietly=TRUE))) 
            stop("terra required for terra output")
# Suggestion https://github.com/rsbivand/rgrass/pull/45#discussion_r816113064 Floris Vanderhaeghe
        srcs <- getMethod("sources", "SpatRaster")(x)
        mems <- getMethod("inMemory", "SpatRaster")(x)
        if (length(srcs) == 1L && !mems[1]) tf <- srcs[1]
        else tf <- ""
        if (!file.exists(tf)) {
            drv <- "RRASTER"
            fxt <- ".grd"
            gdalver <- gsub("[A-Za-z]", "", strsplit(terra::gdal(), "-")[[1]][1])
            if (gdalver < "2.3.0") {
                drv <- "GTiff"
                fxt <- ".tif"
            }
            tf <- tempfile(fileext=fxt)
            res <- getMethod("writeRaster", c("SpatRaster", "character"))(x,
                filename=tf, overwrite=TRUE, filetype=drv)
            tmpfl <- TRUE
        } else {
            res <- x
            tmpfl <- FALSE
        }
        execGRASS("r.in.gdal", flags=flags, input=tf, output=vname)
#        if (tmpfl) unlink(tf)
        if (verbose) cat("SpatRaster read into GRASS using r.in.gdal from",
            ifelse(tmpfl, "memory", "file"), "\n")
        if (getMethod("nlyr", "SpatRaster")(x) == 1L) {
            xcats <- getMethod("cats", "SpatRaster")(x)[[1]]
            if (!is.null(xcats)) {
                tfc <- tempfile()
                write.table(xcats, tfc, sep=":", row.names=FALSE, 
                    col.names=FALSE, quote=FALSE)
                execGRASS("r.category", map=vname, rules=tfc, separator=":")
            }
        }
    } else {
        stop("object neither SpatialGridDataFrame nor SpatRaster")
    }
    if (get.suppressEchoCmdInFuncOption()) {
        tull <- set.echoCmdOption(inEchoCmd)
    }
    invisible(res)
}

writeBinGrid_ng <- function(x, fname, attr = 1, na.value = NULL) { 
	if (!sp::gridded(x))
		stop("can only write SpatialGridDataFrame objects to binary grid")
	x = as(x, "SpatialGridDataFrame")
	gp = sp::gridparameters(x)
	if (length(gp$cells.dim) != 2)
		stop("binary grid only supports 2D grids")
	z = x@data[[attr]]
	if (is.factor(z)) z <- as.numeric(z)
	if (!is.numeric(z)) stop("only numeric values may be exported")
	if (is.null(na.value)) {
		na.value <- floor(min(z, na.rm=TRUE)) - 1
	} else {
		if (!is.finite(na.value) || !is.numeric(na.value))
			stop("invalid NODATA value")
		if (na.value != round(na.value)) 
			warning("NODATA rounded to integer")
		na.value <- round(na.value)
	}
	res <- list()
	res$anull <- formatC(na.value, format="d")
	z[is.na(z)] = as.integer(na.value)
	if (storage.mode(z) == "integer") {
		sz <- 4
	} else if (storage.mode(z) == "double") {
		sz <- 8
		res$flag <- "d"
	} else stop("unknown storage mode")
	res$bytes <- formatC(sz, format="d")
	f = file(fname, open = "wb")
	writeBin(z, con=f, size=sz)
	close(f)
	grd <- slot(x, "grid")
	f = file(paste(fname, "hdr", sep="."), open = "wt")
	writeLines(paste("nrows", grd@cells.dim[2]), f)
	res$rows <- formatC(grd@cells.dim[2], format="d")
	writeLines(paste("ncols", grd@cells.dim[1]), f)
	res$cols <- formatC(grd@cells.dim[1], format="d")
	writeLines(paste("nbands 1"), f)
	writeLines(paste("nbits", 8*sz), f)
	writeLines(paste("byteorder", ifelse(.Platform$endian == "little", 
		"I", "M")), f)
	writeLines(paste("layout bil"), f)
	writeLines(paste("skipbytes 0"), f)
	writeLines(paste("nodata", na.value), f)
	close(f)
	f = file(paste(fname, "wld", sep="."), open = "wt")
	writeLines(formatC(grd@cellsize[1], format="f"), f)
	writeLines("0.0", f)
	writeLines("0.0", f)
	writeLines(formatC(-grd@cellsize[2], format="f"), f)
	writeLines(formatC(grd@cellcentre.offset[1], format="f"), f)
	writeLines(formatC((grd@cellcentre.offset[2] +
		grd@cellsize[2]*(grd@cells.dim[2]-1)), format="f"), f)
	close(f)
	res$north <- formatC((grd@cellcentre.offset[2] +
		grd@cellsize[2]*(grd@cells.dim[2]-1)
		+ 0.5*grd@cellsize[2]), format="f")
	res$south <- formatC(grd@cellcentre.offset[2] - 0.5*grd@cellsize[2], 
		format="f")
	res$east <- formatC((grd@cellcentre.offset[1] + 
		grd@cellsize[1]*(grd@cells.dim[1]-1) + 0.5*grd@cellsize[1]), 
		format="f")
	res$west <- formatC(grd@cellcentre.offset[1] - 0.5*grd@cellsize[1], 
		format="f")
	invisible(res)
}