R/indexReplace.R

In miraisolutions/raster: Geographic Data Analysis and Modeling

# Author: Robert J. Hijmans
# Date :  January 2009
# Version 1.0
# Licence GPL v3


setReplaceMethod("[", c("RasterLayer", "RasterLayer", "missing"),
	function(x, i, j, value) {

		i <- crop(i, x)
		
		if (inherits(value, 'RasterLayer')) {
			value <- getValues(value)
		}
		
		if (! hasValues(i) ) {
			i <- cellsFromExtent(x, i)
			
		} else if (compareRaster(x, i, stopiffalse=FALSE, showwarning=FALSE)) {
			i <- as.logical( getValues(i) )
		
		} else {
			j <- as.logical( getValues(i) )
			i <- cellsFromExtent(x, i)[j]
			x[i] <- value
			return(x)
		}		
	
		.replace(x, i, value=value, recycle=1) 
	}
)



setReplaceMethod("[", c("RasterLayer","missing","missing"),
	function(x, i, j, value) {
	
		if (length(value) == ncell(x)) {
			x <- try( setValues(x, value))
		} else if (length(value) == 1) {
			x <- try( setValues(x, rep(value, times=ncell(x))) )
		} else {
			v <- try( vector(length=ncell(x)) )
			if (class(x) != 'try-error') {
				v[] <- value
				x <- try( setValues(x, v) )
			}
		}
		if (class(x) == 'try-error') {
			stop('cannot replace values on this raster (it is too large')
		}
		return(x)
	
	}
)


.replace <- function(x, i, value, recycle=1) {
	
	if ( is.logical(i) ) {
		i <- which(i)
	} else {
	#	if (! is.numeric(i)) { 
	#		i <- as.integer(i) 
	#	}
		i <- stats::na.omit(i)
	}

	nl <- nlayers(x)
  # recycling
	if (nl > 1 & recycle > 0) {
		rec2 <- ceiling(nl / recycle)
		if (rec2 > 1) {
			add <- ncell(x)*recycle * (0:(rec2-1))
			i <- as.vector(t((matrix(rep(i, rec2), nrow=rec2, byrow=TRUE)) + add))
		}
	}
	j <- i > 0 & i <= (ncell(x)*nl)
	
	if (!all(j)) {
		i <- i[j]
		if (length(value) > 1) {
			value <- value[j]
		}
	}


	if ( inMemory(x) ) {
		if (inherits(x, 'RasterStack')) {
			x <- brick( x, values=TRUE )
			# this may go to disk, hence we check again below
		}	
	}
	
	if ( inMemory(x) ) {
		x@data@values[i] <- value
		x <- setMinMax(x)
		x <- .clearFile(x)
		return(x)
		
	} else if (canProcessInMemory(x)) {
		if (inherits(x, 'RasterStack')) {
			x <- brick( x, values=TRUE )
			if (!inMemory(x)) {
				x <- readAll(x) 
			}
			x <- .clearFile(x)
		} else if ( fromDisk(x) ) {
			x <- readAll(x)
			x <- .clearFile(x)
		} else {
			x <- setValues(x, rep(NA, times=ncell(x)))
		}
		x@data@values[i] <- value
		x <- setMinMax(x)
		return(x)
			
	} else {
	
		tr <- blockSize(x)
		pb <- pbCreate(tr$n, label='replace')
		hv <- hasValues(x)
		if (nl==1) {
			if (! length(value) %in% c(1, length(i))) {
				stop('cannot replace values in large Raster objects if their length is not 1 or the number of cells to be replaced')
			}
			r <- raster(x)
			r <- writeStart(r, filename=rasterTmpFile(), overwrite=TRUE )
			for (k in 1:tr$n) {
				# cells <- cellFromRowCol(x, tr$row[k], 1):cellFromRowCol(x, tr$row[k]+tr$nrows[k]-1, ncol(x))
				cell1 <- cellFromRowCol(x, tr$row[k], 1)
				cell2 <- cell1 + tr$nrows[k] * ncol(x) - 1
				if (hv) {
					v <- getValues(x, row=tr$row[k], nrows=tr$nrows[k])
				} else {
					v <- rep(NA, 1+cell2-cell1)
				}
				j <- which(i >= cell1 & i <= cell2)
				if (length(j) > 0) {
					localcells <- i[j] - (cell1-1)
					if (length(value) == length(i)) {
						v[localcells] <- value[j]
					} else {
						v[localcells] <- value
					}
				}
				r <- writeValues(r, v, tr$row[k])
				pbStep(pb, k) 	
			}
			r <- writeStop(r)
			pbClose(pb)
			return(r)
				
		} else {
			if (! length(value) %in% c(1, length(i))) {
				stop('length of replacement values does not match the length of the index')
			}
			r <- brick(x, values=FALSE)
			r <- writeStart(r, filename=rasterTmpFile(), overwrite=TRUE )
			add <- (0:(nl-1)) * ncell(x)
			for (k in 1:tr$n) {
				cells <- cellFromRowCol(x, tr$row[k], 1):cellFromRowCol(x, tr$row[k]+tr$nrows[k]-1, ncol(x))
				if (hv) {
					v <- getValues(x, row=tr$row[k], nrows=tr$nrows[k])
				} else {
					v <- matrix(NA, nrow=length(cells), ncol=nl)
				}
				cells <- cells + rep(add, each=length(cells))
				j <- cells %in% i
				if (sum(j) > 0) {
					v[j] <- value
				}
				r <- writeValues(r, v, tr$row[k])
				pbStep(pb, k)
			}
			r <- writeStop(r)
			pbClose(pb)
			return(r)
		}	
	}
}