[Raster-commits] r378 - in pkg/raster: . R

[noreply at r-forge.r-project.org]

Author: rhijmans
Date: 2009-03-20 16:04:16 +0100 (Fri, 20 Mar 2009)
New Revision: 378

Modified:
   pkg/raster/DESCRIPTION
   pkg/raster/R/neighborhood.R
Log:


Modified: pkg/raster/DESCRIPTION
===================================================================
--- pkg/raster/DESCRIPTION	2009-03-20 09:41:38 UTC (rev 377)
+++ pkg/raster/DESCRIPTION	2009-03-20 15:04:16 UTC (rev 378)
@@ -1,8 +1,8 @@
 Package: raster
 Type: Package
 Title: Raster data handling for geographic data analysis and modeling
-Version: 0.8.9-8
-Date: 19-March-2009
+Version: 0.8.9-9
+Date: 20-March-2009
 Depends: methods, sp, rgdal (>= 0.5-33), R (>= 2.8.0)
 Author: Robert J. Hijmans & Jacob van Etten
 Maintainer: Robert J. Hijmans <r.hijmans at gmail.com> 

Modified: pkg/raster/R/neighborhood.R
===================================================================
--- pkg/raster/R/neighborhood.R	2009-03-20 09:41:38 UTC (rev 377)
+++ pkg/raster/R/neighborhood.R	2009-03-20 15:04:16 UTC (rev 378)
@@ -1,97 +1,83 @@
 # Author: Robert J. Hijmans, r.hijmans at gmail.com
 # International Rice Research Institute
 # Date :  June 2008
-# Version 0,7
+# Version 0.8
 # Licence GPL v3
 
 
-.calc.ngb <- function(rows, ngb, fun, keepdata) {
-	lim <- floor(ngb / 2)
-	res <- vector(length=length(rows[1,]))
-	lr <- length(rows[1,])
-    for (i in 1:length(rows[1,])) {
-		d <- rows[, max(1,(i-lim)):min((i+lim),lr)]
-		d <- as.vector(d)
-		dd <- as.vector(na.omit(d))
-		if (length(dd) == 0) {
-			res[i] <- NA
-		} else if (keepdata) { 
-			res[i] <- fun(dd)
-		} else {
-			if (length(dd) == length(d)) { 
-				res[i] <- fun(d)
-			} else {
-				res[i] <- NA
-			}
+.calcNGB <- function(rows, colnrs, res, fun, keepdata) {
+	res[] <- NA
+    for (i in 1:dim(rows)[2]) {
+		d <- as.vector(rows[, colnrs[i, ]])
+		if (keepdata) {
+			d <- na.omit(d)
 		}
+		res[i] <- fun(d)
 	}	
 	return(res)
 }
 
 
-.calc.ngb2 <- function(rows, ngb, fun, keepdata) {
-#TODO, current function is very slow
-	lim <- floor(ngb / 2)
-	res <- array(dim=length(rows[1,]))
-	addNA <- (matrix(ncol=lim, nrow=ngb))
-	rows <- cbind(addNA, rows, addNA)
-#	d <- rows[, max(1,(i-lim)):min((i+lim),lr)]
-#	if (rm.NA) {outraster at data@values <- as.vector(tapply(raster at data,cell.index,function(x){fun(na.omit(x))}))}
-#		else {outraster at data@values <- as.vector(tapply(raster at data,cell.index,fun))}
-}
-	
 
 neighborhood <- function(raster, fun=mean, filename="", ngb=3, keepdata=TRUE, overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
-	ngb <- round(ngb)
-	if ((ngb / 2) == floor(ngb/2)) { stop("only odd neighborhoods are supported") }
-	if (ngb == 1) { stop("ngb should be 3 or larger")  } 
+	ngb <- as.integer(round(ngb))
+	if (ngb < 3) { stop("ngb should be 3 or larger") } 
+	if (ngb %% 2 == 0) { stop("only odd neighborhoods are supported") }
 	lim <- floor(ngb / 2)
-	
+
 	filename <- trim(filename)
-	ngbgrid <- raster(raster, filename)
+	ngbgrid <- raster(raster, filename=filename)
 	dataType(ngbgrid) <- datatype
 
 # first create an empty matrix with nrows = ngb and ncols = raster at ncols
-	ngbdata1 <- array(data = NA, dim = c(ngb, ncol(raster)))
-	ngbdata <- ngbdata1
-	
-	rr <- 1
-	v <- vector(length=0)
-	starttime <- proc.time()
+	res <- vector(length=length(ncol(ngbgrid)))
+	lim <- floor(ngb / 2)
+	colnrs <- (-lim+1):(ncol(ngbgrid)+lim)
+	colnrs <- embed(colnrs, ngb)
+	colnrs[colnrs > ncol(ngbgrid) | colnrs < 0] <- 0
 
-	for (r in 1:nrow(raster)) {		
+	ngbdata <- matrix(NA, nrow=0, ncol=ncol(ngbgrid))
+# add all rows needed for first ngb, minus 1 that will be read in first loop	
+	for (r in 1:lim) {
 		if (dataContent(raster)=='all') {
 			rowdata <- valuesRow(raster, r)
 		} else {	
 			rowdata <- values(readRow(raster, r))
-		}	
-		ngbdata <- rbind(ngbdata[2:ngb,], t(rowdata))
-		if (r > lim) {
-			ngbvals <- .calc.ngb(ngbdata, ngb, fun, keepdata)
-			if (filename != "") {
-				ngbgrid <- setValues(ngbgrid, ngbvals, rr)
-				ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
+		}
+		ngbdata <- rbind(ngbdata, rowdata)
+	}
+
+	res <- vector(length=ncol(ngbdata))
+
+	v <- vector(length=0)
+	starttime <- proc.time()
+
+	for (r in 1:nrow(ngbgrid)) {		
+		rr <- r + lim
+		if (rr <= nrow(ngbgrid)) {
+			if (dataContent(raster)=='all') {
+				rowdata <- valuesRow(raster, r)
+			} else {	
+				rowdata <- values(readRow(raster, r))
+			}
+			if (dim(ngbdata)[1] == ngb) {
+				ngbdata <- rbind(ngbdata[2:ngb,], rowdata)
 			} else {
-				v <- c(v, ngbvals)
+				ngbdata <- rbind(ngbdata, rowdata)			
 			}
-			rr <- rr + 1
+		} else {
+			ngbdata <- ngbdata[-1, ,drop=FALSE]
 		}
 
-		if (r %in% track) { .showTrack(r, raster at nrows, track, starttime) }
-
-	}
-	
-	ngbdata1 <- array(data = NA, dim = c(ngb, raster at ncols))
-	for (r in (nrow(raster)+1):(nrow(raster)+lim)) {
-		ngbdata <- rbind(ngbdata[2:ngb,], t(ngbdata1[1,]))
-		ngbvals <- .calc.ngb(ngbdata, ngb, fun, keepdata)
+		
+		ngbvals <- .calcNGB(ngbdata, colnrs, res, fun, keepdata)
 		if (filename != "") {
-			ngbgrid <- setValues(ngbgrid, ngbvals, rr)
+			ngbgrid <- setValues(ngbgrid, ngbvals, r)
 			ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
 		} else {
 			v <- c(v, ngbvals)
 		}
-		rr <- rr + 1
+#		if (r %in% track) { .showTrack(r, ngbgrid at nrows, track, starttime) }
 	}
 	if (filename == "") { 
 		ngbgrid <- setValues(ngbgrid, v)