[Raster-commits] r380 - in pkg/raster: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 21 15:43:37 CET 2009
Author: rhijmans
Date: 2009-03-21 15:43:36 +0100 (Sat, 21 Mar 2009)
New Revision: 380
Added:
pkg/raster/R/focalFilter.R
pkg/raster/R/focalFun.R
pkg/raster/R/zonal.R
pkg/raster/man/focal.Rd
pkg/raster/man/zonal.Rd
Removed:
pkg/raster/R/neighborhood.R
pkg/raster/man/neighborhood.Rd
Modified:
pkg/raster/DESCRIPTION
pkg/raster/R/plot.R
pkg/raster/R/readRaster.R
Log:
Modified: pkg/raster/DESCRIPTION
===================================================================
--- pkg/raster/DESCRIPTION 2009-03-20 16:25:53 UTC (rev 379)
+++ pkg/raster/DESCRIPTION 2009-03-21 14:43:36 UTC (rev 380)
@@ -1,8 +1,8 @@
Package: raster
Type: Package
Title: Raster data handling for geographic data analysis and modeling
-Version: 0.8.9-9
-Date: 20-March-2009
+Version: 0.8.9-10
+Date: 21-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>
Added: pkg/raster/R/focalFilter.R
===================================================================
--- pkg/raster/R/focalFilter.R (rev 0)
+++ pkg/raster/R/focalFilter.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -0,0 +1,87 @@
+# Author: Robert J. Hijmans, r.hijmans at gmail.com
+# International Rice Research Institute
+# Date : March 2009
+# Version 0.8
+# Licence GPL v3
+
+
+
+.calcFilter <- function(rows, colnrs, res, filter) {
+ res[] <- NA
+ for (i in 1:dim(rows)[2]) {
+ d <- rows[, colnrs[i, ]]
+ if (!all(dim(d) == dim(filter))) {
+ res[i] <- NA
+ } else {
+ res[i] <- sum(d * filter)
+ }
+ }
+ return(res)
+}
+
+
+focalFilter <- function(raster, filter, filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
+ if (!is.matrix(filter)) {stop('filter must be a matrix')}
+ ngb <- dim(filter)
+
+ filename <- trim(filename)
+ ngbgrid <- raster(raster, filename=filename)
+ dataType(ngbgrid) <- datatype
+
+ res <- vector(length=length(ncol(ngbgrid)))
+
+ limcol <- floor(ngb[2] / 2)
+ colnrs <- (-limcol+1):(ncol(ngbgrid)+limcol)
+ colnrs <- .embed(colnrs, ngb[2])
+ colnrs[colnrs > ncol(ngbgrid) | colnrs < 0] <- 0
+
+ limrow <- floor(ngb[1] / 2)
+ 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:limrow) {
+ if (dataContent(raster)=='all') {
+ rowdata <- valuesRow(raster, r)
+ } else {
+ rowdata <- values(readRow(raster, r))
+ }
+ ngbdata <- rbind(ngbdata, rowdata)
+ }
+
+ res <- vector(length=ncol(ngbdata))
+
+ v <- vector(length=0)
+ starttime <- proc.time()
+
+ for (r in 1:nrow(ngbgrid)) {
+ rr <- r + limrow
+ if (rr <= nrow(ngbgrid)) {
+ if (dataContent(raster)=='all') {
+ rowdata <- valuesRow(raster, rr)
+ } else {
+ rowdata <- values(readRow(raster, rr))
+ }
+ if (dim(ngbdata)[1] == ngb[1]) {
+ ngbdata <- rbind(ngbdata[2:ngb[1],], rowdata)
+ } else {
+ ngbdata <- rbind(ngbdata, rowdata)
+ }
+ } else {
+ ngbdata <- ngbdata[-1, ,drop=FALSE]
+ }
+
+
+ ngbvals <- .calcFilter(ngbdata, colnrs, res, filter)
+ if (filename != "") {
+ ngbgrid <- setValues(ngbgrid, ngbvals, r)
+ ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
+ } else {
+ v <- c(v, ngbvals)
+ }
+ if (r %in% track) { .showTrack(r, ngbgrid at nrows, track, starttime) }
+ }
+ if (filename == "") {
+ ngbgrid <- setValues(ngbgrid, v)
+ }
+ return(ngbgrid)
+}
+
Added: pkg/raster/R/focalFun.R
===================================================================
--- pkg/raster/R/focalFun.R (rev 0)
+++ pkg/raster/R/focalFun.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -0,0 +1,100 @@
+# Author: Robert J. Hijmans, r.hijmans at gmail.com
+# International Rice Research Institute
+# Date : June 2008
+# Version 0.8
+# Licence GPL v3
+
+.embed <- function(x, dimension) {
+ n <- length(x)
+ m <- n - dimension + 1
+ data <- x[1:m + rep.int(1:dimension, rep.int(m, dimension)) - 1]
+ dim(data) <- c(m, dimension)
+ return(data)
+}
+
+
+.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)
+}
+
+
+focal <- function(raster, fun=mean, filename="", ngb=3, keepdata=TRUE, overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
+ ngb <- as.integer(round(ngb))
+ if (length(ngb) == 1) {
+ ngb <- c(ngb, ngb)
+ } else if (length(ngb) > 2) {
+ stop('ngb should be a single value or two values')
+ }
+ if (min(ngb) < 3) { stop("ngb should be 3 or larger") }
+
+# if (ngb[1] %% 2 == 0 | ngb[2] %% 2 == 0) { stop("only odd neighborhoods are supported") }
+
+ filename <- trim(filename)
+ ngbgrid <- raster(raster, filename=filename)
+ dataType(ngbgrid) <- datatype
+
+# first create an empty matrix with nrows = ngb and ncols = raster at ncols
+ res <- vector(length=length(ncol(ngbgrid)))
+ limcol <- floor(ngb[2] / 2)
+ colnrs <- (-limcol+1):(ncol(ngbgrid)+limcol)
+ colnrs <- .embed(colnrs, ngb[2])
+ colnrs[colnrs > ncol(ngbgrid) | colnrs < 0] <- 0
+
+ limrow <- floor(ngb[1] / 2)
+ 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:limrow) {
+ if (dataContent(raster)=='all') {
+ rowdata <- valuesRow(raster, r)
+ } else {
+ rowdata <- values(readRow(raster, r))
+ }
+ ngbdata <- rbind(ngbdata, rowdata)
+ }
+
+ res <- vector(length=ncol(ngbdata))
+
+ v <- vector(length=0)
+ starttime <- proc.time()
+
+ for (r in 1:nrow(ngbgrid)) {
+ rr <- r + limrow
+ if (rr <= nrow(ngbgrid)) {
+ if (dataContent(raster)=='all') {
+ rowdata <- valuesRow(raster, rr)
+ } else {
+ rowdata <- values(readRow(raster, rr))
+ }
+ if (dim(ngbdata)[1] == ngb[1]) {
+ ngbdata <- rbind(ngbdata[2:ngb[1],], rowdata)
+ } else {
+ ngbdata <- rbind(ngbdata, rowdata)
+ }
+ } else {
+ ngbdata <- ngbdata[-1, ,drop=FALSE]
+ }
+
+
+ ngbvals <- .calcNGB(ngbdata, colnrs, res, fun, keepdata)
+ if (filename != "") {
+ ngbgrid <- setValues(ngbgrid, ngbvals, r)
+ ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
+ } else {
+ v <- c(v, ngbvals)
+ }
+ if (r %in% track) { .showTrack(r, ngbgrid at nrows, track, starttime) }
+ }
+ if (filename == "") {
+ ngbgrid <- setValues(ngbgrid, v)
+ }
+ return(ngbgrid)
+}
+
Deleted: pkg/raster/R/neighborhood.R
===================================================================
--- pkg/raster/R/neighborhood.R 2009-03-20 16:25:53 UTC (rev 379)
+++ pkg/raster/R/neighborhood.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -1,87 +0,0 @@
-# Author: Robert J. Hijmans, r.hijmans at gmail.com
-# International Rice Research Institute
-# Date : June 2008
-# Version 0.8
-# Licence GPL v3
-
-
-.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)
-}
-
-
-
-neighborhood <- function(raster, fun=mean, filename="", ngb=3, keepdata=TRUE, overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
- 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=filename)
- dataType(ngbgrid) <- datatype
-
-# first create an empty matrix with nrows = ngb and ncols = raster at ncols
- 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
-
- 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, 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 {
- ngbdata <- rbind(ngbdata, rowdata)
- }
- } else {
- ngbdata <- ngbdata[-1, ,drop=FALSE]
- }
-
-
- ngbvals <- .calcNGB(ngbdata, colnrs, res, fun, keepdata)
- if (filename != "") {
- ngbgrid <- setValues(ngbgrid, ngbvals, r)
- ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
- } else {
- v <- c(v, ngbvals)
- }
-# if (r %in% track) { .showTrack(r, ngbgrid at nrows, track, starttime) }
- }
- if (filename == "") {
- ngbgrid <- setValues(ngbgrid, v)
- }
- return(ngbgrid)
-}
-
Modified: pkg/raster/R/plot.R
===================================================================
--- pkg/raster/R/plot.R 2009-03-20 16:25:53 UTC (rev 379)
+++ pkg/raster/R/plot.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -53,7 +53,7 @@
if (length(x) < nc) {
warning(paste('plot used a sample of ', round(100*length(x)/nc), "% of the cells", sep=""))
}
- plot(x, y, ...)
+ plot(x, y, cex=cex, ...)
}
)
Modified: pkg/raster/R/readRaster.R
===================================================================
--- pkg/raster/R/readRaster.R 2009-03-20 16:25:53 UTC (rev 379)
+++ pkg/raster/R/readRaster.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -23,7 +23,14 @@
}
if (dataSource(raster)=='ram') {
- result <- valuesRow(raster, rownr)[startcol:endcol]
+ if (rownr < 1) {
+ if (dataContent(raster) != 'all') {
+ stop('cannot read data for this RasterLayer')
+ }
+ return(raster)
+ } else {
+ result <- valuesRow(raster, rownr)[startcol:endcol]
+ }
} else if (.driver(raster) == 'raster') {
Added: pkg/raster/R/zonal.R
===================================================================
--- pkg/raster/R/zonal.R (rev 0)
+++ pkg/raster/R/zonal.R 2009-03-21 14:43:36 UTC (rev 380)
@@ -0,0 +1,69 @@
+# Author: Robert J. Hijmans, r.hijmans at gmail.com
+# International Rice Research Institute
+# Date : March 2009
+# Version 0.8
+# Licence GPL v3
+
+zonal <- function(raster, zones, stat='mean', keepdata=TRUE, track=-1) {
+
+ if (class(stat) != 'character') {
+ if (canProcessInMemory(raster, 4)) {
+ d <- cbind(values(readAll(raster)), as.integer(values(readAll(zones))))
+ if (keepdata) {
+ d <- na.omit(d)
+ }
+ alltab <- tapply(d[,1], d[,2], stat)
+ stat <- deparse(substitute(stat))
+ } else {
+ stop("RasterLayers are too large. You an use fun='sum', 'mean', 'min', or 'max', but not a function")
+ }
+ } else {
+
+ counts <- FALSE
+ if (stat == 'sum') {
+ fun <- sum
+ } else if (stat == 'min') {
+ fun <- min
+ } else if (stat == 'max') {
+ fun <- max
+ } else if (stat == 'mean') {
+ fun <- sum
+ counts <- TRUE
+ } else {
+ stop("invalid 'stat', should be 'sum', 'min', 'max', or 'mean'")
+ }
+
+ alltab <- array(dim=0)
+ cnttab <- alltab
+ starttime <- proc.time()
+ for (r in 1:nrow(raster)) {
+ d <- cbind(valuesRow(raster, r), as.integer(valuesRow(zones, r)))
+ if (keepdata) {
+ d <- na.omit(d)
+ }
+ alltab <- c(alltab, tapply(d[,1], d[,2], fun))
+ if (counts) {
+ cnttab <- c(cnttab, tapply(d[,1], d[,2], length))
+ }
+ if (length(alltab) > 10000) {
+ groups <- as.integer(names(alltab))
+ alltab <- tapply(as.vector(alltab), groups, fun)
+ if (counts) {
+ cnttab <- tapply(as.vector(cnttab), groups, sum)
+ }
+ }
+ if (r %in% track) { .showTrack(r, raster at nrows, track, starttime) }
+ }
+ groups <- as.integer(names(alltab))
+ alltab <- tapply(as.vector(alltab), groups, fun)
+ if (counts) {
+ cnttab <- tapply(as.vector(cnttab), groups, sum)
+ alltab <- alltab / cnttab
+ }
+ }
+ zone <- as.integer(names(alltab))
+ alltab <- data.frame(zone, alltab)
+ colnames(alltab) <- c('zone', stat)
+ return(alltab)
+}
+
Added: pkg/raster/man/focal.Rd
===================================================================
--- pkg/raster/man/focal.Rd (rev 0)
+++ pkg/raster/man/focal.Rd 2009-03-21 14:43:36 UTC (rev 380)
@@ -0,0 +1,70 @@
+\name{focal}
+
+\alias{focal}
+\alias{focalFilter}
+
+\title{Focal}
+
+\description{
+Calculate values for the neighborhood of focal cells
+}
+
+\usage{
+focal(raster, fun=mean, filename="", ngb=3, keepdata=TRUE, overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
+focalFilter(raster, filter, filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
+}
+
+\arguments{
+ \item{raster}{A RasterLayer object}
+ \item{fun}{The function to be applied}
+ \item{filename}{Output filename for a new raster}
+ \item{ngb}{Neighborhood size. See Details.}
+ \item{keepdata}{Logical. If \code{TRUE}, NA will be removed for neighborhood computations. The result will only be NA if all cells are NA}
+ \item{overwrite}{Logical to indicate whether an existing output file should be overwritten}
+ \item{filetype}{output file type. Either 'raster', 'ascii' or a supported GDAL 'driver' name see \code{\link[raster]{writeRaster}}}
+ \item{datatype}{output data type; see \code{\link[raster]{setDatatype}}}
+ \item{track}{vector of row numbers for which the function will report that they have been processed}
+ \item{filter}{a matrix of weights. See Details}
+}
+
+\details{
+\code{focal} computes values based on those in the neighborhood of each cell. The values are computed with function 'fun',
+using the values of the cells in the neigborhood of a focal cell (including the focal cell).
+
+A neighborhood is expressed in number of cells in a single direction or in two direction from the focal cell. I.e. \code{ngb=3} refers to 2 cells
+at each side of the focal cell. Queen's case, 8 cells in total. This is equivalent to \code{ngb=c(3,3)}.
+You can also specify a rectangular neighborhood, e.g. \code{ngb=c(3,5)}. Neighborhoods are normally odd numbers so that the focal
+cell is centered, but even numbers are permitted in this function.
+
+\code{focalFilter} does not use a function but a matrix of weights for the neigbhorhood of the focal cells.
+For example, \code{filter=matrix(1/9, nrow=3, ncol=3)} would be equivalent to \code{mean} with \code{ngb=3} in the \code{focal} function.
+
+Gaussian filter:
+\code{filter=matrix(c(1,2,3,2,1,2,3,4,3,2,3,4,5,4,3,2,3,4,3,2,1,2,3,2,1), nrow=5)/65}
+
+Laplacian filter:
+\code{filter=matrix(c(0,1,0,1,-4,1,0,1,0), nrow=3)}
+
+Sobel filter:
+\code{filter=matrix(c(1,2,1,0,0,0,-1,-2,-1) / 4, nrow=3)}
+
+}
+
+\value{
+A new RasterLayer object (in the R environment), and in some cases the side effect of a new file on disk.
+}
+
+\author{Robert J. Hijmans}
+
+\examples{
+
+r <- raster(ncols=36, nrows=18)
+r[] <- runif(ncell(r))
+rf <- focal(r, fun=sum, ngb=3)
+
+gf=matrix(c(1,2,3,2,1,2,3,4,3,2,3,4,5,4,3,2,3,4,3,2,1,2,3,2,1), nrow=5)/65
+rff <- focalFilter(r, filter = gf)
+
+}
+\keyword{spatial}
+
Deleted: pkg/raster/man/neighborhood.Rd
===================================================================
--- pkg/raster/man/neighborhood.Rd 2009-03-20 16:25:53 UTC (rev 379)
+++ pkg/raster/man/neighborhood.Rd 2009-03-21 14:43:36 UTC (rev 380)
@@ -1,47 +0,0 @@
-\name{neighborhood}
-
-\alias{neighborhood}
-
-\title{RasterLayer calculate}
-
-\description{
-Neighborhood calculations
-}
-
-\usage{
-neighborhood(raster, fun=mean, filename="", ngb=3, keepdata=TRUE, overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
-}
-
-\arguments{
- \item{raster}{A RasterLayer object}
- \item{fun}{The function to be applied}
- \item{filename}{Output filename for a new raster}
- \item{ngb}{Neighborhood size. Expressed in number of cells in a single direction. I.e. 3 refers to 2 cells at each side of the focal cell. Queen's case, 8 cells in total}
- \item{keepdata}{Logical. If \code{TRUE}, NA will be removed for neighborhood computations. The result will only be NA if all cells are NA}
- \item{overwrite}{Logical to indicate whether an existing output file should be overwritten}
- \item{filetype}{output file type. Either 'raster', 'ascii' or a supported GDAL 'driver' name see \code{\link[raster]{writeRaster}}}
- \item{datatype}{output data type; see \code{\link[raster]{setDatatype}}}
- \item{track}{vector of row numbers for which the function will report that they have been processed}
-}
-
-\details{
-\code{neighborhood} computes a neighborhood value. I.e. a value computed with function 'fun', for all cells in the square neigborhood (of size ngb * ngb) around each cell
-If the input RasterLayer object has all values in memory (e.g. after readAll(raster)), the function will also return the new values in memory. If a filename is provided, the values will also be saved to that file.
-If the values are not in memory the new values will be written to file.
-}
-
-\value{
-A new RasterLayer object (in the R environment), and in some cases the side effect of a new file on disk.
-}
-
-\author{Robert J. Hijmans}
-
-\examples{
-
-r <- raster(ncols=36, nrows=18)
-r[] <- runif(ncell(r))
-ngb <- neighborhood(r, fun=sum)
-
-}
-\keyword{spatial}
-
Added: pkg/raster/man/zonal.Rd
===================================================================
--- pkg/raster/man/zonal.Rd (rev 0)
+++ pkg/raster/man/zonal.Rd 2009-03-21 14:43:36 UTC (rev 380)
@@ -0,0 +1,45 @@
+\name{zonal}
+
+\alias{zonal}
+
+\title{Zonal statistics}
+
+\description{
+Compute zonal statistics. I.e. the values of the cells of a RasterLayer based on the values of the "zones" RasterLayer.
+}
+
+\usage{
+zonal(raster, zones, stat='mean', keepdata=TRUE, track=-1)
+}
+
+\arguments{
+ \item{raster}{A RasterLayer object}
+ \item{zones}{A RasterLayer object; with codes representing zones that are treated as integer values}
+ \item{stat}{The function to be applied. Either as character: 'mean', 'min', 'max', 'sum'; or a function (see Details) }
+ \item{keepdata}{Logical. If \code{TRUE}, NA will be removed for neighborhood computations. The result will only be NA if all cells are NA}
+ \item{track}{vector of row numbers for which the function will report that they have been processed}
+}
+
+\details{
+If \code{stat} is a \code{function}, \code{zonal} will fail for very large RasterLayers
+}
+
+\value{
+A data frame with a value for each zone
+}
+
+\author{Robert J. Hijmans}
+
+\examples{
+r <- raster(ncols=10, nrows=10)
+r[] <- runif(ncell(r))
+z <- r
+z[] <- rep(1:5, each=20)
+zonal(r, z, 'mean')
+zonal(r, z, 'min')
+zonal(r, z, 'sum')
+# using a function, rather than a character value
+zonal(r, z, sum)
+}
+
+\keyword{spatial}
More information about the Raster-commits
mailing list