Nothing
R/merge.R
In raster: Geographic Data Analysis and Modeling
# Author: Robert J. Hijmans
# Date : October 2008
# Version 0.9
# Licence GPL v3
# redesinged for multiple row processing
# and arguments ext and overlap
# October 2011
# version 1
setMethod('merge', signature(x='Extent', y='ANY'),
function(x, y, ...) {
x <- c(x, y, list(...))
x <- sapply(x, extent)
x <- x[sapply(x, function(x) inherits(x, 'Extent'))]
x <- lapply(x, function(e) t(bbox(e)))
x <- do.call(rbind, x)
x <- apply(x, 2, range)
extent(as.vector(x))
}
)
setMethod('merge', signature(x='RasterStackBrick', y='missing'),
function(x, ..., tolerance=0.05, filename="", ext=NULL) {
nl <- nlayers(x)
if (nl < 2) {
return(x)
} else if (nl == 2) {
merge(x[[1]], x[[2]], tolerance=tolerance, filename=filename, overlap=TRUE, ext=ext)
} else {
do.call(merge, c(x=x[[1]], y=x[[2]], .makeRasterList(x[[3:nl]]), tolerance=tolerance, filename=filename, overlap=TRUE, ext=ext))
}
}
)
setMethod('merge', signature(x='Raster', y='Raster'),
function(x, y, ..., tolerance=0.05, filename="", overlap=TRUE, ext=NULL) {
x <- c(x, y, list(...))
isRast <- sapply(x, function(x) inherits(x, 'Raster'))
dotargs <- x[ !isRast ]
x <- x[ isRast ]
compareRaster(x, extent=FALSE, rowcol=FALSE, orig=TRUE, res=TRUE, tolerance=tolerance)
if (is.null(dotargs$datatype)) {
dotargs$datatype <- .commonDataType(sapply(x, dataType))
}
filename <- trim(filename)
dotargs$filename <- filename
nl <- max(unique(sapply(x, nlayers)))
bb <- .unionExtent(x)
if (nl > 1) {
out <- brick(x[[1]], values=FALSE, nl=nl)
} else {
out <- raster(x[[1]])
}
out <- setExtent(out, bb, keepres=TRUE, snap=FALSE)
hasV <- sapply(x, hasValues)
if (!any(hasV)) {
return(out)
}
if (!is.null(ext)) {
ext <- extent(ext)
out1 <- extend(out, union(ext, extent(out)))
out1 <- crop(out1, ext)
test <- try( intersect(extent(out), extent(out1)) )
if (inherits(test, "try-error")) {
stop('"ext" does not overlap with any of the input data')
}
out <- out1
ext <- extent(out)
}
if ( canProcessInMemory(out, 3) ) {
if (!is.null(ext)) {
if (overlap) {
if (nl > 1) {
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells, ]
dat <- extract(x[[i]], ext)
if (!is.matrix(dat)) {
dat <- matrix(dat, ncol=1)
}
na <- ! rowSums(dat)==nl
vv[na, ] <- dat[na, ]
v[cells, ] <- vv
}
}
} else {
v <- rep(NA, ncell(out))
for (i in length(x):1) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
xy <- xyFromCell(out, cells)
d <- extract(x[[i]], xy)
j <- !is.na(d)
v[cells[j]] <- d[j]
}
}
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
} else { # ignore overlap (if any)
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in length(x):1 ) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
v[cells, ] <- extract(x[[i]], ext)
}
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
}
} else {
if (overlap) {
if (nl > 1) {
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in 1:length(x)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells, ]
dat <- getValues(x[[i]])
if (!is.matrix(dat)) {
dat <- matrix(dat, ncol=1)
}
na <- ! rowSums(is.na(dat)) == nl
vv[na, ] <- dat[na, ]
v[cells, ] <- vv
}
} else {
v <- rep(NA, ncell(out))
for (i in 1:length(x)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells]
vv[is.na(vv)] <- getValues(x[[i]])[is.na(vv)]
v[cells] <- vv
}
}
rm(vv)
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
} else { # no overlap (or ignore overlap)
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in length(x):1) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
v[cells, ] <- getValues(x[[i]])
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
}
}
}
if (is.null(ext)) {
rowcol <- matrix(NA, ncol=6, nrow=length(x))
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1) # first row/col on old raster[[i]]
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) ) # last row/col on old raster[[i]]
rowcol[i,1] <- rowFromY(out, xy1[2]) # start row on new raster
rowcol[i,2] <- rowFromY(out, xy2[2]) # end row
rowcol[i,3] <- colFromX(out, xy1[1]) # start col
rowcol[i,4] <- colFromX(out, xy2[1]) # end col
rowcol[i,5] <- i # layer
rowcol[i,6] <- nrow(x[[i]])
}
tr <- blockSize(out)
# tr$row <- sort(unique(c(tr$row, rowcol[,1], rowcol[,2]+1)))
# tr$row <- subset(tr$row, tr$row <= nrow(out))
# tr$nrows <- c(tr$row[-1], nrow(out)+1) - c(tr$row)
# tr$n <- length(tr$row)
pb <- pbCreate(tr$n, dotargs$progress, label='merge')
dotargs$x <- out
out <- do.call(writeStart, dotargs)
if (overlap) {
if (nl == 1) {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i], ncol=ncol(out))
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) { #reverse order so that the first raster covers the second etc.
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
vv[z1:z2, rc[j,3]:rc[j,4]] <- matrix(getValues(x[[ rc[j,5] ]], r1, nr), nrow=nr, byrow=TRUE)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, as.vector(t(v)), tr$row[i])
pbStep(pb, i)
}
} else {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) {
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
vv[cells, ] <- getValues(x[[ rc[j,5] ]], r1, nr)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
} else { # not overlap
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
for (j in nrow(rc):1) {
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
v[cells, ] <- getValues(x[[ rc[j,5] ]], r1, nr)
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
}
} else { # ext is not null
rowcol <- matrix(NA, ncol=10, nrow=length(x))
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1) # first row/col on old raster[[i]]
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) ) # last row/col on old raster[[i]]
xyout1 <- xyFromCell(out, 1)
xyout2 <- xyFromCell(out, ncell(out))
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
j <- rowFromY(out, xy1[2])
rowcol[i,1] <- ifelse(is.na(j), 1, j) # start row on new raster
j <- rowFromY(out, xy2[2])
rowcol[i,2] <- ifelse(is.na(j), nrow(out), j) # end row
j <- colFromX(out, xy1[1])
rowcol[i,3] <- ifelse(is.na(j), 1, j) # start col
j <- colFromX(out, xy2[1])
rowcol[i,4] <- ifelse(is.na(j), ncol(out), j) # end col
rowcol[i,5] <- nrow(x[[i]])
j <- rowFromY(x[[i]], xyout1[2])
rowcol[i,6] <- ifelse(is.na(j), 1, j)
j <- rowFromY(x[[i]], xyout2[2])
rowcol[i,7] <- ifelse(is.na(j), nrow(x[[i]]), j) - rowcol[i,6] + 1
j <- colFromX(x[[i]], xyout1[1])
rowcol[i,8] <- ifelse(is.na(j), 1, j)
j <- colFromX(x[[i]], xyout2[1])
rowcol[i,9] <- ifelse(is.na(j), ncol(x[[i]]), j) - rowcol[i,8] + 1
rowcol[i,10] <- i # layer
}
}
rowcol <- subset(rowcol, !is.na(rowcol[,1]))
tr <- blockSize(out)
# tr$row <- sort(unique(c(tr$row, rowcol[,1], rowcol[,2]+1)))
# tr$row <- subset(tr$row, tr$row <= nrow(out))
# tr$nrows <- c(tr$row[-1], nrow(out)+1) - c(tr$row)
# tr$n <- length(tr$row)
pb <- pbCreate(tr$n, dotargs$progress, label='merge')
dotargs$x <- out
out <- do.call(writeStart, dotargs)
if (overlap) {
if (nl == 1) {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i], ncol=ncol(out))
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) { #reverse order so that the first raster covers the second etc.
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
vv[z1:z2, rc[j,3]:rc[j,4]] <- matrix(getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9]), nrow=nr, byrow=TRUE)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, as.vector(t(v)), tr$row[i])
pbStep(pb, i)
}
} else {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) {
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
vv[cells, ] <- getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9])
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
} else { # no overlap
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
for (j in nrow(rc):1) {
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
v[cells, ] <- getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9])
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
}
pbClose(pb)
writeStop(out)
}
)
Try the raster package in your browser
Any scripts or data that you put into this service are public.
raster documentation built on Jan. 5, 2021, 3:01 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Author: Robert J. Hijmans
# Date : October 2008
# Version 0.9
# Licence GPL v3
# redesinged for multiple row processing
# and arguments ext and overlap
# October 2011
# version 1
setMethod('merge', signature(x='Extent', y='ANY'),
function(x, y, ...) {
x <- c(x, y, list(...))
x <- sapply(x, extent)
x <- x[sapply(x, function(x) inherits(x, 'Extent'))]
x <- lapply(x, function(e) t(bbox(e)))
x <- do.call(rbind, x)
x <- apply(x, 2, range)
extent(as.vector(x))
}
)
setMethod('merge', signature(x='RasterStackBrick', y='missing'),
function(x, ..., tolerance=0.05, filename="", ext=NULL) {
nl <- nlayers(x)
if (nl < 2) {
return(x)
} else if (nl == 2) {
merge(x[[1]], x[[2]], tolerance=tolerance, filename=filename, overlap=TRUE, ext=ext)
} else {
do.call(merge, c(x=x[[1]], y=x[[2]], .makeRasterList(x[[3:nl]]), tolerance=tolerance, filename=filename, overlap=TRUE, ext=ext))
}
}
)
setMethod('merge', signature(x='Raster', y='Raster'),
function(x, y, ..., tolerance=0.05, filename="", overlap=TRUE, ext=NULL) {
x <- c(x, y, list(...))
isRast <- sapply(x, function(x) inherits(x, 'Raster'))
dotargs <- x[ !isRast ]
x <- x[ isRast ]
compareRaster(x, extent=FALSE, rowcol=FALSE, orig=TRUE, res=TRUE, tolerance=tolerance)
if (is.null(dotargs$datatype)) {
dotargs$datatype <- .commonDataType(sapply(x, dataType))
}
filename <- trim(filename)
dotargs$filename <- filename
nl <- max(unique(sapply(x, nlayers)))
bb <- .unionExtent(x)
if (nl > 1) {
out <- brick(x[[1]], values=FALSE, nl=nl)
} else {
out <- raster(x[[1]])
}
out <- setExtent(out, bb, keepres=TRUE, snap=FALSE)
hasV <- sapply(x, hasValues)
if (!any(hasV)) {
return(out)
}
if (!is.null(ext)) {
ext <- extent(ext)
out1 <- extend(out, union(ext, extent(out)))
out1 <- crop(out1, ext)
test <- try( intersect(extent(out), extent(out1)) )
if (inherits(test, "try-error")) {
stop('"ext" does not overlap with any of the input data')
}
out <- out1
ext <- extent(out)
}
if ( canProcessInMemory(out, 3) ) {
if (!is.null(ext)) {
if (overlap) {
if (nl > 1) {
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells, ]
dat <- extract(x[[i]], ext)
if (!is.matrix(dat)) {
dat <- matrix(dat, ncol=1)
}
na <- ! rowSums(dat)==nl
vv[na, ] <- dat[na, ]
v[cells, ] <- vv
}
}
} else {
v <- rep(NA, ncell(out))
for (i in length(x):1) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
xy <- xyFromCell(out, cells)
d <- extract(x[[i]], xy)
j <- !is.na(d)
v[cells[j]] <- d[j]
}
}
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
} else { # ignore overlap (if any)
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in length(x):1 ) {
xy1 <- xyFromCell(x[[i]], 1)
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) )
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
v[cells, ] <- extract(x[[i]], ext)
}
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
}
} else {
if (overlap) {
if (nl > 1) {
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in 1:length(x)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells, ]
dat <- getValues(x[[i]])
if (!is.matrix(dat)) {
dat <- matrix(dat, ncol=1)
}
na <- ! rowSums(is.na(dat)) == nl
vv[na, ] <- dat[na, ]
v[cells, ] <- vv
}
} else {
v <- rep(NA, ncell(out))
for (i in 1:length(x)) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
vv <- v[cells]
vv[is.na(vv)] <- getValues(x[[i]])[is.na(vv)]
v[cells] <- vv
}
}
rm(vv)
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
} else { # no overlap (or ignore overlap)
v <- matrix(NA, nrow=ncell(out), ncol=nl)
for (i in length(x):1) {
cells <- cellsFromExtent( out, extent(x[[i]]) )
v[cells, ] <- getValues(x[[i]])
}
out <- setValues(out, v)
if (filename != '') {
dotargs$x <- out
out <- do.call(writeRaster, dotargs)
}
return(out)
}
}
}
if (is.null(ext)) {
rowcol <- matrix(NA, ncol=6, nrow=length(x))
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1) # first row/col on old raster[[i]]
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) ) # last row/col on old raster[[i]]
rowcol[i,1] <- rowFromY(out, xy1[2]) # start row on new raster
rowcol[i,2] <- rowFromY(out, xy2[2]) # end row
rowcol[i,3] <- colFromX(out, xy1[1]) # start col
rowcol[i,4] <- colFromX(out, xy2[1]) # end col
rowcol[i,5] <- i # layer
rowcol[i,6] <- nrow(x[[i]])
}
tr <- blockSize(out)
# tr$row <- sort(unique(c(tr$row, rowcol[,1], rowcol[,2]+1)))
# tr$row <- subset(tr$row, tr$row <= nrow(out))
# tr$nrows <- c(tr$row[-1], nrow(out)+1) - c(tr$row)
# tr$n <- length(tr$row)
pb <- pbCreate(tr$n, dotargs$progress, label='merge')
dotargs$x <- out
out <- do.call(writeStart, dotargs)
if (overlap) {
if (nl == 1) {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i], ncol=ncol(out))
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) { #reverse order so that the first raster covers the second etc.
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
vv[z1:z2, rc[j,3]:rc[j,4]] <- matrix(getValues(x[[ rc[j,5] ]], r1, nr), nrow=nr, byrow=TRUE)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, as.vector(t(v)), tr$row[i])
pbStep(pb, i)
}
} else {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) {
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
vv[cells, ] <- getValues(x[[ rc[j,5] ]], r1, nr)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
} else { # not overlap
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
for (j in nrow(rc):1) {
r1 <- tr$row[i]-rc[j,1]+1
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,6], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
v[cells, ] <- getValues(x[[ rc[j,5] ]], r1, nr)
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
}
} else { # ext is not null
rowcol <- matrix(NA, ncol=10, nrow=length(x))
for (i in 1:length(x)) {
xy1 <- xyFromCell(x[[i]], 1) # first row/col on old raster[[i]]
xy2 <- xyFromCell(x[[i]], ncell(x[[i]]) ) # last row/col on old raster[[i]]
xyout1 <- xyFromCell(out, 1)
xyout2 <- xyFromCell(out, ncell(out))
if (xy1[2] > ymin(out) & xy2[2] < ymax(out) & xy1[1] < xmax(out) & xy2[1] > xmin(out)) {
j <- rowFromY(out, xy1[2])
rowcol[i,1] <- ifelse(is.na(j), 1, j) # start row on new raster
j <- rowFromY(out, xy2[2])
rowcol[i,2] <- ifelse(is.na(j), nrow(out), j) # end row
j <- colFromX(out, xy1[1])
rowcol[i,3] <- ifelse(is.na(j), 1, j) # start col
j <- colFromX(out, xy2[1])
rowcol[i,4] <- ifelse(is.na(j), ncol(out), j) # end col
rowcol[i,5] <- nrow(x[[i]])
j <- rowFromY(x[[i]], xyout1[2])
rowcol[i,6] <- ifelse(is.na(j), 1, j)
j <- rowFromY(x[[i]], xyout2[2])
rowcol[i,7] <- ifelse(is.na(j), nrow(x[[i]]), j) - rowcol[i,6] + 1
j <- colFromX(x[[i]], xyout1[1])
rowcol[i,8] <- ifelse(is.na(j), 1, j)
j <- colFromX(x[[i]], xyout2[1])
rowcol[i,9] <- ifelse(is.na(j), ncol(x[[i]]), j) - rowcol[i,8] + 1
rowcol[i,10] <- i # layer
}
}
rowcol <- subset(rowcol, !is.na(rowcol[,1]))
tr <- blockSize(out)
# tr$row <- sort(unique(c(tr$row, rowcol[,1], rowcol[,2]+1)))
# tr$row <- subset(tr$row, tr$row <= nrow(out))
# tr$nrows <- c(tr$row[-1], nrow(out)+1) - c(tr$row)
# tr$n <- length(tr$row)
pb <- pbCreate(tr$n, dotargs$progress, label='merge')
dotargs$x <- out
out <- do.call(writeStart, dotargs)
if (overlap) {
if (nl == 1) {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i], ncol=ncol(out))
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) { #reverse order so that the first raster covers the second etc.
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
vv[z1:z2, rc[j,3]:rc[j,4]] <- matrix(getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9]), nrow=nr, byrow=TRUE)
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, as.vector(t(v)), tr$row[i])
pbStep(pb, i)
}
} else {
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
vv <- v
for (j in nrow(rc):1) {
vv[] <- NA
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
vv[cells, ] <- getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9])
v[!is.na(vv)] <- vv[!is.na(vv)]
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
} else { # no overlap
for (i in 1:tr$n) {
v <- matrix(NA, nrow=tr$nrow[i]*ncol(out), ncol=nl)
rc <- subset(rowcol, (tr$row[i]+tr$nrow[i]-1) > rowcol[,1] & tr$row[i] < rowcol[,2])
if (nrow(rc) > 0) {
for (j in nrow(rc):1) {
r1 <- tr$row[i]-rc[j,1]+rc[j,6]
r2 <- r1 + tr$nrow[i]-1
z1 <- abs(min(1,r1)-1)+1
r1 <- max(1, r1)
r2 <- min(rc[j,5], r2)
nr <- r2 - r1 + 1
z2 <- z1 + nr - 1
cells <- cellFromRowColCombine(out, z1:z2, rc[j,3]:rc[j,4])
v[cells, ] <- getValuesBlock(x[[ rc[j,10] ]], r1, nr, rc[j,8], rc[j,9])
}
}
out <- writeValues(out, v, tr$row[i])
pbStep(pb, i)
}
}
}
pbClose(pb)
writeStop(out)
}
)
Try the raster package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.