Nothing
# Interpreted GRASS 6+ interface functions
# Copyright (c) 2005-2013 Roger S. Bivand
#
readRAST <- readRAST6 <- function(vname, cat=NULL, ignore.stderr = NULL,
NODATA=NULL, plugin=NULL, mapset=NULL, useGDAL=NULL, close_OK=TRUE,
drivername="GTiff", driverFileExt=NULL, return_SGDF=TRUE) {
if (!is.null(cat))
if(length(vname) != length(cat))
stop("vname and cat not same length")
if (get.suppressEchoCmdInFuncOption()) {
inEchoCmd <- get.echoCmdOption()
tull <- set.echoCmdOption(FALSE)
}
if (is.null(plugin))
plugin <- get.pluginOption()
stopifnot(is.logical(plugin)|| is.null(plugin))
if (!is.null(plugin) && plugin && length(vname) > 1) plugin <- FALSE
if (is.null(ignore.stderr))
ignore.stderr <- get.ignore.stderrOption()
stopifnot(is.logical(ignore.stderr))
if (is.null(useGDAL))
useGDAL <- get.useGDALOption()
stopifnot(is.logical(useGDAL))
if (useGDAL) {
if (requireNamespace("rgdal", quietly = TRUE)) {
gdalD <- rgdal::gdalDrivers()$name
} else {
stop("rgdal not available")
}
}
if (!useGDAL && is.null(plugin)) plugin <- FALSE
if (close_OK) {
openedConns <- as.integer(row.names(showConnections()))
}
if (is.null(plugin)) plugin <- "GRASS" %in% gdalD
if (length(vname) > 1) plugin <- FALSE
if (plugin) {
resa <- read_plugin(vname, mapset=NULL, ignore.stderr=ignore.stderr)
} else {
pid <- as.integer(round(runif(1, 1, 1000)))
p4 <- CRS(getLocationProj())
Gver <- execGRASS("g.version", intern=TRUE,
ignore.stderr=ignore.stderr)
G63 <- !(Gver < "GRASS 6.3") #Gver > "GRASS 6.2"
# 090311 fix for -c flag
Cflag <- "c" %in% parseGRASS("r.out.gdal")$fnames
reslist <- vector(mode="list", length=length(vname))
names(reslist) <- vname
for (i in seq(along=vname)) {
whCELL <- execGRASS("r.info", flags="t",
map=vname[i], intern=TRUE,
ignore.stderr=ignore.stderr)
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)
if (!is.null(driverFileExt)) {
gtmpfl11 <- paste(gtmpfl1, driverFileExt, sep=".")
rtmpfl11 <- paste(rtmpfl1, driverFileExt, sep=".")
}
if (!is.null(NODATA)) {
if (!is.finite(NODATA) || !is.numeric(NODATA))
stop("invalid NODATA value")
if (NODATA != round(NODATA))
warning("NODATA rounded to integer")
NODATA <- round(NODATA)
}
if (useGDAL && G63) {
if (requireNamespace("rgdal", quietly = TRUE)) {
gdalDGRASS <- execGRASS("r.out.gdal", flags="l",
intern=TRUE, ignore.stderr=ignore.stderr)
if (!(drivername %in% gdalD))
stop(paste("Requested driver", drivername,
"not available in rgdal"))
if (length(grep(drivername, gdalDGRASS)) == 0)
stop(paste("Requested driver", drivername,
"not available in GRASS"))
type <- ifelse (to_int, "Int32",
ifelse(Dcell, "Float64", "Float32"))
flags <- c("quiet")
if (Cflag) flags <- c("c", flags)
paras <- list(input=vname[i], output=gtmpfl11, type=type,
format=drivername)
if (!is.null(NODATA)) paras <- c(paras, nodata=NODATA)
execGRASS("r.out.gdal", flags=flags,
parameters=paras, ignore.stderr=ignore.stderr)
# res <- readGDAL(rtmpfl11, p4s=getLocationProj(),
# silent=ignore.stderr)
# names(res) <- vname[i]
if (i == 1) gdal_info <- rgdal::GDALinfo(rtmpfl11,
silent=ignore.stderr)
DS <- rgdal::GDAL.open(rtmpfl11, read.only=FALSE)
reslist[[i]] <- as.vector(rgdal::getRasterData(DS, band=1))
rgdal::deleteDataset(DS)
} else {
stop("rgdal not available")
}
# 130422 at rgdal 0.8-8 GDAL.close(DS)
} else {
# 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)))
NODATA <- as.integer(999)
else {
lres$min <- floor(as.double(lres$min))
NODATA <- floor(lres$min) - 1
}
}
rOutBinFlags <- "b"
tryCatch(
{
# 120118 Rainer Krug
if (Gver < "GRASS 6.4.2") {
if (to_int) rOutBinFlags <- c(rOutBinFlags, "i")
execGRASS("r.out.bin", flags=rOutBinFlags,
input=vname[i], output=gtmpfl11,
null=as.integer(NODATA),
ignore.stderr=ignore.stderr)
} else {
if (to_int) rOutBinFlags <- c(rOutBinFlags, "i")
else rOutBinFlags <- c(rOutBinFlags, "f")
rOutBinBytes <- 4L
if (Dcell) rOutBinBytes <- 8L
execGRASS("r.out.bin", flags=rOutBinFlags,
input=vname[i], output=gtmpfl11, bytes=rOutBinBytes,
null=as.integer(NODATA),
ignore.stderr=ignore.stderr)
}
gdal_info <- bin_gdal_info(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))
}
)
}
# if (i == 1) resa <- res
# else {
# grida <- getGridTopology(resa)
# grid <- getGridTopology(res)
# if (!isTRUE(all.equal(grida, grid)))
# stop("topology is not equal")
# onames <- c(names(resa@data), names(res@data))
# ncols <- dim(resa@data)[2]
# lst <- vector(mode="list", length=ncols+1)
# names(lst) <- onames
# for (i in 1:ncols) lst[[i]] <- resa@data[[i]]
# lst[[ncols+1]] <- res@data[[1]]
# df <- data.frame(lst)
# resa <- SpatialGridDataFrame(grid=grida,
# data=df, proj4string=p4)
# }
}
co <- unname(c((gdal_info[4] + (gdal_info[6]/2)),
(gdal_info[5] + (gdal_info[7]/2))))
grid <- GridTopology(co, unname(c(gdal_info[6], gdal_info[7])),
unname(c(gdal_info[2], gdal_info[1])))
if (close_OK) { #closeAllConnections()
openConns_now <- as.integer(row.names(showConnections()))
toBeClosed <- openConns_now[!(openConns_now %in% openedConns)]
for (bye in toBeClosed) close(bye)
}
if (!return_SGDF) {
res <- list(grid=grid, dataList=reslist, proj4string=p4)
class(res) <- "gridList"
if (get.suppressEchoCmdInFuncOption()) {
tull <- set.echoCmdOption(inEchoCmd)
}
return(res)
}
if (length(unique(sapply(reslist, length))) != 1)
stop ("bands differ in length")
df <- as.data.frame(reslist)
resa <- 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]
}
resa@data[[i]] <- factor(resa@data[[i]],
levels=catnos, labels=catlabs)
}
}
}
}
if (get.suppressEchoCmdInFuncOption()) {
tull <- set.echoCmdOption(inEchoCmd)
}
resa
}
bin_gdal_info <- 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
}
#readBinGrid <- function(fname, colname=basename(fname),
# proj4string=CRS(as.character(NA)), integer) {
# if (missing(integer)) stop("integer TRUE/FALSE required")
# 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(integer, 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)
#
# what <- ifelse(integer, "integer", "double")
# n <- lres$nrows * lres$ncols
# size <- lres$nbits/8
# map <- readBin(fname, what=what, n=n, size=size, signed=TRUE,
# endian=endian)
# is.na(map) <- map == lres$nodata
# grid = GridTopology(c(lres$w_cc, lres$s_cc),
# c(lres$ewres, lres$nsres), c(lres$ncols,lres$nrows))
# df <- list(var1=map)
# names(df) <- colname
# df1 <- data.frame(df)
#
# res <- SpatialGridDataFrame(grid, data = df1, proj4string=proj4string)
# res
#}
read_plugin <- function(vname, mapset=NULL, ignore.stderr=NULL) {
if (is.null(ignore.stderr))
ignore.stderr <- get.ignore.stderrOption()
stopifnot(is.logical(ignore.stderr))
if (length(vname) > 1) stop("single raster required for plugin")
gg <- gmeta()
if (is.null(mapset)) {
c_at <- strsplit(vname[1], "@")[[1]]
if (length(c_at) == 1) {
mapset <- .g_findfile(vname[1], type="cell")
} else if (length(c_at) == 2) {
mapset <- c_at[2]
vname[1] <- c_at[1]
} else stop("malformed raster name")
}
fname <- paste(gg$GISDBASE, gg$LOCATION_NAME, mapset,
"cellhd", vname[1], sep="/")
if (requireNamespace("rgdal", quietly = TRUE)) {
fninfo <- rgdal::GDALinfo(fname, silent=ignore.stderr)
chks <- logical(4)
names(chks) <- c("cols", "rows", "origin.northing",
"origin.easting")
chks[1] <- isTRUE(all.equal(abs((gg$w-gg$e)/gg$ewres), fninfo[2],
tol=2e-7, check.attributes=FALSE))
chks[2] <- isTRUE(all.equal(abs((gg$n-gg$s)/gg$nsres), fninfo[1],
tol=2e-7, check.attributes=FALSE))
# changed from gg$n 100129, thanks to Rainer Krug
chks[3] <- isTRUE(all.equal(gg$s, fninfo[5],
check.attributes=FALSE))
chks[4] <- isTRUE(all.equal(gg$w, fninfo[4],
check.attributes=FALSE))
if (any(!chks)) {
cat("raster map/current region mismatch detected in components:\n")
print(chks)
stop("rerun with plugin=FALSE\n")
}
resa <- rgdal::readGDAL(fname, silent=ignore.stderr)
names(resa) <- make.names(vname)
resa
} else {
stop("rgdal not available")
}
}
writeRAST <- writeRAST6 <- function(x, vname, zcol = 1, NODATA=NULL,
ignore.stderr = NULL, useGDAL=NULL, overwrite=FALSE, flags=NULL,
drivername="GTiff") {
if (get.suppressEchoCmdInFuncOption()) {
inEchoCmd <- get.echoCmdOption()
tull <- set.echoCmdOption(FALSE)
}
if (is.null(ignore.stderr))
ignore.stderr <- get.ignore.stderrOption()
stopifnot(is.logical(ignore.stderr))
if (is.null(useGDAL))
useGDAL <- get.useGDALOption()
stopifnot(is.logical(useGDAL))
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")
Gver <- execGRASS("g.version", intern=TRUE,
ignore.stderr=ignore.stderr)
G63 <- !(Gver < "GRASS 6.3") #Gver > "GRASS 6.2"
if (useGDAL && G63) {
if (requireNamespace("rgdal", quietly = TRUE)) {
gdalD <- rgdal::gdalDrivers()$name
gdalDGRASS <- execGRASS("r.out.gdal", flags="l",
intern=TRUE, ignore.stderr=ignore.stderr)
if (!(drivername %in% gdalD))
stop(paste("Requested driver", drivername,
"not available in rgdal"))
if (length(grep(drivername, gdalDGRASS)) == 0)
stop(paste("Requested driver", drivername,
"not available in GRASS"))
if (is.factor(x@data[[zcol]]))
x@data[[zcol]] <- as.numeric(x@data[[zcol]])
if (!is.numeric(x@data[[zcol]]))
stop("only numeric values may be exported")
if (is.null(NODATA)) {
NODATA <- floor(min(x@data[[zcol]], na.rm=TRUE)) - 1
} else {
if (!is.finite(NODATA) || !is.numeric(NODATA))
stop("invalid NODATA value")
if (NODATA != round(NODATA))
warning("NODATA rounded to integer")
NODATA <- round(NODATA)
}
sm <- storage.mode(x[[zcol]])
type <- ifelse(sm == "integer", "Int32", "Float32")
res <- rgdal::writeGDAL(x[zcol], fname=rtmpfl11, type=type,
drivername=drivername, mvFlag = NODATA)
execGRASS("r.in.gdal", flags=flags, input=gtmpfl11,
output=vname, ignore.stderr=ignore.stderr)
DS <- rgdal::GDAL.open(rtmpfl11, read.only=FALSE)
rgdal::deleteDataset(DS)
} else {
stop("rgdal not available")
}
} else {
res <- writeBinGrid(x, rtmpfl11, attr = zcol, na.value = NODATA)
flags <- c(res$flag, flags)
tryCatch(
{
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)
},
finally = {
unlink(paste(rtmpfl1, list.files(rtmpfl1, pattern=fid),
sep=.Platform$file.sep))
}
)
}
if (get.suppressEchoCmdInFuncOption()) {
tull <- set.echoCmdOption(inEchoCmd)
}
invisible(res)
}
writeBinGrid <- function(x, fname, attr = 1, na.value = NULL) {
if (!gridded(x))
stop("can only write SpatialGridDataFrame objects to binary grid")
x = as(x, "SpatialGridDataFrame")
gp = 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.