[Raster-commits] r398 - in pkg/raster: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 6 18:16:34 CEST 2009
Author: rhijmans
Date: 2009-04-06 18:16:33 +0200 (Mon, 06 Apr 2009)
New Revision: 398
Added:
pkg/raster/man/focalFilter.Rd
Modified:
pkg/raster/R/focalFilter.R
pkg/raster/R/polygonToRaster.R
pkg/raster/man/LinesToRaster.Rd
pkg/raster/man/PolygonsToRaster.Rd
pkg/raster/man/focal.Rd
pkg/raster/man/setDatatype.Rd
Log:
Modified: pkg/raster/R/focalFilter.R
===================================================================
--- pkg/raster/R/focalFilter.R 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/R/focalFilter.R 2009-04-06 16:16:33 UTC (rev 398)
@@ -6,21 +6,21 @@
-.calcFilter <- function(rows, colnrs, res, filter) {
+.calcFilter <- function(rows, colnrs, res, filter, fun) {
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)
+ res[i] <- fun(d * filter)
}
}
return(res)
}
-focalFilter <- function(raster, filter, filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
+focalFilter <- function(raster, filter, fun=sum, filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
if (!is.matrix(filter)) {stop('filter must be a matrix')}
ngb <- dim(filter)
@@ -70,7 +70,7 @@
}
- ngbvals <- .calcFilter(ngbdata, colnrs, res, filter)
+ ngbvals <- .calcFilter(ngbdata, colnrs, res, filter, fun)
if (filename != "") {
ngbgrid <- setValues(ngbgrid, ngbvals, r)
ngbgrid <- writeRaster(ngbgrid, overwrite=overwrite, filetype=filetype)
Modified: pkg/raster/R/polygonToRaster.R
===================================================================
--- pkg/raster/R/polygonToRaster.R 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/R/polygonToRaster.R 2009-04-06 16:16:33 UTC (rev 398)
@@ -60,9 +60,14 @@
-polygonsToRaster <- function(spPolys, raster, field=0, updateRaster=FALSE, updateValue="NA", filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
+polygonsToRaster <- function(spPolys, raster, field=0, overlap='last', updateRaster=FALSE, updateValue="NA",
+ filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1) {
+
filename <- trim(filename)
+ if (!(overlap %in% c('first', 'last', 'sum', 'min', 'max', 'count'))) {
+ stop('invalid value for overlap')
+ }
if (updateRaster) {
oldraster <- raster
if (!(updateValue == 'NA' | updateValue == '!NA' | updateValue == 'all' | updateValue == 'zero')) {
@@ -80,11 +85,15 @@
if (spbb[1,1] > rsbb[1,2] | spbb[2,1] > rsbb[2,2]) {
stop('polygon and raster have no overlapping areas')
}
+
npol <- length(spPolys at polygons)
- if (class(spPolys) == 'SpatialPolygons' | field == 0) {
+
+ if (field < 0) {
+ putvals <- rep(1, length=npol)
+ } else if (class(spPolys) == 'SpatialPolygons' | field == 0) {
putvals <- as.integer(1:npol)
} else {
- putvals <- as.vector(spPolys at data[,field])
+ putvals <- as.vector(spPolys at data[[field]])
if (class(putvals) == 'character') {
stop('selected field is charater type')
}
@@ -144,6 +153,7 @@
intersection <- .intersectLinePolygon(myline, mypoly at coords)
x <- sort(intersection[,1])
if (length(x) > 0) {
+ rvtmp <- rv1
if ( sum(x[-length(x)] == x[-1]) > 0 ) {
# single node intersection going out of polygon ....
spPnts <- xyFromCell(raster, cellFromRowCol(raster, rep(r, ncol(raster)), 1:ncol(raster)), TRUE)
@@ -152,7 +162,7 @@
if ( subpol[i, 4] == 1 ) {
holes[over] <- TRUE
} else {
- rv[over] <- subpol[i,4]
+ rvtmp[over] <- subpol[i,4]
}
# print(paste('exit node intersection on row:', r))
} else {
@@ -177,10 +187,24 @@
if ( subpol[i, 5] == 1 ) {
holes[col1:col2] <- TRUE
} else {
- rv[col1:col2] <- subpol[i,4]
+ rvtmp[col1:col2] <- subpol[i,4]
}
}
}
+ if (overlap=='last') {
+ rv[!is.na(rvtmp)] <- rvtmp[!is.na(rvtmp)]
+ } else if (overlap=='first') {
+ rv[is.na(rv)] <- rvtmp[is.na(rv)]
+ } else if (overlap=='sum') {
+ rv[!is.na(rv) & !is.na(rvtmp)] <- rv[!is.na(rv) & !is.na(rvtmp)] + rvtmp[!is.na(rv) & !is.na(rvtmp)]
+ rv[is.na(rv)] <- rvtmp[is.na(rv)]
+ } else if (overlap=='min') {
+ rv[!is.na(rv) & !is.na(rvtmp)] <- pmin(rv[!is.na(rv) & !is.na(rvtmp)], rvtmp[!is.na(rv) & !is.na(rvtmp)])
+ rv[is.na(rv)] <- rvtmp[is.na(rv)]
+ } else if (overlap=='max') {
+ rv[!is.na(rv) & !is.na(rvtmp)] <- pmax(rv[!is.na(rv) & !is.na(rvtmp)], rvtmp[!is.na(rv) & !is.na(rvtmp)])
+ rv[is.na(rv)] <- rvtmp[is.na(rv)]
+ }
}
}
}
Modified: pkg/raster/man/LinesToRaster.Rd
===================================================================
--- pkg/raster/man/LinesToRaster.Rd 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/man/LinesToRaster.Rd 2009-04-06 16:16:33 UTC (rev 398)
@@ -2,7 +2,7 @@
\alias{linesToRaster}
-\title{Transfer lines to a RasterLayer }
+\title{Lines to raster}
\description{
Lines to raster conversion. Lines are raserized.
Modified: pkg/raster/man/PolygonsToRaster.Rd
===================================================================
--- pkg/raster/man/PolygonsToRaster.Rd 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/man/PolygonsToRaster.Rd 2009-04-06 16:16:33 UTC (rev 398)
@@ -10,15 +10,16 @@
}
\usage{
-polygonsToRaster(spPolys, raster, field=0, updateRaster=FALSE, updateValue="NA", filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
+polygonsToRaster(spPolys, raster, field=0, overlap='last', updateRaster=FALSE, updateValue="NA", filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
}
\arguments{
- \item{spPolys}{ a SpatialPolygons or a SpatialPolygonsDataFrame object }
- \item{raster}{ RasterLayer }
- \item{field}{ The index of the column in the SpatialPolygonsDataFrame to be transfered to the RasterLayer }
- \item{filename}{ output filename }
- \item{overwrite}{ logical. if \code{TRUE} ouput file will be overwritten if it exists }
+ \item{spPolys}{a SpatialPolygons or a SpatialPolygonsDataFrame object }
+ \item{raster}{a RasterLayer}
+ \item{field}{The index (or colname) of the column in the SpatialPolygonsDataFrame to be transfered to the RasterLayer (see Details)}
+ \item{overlap}{character. Determines what to do with cells with overlapping polygons. Either 'first', 'last', 'sum', 'min', or 'max'}
+ \item{filename}{output filename }
+ \item{overwrite}{logical. if \code{TRUE} ouput file will be overwritten if it exists }
\item{updateRaster}{logical. If \code{TRUE} the values of the input RasterLayer are updated where the polygons overlap cells }
\item{updateValue}{character. Select cells to be updated (if \code{updateRaster == TRUE}) by their current values. Either \code{'all'}, \code{'NA'}, \code{'!NA'}, or \code{'zero'} }
\item{filetype}{output file type. Either 'raster', 'ascii' or a supported GDAL 'driver' name see \code{\link[raster]{writeRaster}}}
@@ -27,7 +28,8 @@
}
\details{
-For SpatialPolygons, the 'field' argument is ignored. The polygon index is tranferred.
+For SpatialPolygons, the \code{field} argument is ignored. The polygon index is used (i.e. numbers from 1 to the number of polygons. This is also
+done for SpatialPolygonDataFrame objects when \code{field}==0. If \code{field} < 0, all polygons get the value 1.
}
\author{Robert J. Hijmans}
@@ -35,18 +37,20 @@
\seealso{ \code{\link[raster]{linesToRaster}}, \code{\link[raster]{pointsToRaster}}}
\examples{
-cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25), c(-50,0))
-cds2 <- rbind(c(80,20), c(140,60), c(160,0), c(140,-55), c(80,20))
-cds3 <- rbind(c(-180,-20), c(-140,55), c(-60,-20), c(-140,-60), c(-180,-20))
+
+cds1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
+cds2 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25), c(-50,0))
+cds3 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), Polygons(list(Polygon(cds2)), 2), Polygons(list(Polygon(cds3)), 3), Polygons(list(Polygon(hole)), 4) ))
polys at polygons[[4]]@Polygons[[1]]@hole <- TRUE
+
r <- raster()
-r <- polygonsToRaster(polys, r)
+r <- polygonsToRaster(polys, r, overlap='sum')
#plot(r)
#plot(polys, add=T)
-cds4 <- rbind(c(-180,10), c(0,40), c(40,45), c(145,-10), c(-25, -15), c(-180,0), c(-180,10))
+cds4 <- rbind(c(-180,10), c(0,90), c(40,90), c(145,-10), c(-25, -15), c(-180,0), c(-180,10))
addpoly <- SpatialPolygons(list(Polygons(list(Polygon(cds4)), 1)))
addpoly <- as(addpoly, "SpatialPolygonsDataFrame")
addpoly at data[1,1] <- 10
Modified: pkg/raster/man/focal.Rd
===================================================================
--- pkg/raster/man/focal.Rd 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/man/focal.Rd 2009-04-06 16:16:33 UTC (rev 398)
@@ -1,17 +1,15 @@
\name{focal}
\alias{focal}
-\alias{focalFilter}
\title{Focal}
\description{
-Calculate values for the neighborhood of focal cells
+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{
@@ -24,7 +22,6 @@
\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{
@@ -35,25 +32,15 @@
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.
}
+\seealso{ \code{\link[raster]{focalFilter}} }
+
+
\author{Robert J. Hijmans}
\examples{
Added: pkg/raster/man/focalFilter.Rd
===================================================================
--- pkg/raster/man/focalFilter.Rd (rev 0)
+++ pkg/raster/man/focalFilter.Rd 2009-04-06 16:16:33 UTC (rev 398)
@@ -0,0 +1,64 @@
+\name{focalFilter}
+
+\alias{focalFilter}
+
+\title{Focal filter}
+
+\description{
+Calculate values for the neighborhood of focal cells using a matrix of weights
+}
+
+\usage{
+focalFilter(raster, filter, fun=sum, filename="", overwrite=FALSE, filetype='raster', datatype='FLT4S', track=-1)
+}
+
+\arguments{
+ \item{raster}{A RasterLayer object}
+ \item{filter}{a matrix of weights. See Details}
+ \item{fun}{a function to apply to the product of the matrix of weights and the values}
+ \item{filename}{Output filename for a new raster}
+ \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{focalFilter} uses a matrix of weights for the neigbhorhood of the focal cells, together with a function (normally sum).
+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)}
+
+\code{focalFilter} can also be used to use non-rectangular areas. For example a combination of
+
+\code{filter=matrix(c(0,0,0,0,1,1,0,1,1), nrow=3)} and \code{fun=max} returns the max value for the lower-rigth corner of a 3x3 matrix
+around a focal cell
+
+}
+
+\value{
+A new RasterLayer object (in the R environment), and in some cases the side effect of a new file on disk.
+}
+
+\seealso{ \code{\link[raster]{focal}} }
+
+
+\author{Robert J. Hijmans}
+
+\examples{
+
+r <- raster(ncols=36, nrows=18)
+r[] <- runif(ncell(r))
+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}
+
Modified: pkg/raster/man/setDatatype.Rd
===================================================================
--- pkg/raster/man/setDatatype.Rd 2009-04-06 12:15:35 UTC (rev 397)
+++ pkg/raster/man/setDatatype.Rd 2009-04-06 16:16:33 UTC (rev 398)
@@ -28,18 +28,18 @@
The following datatypes are available.
\tabular{rll}{
-datatype \tab min \tab max \cr
-LOG1S \tab FALSE (0)\tab TRUE (1) \cr
-INT1S \tab -127 \tab 127 \cr
-INT1U \tab 0 \tab 255 \cr
-INT2S \tab -32,767\tab 32,767 \cr
-INT2U \tab 0 \tab 65,534 \cr
-INT4S \tab -2,147,483,647 \tab 2,147,483,647 \cr
-INT4U \tab 0 \tab 4,294,967,294 \cr
-INT8S \tab -9,223,372,036,854,775,807 \tab 9,223,372,036,854,775,807 \cr
-INT4U \tab 0 \tab 18,446,744,073,709,551,614 \cr
-FLT4S \tab -3.4E38 \tab 3.4E38 \cr
-FLT8S \tab -1.7E308 \tab 1.7E308 \cr
+\bold{datatype \tab min \tab max} \cr
+\code{LOG1S} \tab FALSE (0)\tab TRUE (1) \cr
+\code{INT1S} \tab -127 \tab 127 \cr
+\code{INT1U} \tab 0 \tab 255 \cr
+\code{INT2S} \tab -32,767\tab 32,767 \cr
+\code{INT2U} \tab 0 \tab 65,534 \cr
+\code{INT4S} \tab -2,147,483,647 \tab 2,147,483,647 \cr
+\code{INT4U} \tab 0 \tab 4,294,967,294 \cr
+\code{INT8S} \tab -9,223,372,036,854,775,807 \tab 9,223,372,036,854,775,807 \cr
+\code{INT4U} \tab 0 \tab 18,446,744,073,709,551,614 \cr
+\code{FLT4S} \tab -3.4E38 \tab 3.4E38 \cr
+\code{FLT8S} \tab -1.7E308 \tab 1.7E308 \cr
}
For all integer types, except the single byte types, the lowest (signed) or highest (unsigned) value is used to store NA.
More information about the Raster-commits
mailing list