[Raster-commits] r217 - in pkg/raster: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 30 04:31:39 CET 2009
Author: rhijmans
Date: 2009-01-30 04:31:38 +0100 (Fri, 30 Jan 2009)
New Revision: 217
Modified:
pkg/raster/DESCRIPTION
pkg/raster/R/polygonToRaster.R
pkg/raster/R/read.raster.R
pkg/raster/R/values.R
pkg/raster/man/PolygonsToRaster.Rd
Log:
Modified: pkg/raster/DESCRIPTION
===================================================================
--- pkg/raster/DESCRIPTION 2009-01-29 15:11:08 UTC (rev 216)
+++ pkg/raster/DESCRIPTION 2009-01-30 03:31:38 UTC (rev 217)
@@ -1,11 +1,11 @@
Package: raster
Type: Package
Title: Raster data handling for geographic data analysis and modeling
-Version: 0.8.7-4
+Version: 0.8.7-5
Date: 30-Jan-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>
Description: Package for reading and manipulating (geographic, spatial) raster (grid) data
License: GPL3
-URL: http://www.r-project.org
\ No newline at end of file
+URL: http://raster.r-forge.r-project.org/
\ No newline at end of file
Modified: pkg/raster/R/polygonToRaster.R
===================================================================
--- pkg/raster/R/polygonToRaster.R 2009-01-29 15:11:08 UTC (rev 216)
+++ pkg/raster/R/polygonToRaster.R 2009-01-30 03:31:38 UTC (rev 217)
@@ -6,37 +6,35 @@
.intersectSegments <- function(x1, y1, x2, y2, x3, y3, x4, y4) {
-#From LISP code by Paul Reiners
+# Translated by RH from LISP code by Paul Reiners
# http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/linesegments.lisp
-# From algorithm by Paul Bourke given here: http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
+# Which was tranlated from the algorithm by Paul Bourke given here: http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
denom <- ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1))
ua_num <- ((x4 - x3) *(y1 - y3)) - ((y4 - y3) * (x1 - x3))
ub_num <- ((x2 - x1) *(y1 - y3)) - ((y2 - y1) * (x1 - x3))
# If the denominator and numerator for the equations for ua and ub are 0 then the two lines are coincident.
if ( denom == 0 & ua_num == 0 & ub_num == 0) {
-# print("A")
# return(c(x1, y1))
xmin <- max(x1, x3)
if (xmin==x1) {ymin <- y1} else {ymin <- y3}
xmax <- min(x2, x4)
if (xmax==x2) {ymax <- y2} else {ymax <- y4}
- return(rbind(c(xmin, ymin), c(xmax, ymax)))
+# RH: for coincident line (segments) returning two intersections : start and end
+ return(rbind(c(xmin, ymin),
+ c(xmax, ymax)))
}
# If the denominator for the equations for ua and ub is 0 then the two lines are parallel.
if (denom == 0) {
-# print("B")
- return(NA)
+ return(c(NA, NA))
}
ua <- ua_num / denom
ub <- ub_num / denom
if ((ua >= 0 & ua <= 1) & (ub >= 0 & ub <= 1) ) {
x <- x1 + ua * (x2 - x1)
y <- y1 + ua * (y2 - y1)
-# print("C")
return(c(x, y))
} else {
-# print("D")
- return(NA)
+ return(c(NA, NA))
}
}
@@ -52,12 +50,12 @@
for (i in 2:length(poly[,1])) {
#compute intersection
xy <- .intersectSegments(poly[i,1], poly[i,2], poly[i-1,1], poly[i-1,2], line[1,1], line[1,2], line[2,1], line[2,2] )
- if (length(xy) > 1) {
+ if (!is.na(xy[1])) {
resxy <- rbind(resxy, xy)
}
}
xy <- .intersectSegments(poly[1,1], poly[1,2], poly[length(poly[,1]),1], poly[length(poly[,1]),2], line[1,1], line[1,2], line[2,1], line[2,2] )
- if (length(xy) > 1) {
+ if (!is.na(xy[1])) {
resxy <- rbind(resxy, xy)
}
return(resxy)
@@ -66,9 +64,15 @@
-polygonsToRaster <- function(sppoly, raster, field=0, filename="", overwrite=FALSE) {
+polygonsToRaster <- function(sppoly, raster, field=0, filename="", overwrite=FALSE, updateRaster=FALSE, updateValue="NA") {
# check if bbox of raster and sppoly overlap
filename <- trim(filename)
+ if (updateRaster) {
+ oldraster <- raster
+ if (!(updateValue == 'NA' | updateValue == '!NA' | updateValue == 'all' | updateValue == 'zero')) {
+ stop('updateValue should be either "all", "NA", "!NA", or "zero"')
+ }
+ }
raster <- setRaster(raster, filename)
spbb <- bbox(sppoly)
@@ -146,6 +150,21 @@
}
}
}
+ if (updateRaster) {
+ oldvals <- values(readRow(oldraster, r))
+ if (updateValue == "all") {
+ ind <- which(!is.na(rv))
+ } else if (updateValue == "zero") {
+ ind <- which(oldvals==0 & !is.na(rv))
+ } else if (updateValue == "NA") {
+ ind <- which(is.na(oldvals))
+ } else {
+ ind <- which(!is.na(oldvals) & !is.na(rv))
+ }
+ oldvals[ind] <- rv[ind]
+ rv <- oldvals
+ }
+
if (filename == "") {
v <- c(v, rv)
} else {
Modified: pkg/raster/R/read.raster.R
===================================================================
--- pkg/raster/R/read.raster.R 2009-01-29 15:11:08 UTC (rev 216)
+++ pkg/raster/R/read.raster.R 2009-01-30 03:31:38 UTC (rev 217)
@@ -42,7 +42,7 @@
#read part of a single row
.rasterRead <- function(raster, rownr, startcol=1, ncolumns=(ncol(raster)-startcol+1)) {
rownr <- round(rownr)
- if (rownr == 0) { stop("rownr == 0. It should be between 1 and nrow(raster), or -1 for all rows") }
+ if (rownr == 0) { stop("rownr == 0. It should be between 1 and nrow(raster), or < 0 for all rows") }
if (rownr > nrow(raster)) { stop("rownr too high") }
if (startcol < 1) { stop("startcol < 1") }
if (startcol > ncol(raster)) { stop("startcol > ncol(raster)") }
@@ -53,8 +53,13 @@
endcol <- ncol(raster)
ncolumns <- ncol(raster) - startcol + 1
}
-
- if (.driver(raster) == 'raster') {
+
+ if (dataSource(raster)=='ram') {
+
+ result <- valuesRow(raster, rownr)[startcol:endcol]
+
+ } else if (.driver(raster) == 'raster') {
+
rastergri <- .setFileExtensionValues(filename(raster))
if (!file.exists( filename(raster))) {
stop(paste(filename(raster)," does not exist"))
@@ -83,9 +88,6 @@
else { #use GDAL
if (is.na(raster at file@band)) { result <- NA }
else {
- if (rownr > nrow(raster)) {
- stop("rownr too high")
- }
if (rownr <= 0) {
offs <- c(0, 0)
reg <- c(nrow(raster), ncol(raster)) # reg <- dim(raster at file@gdalhandle[[1]])
@@ -98,6 +100,7 @@
result <- getRasterData(raster at file@gdalhandle[[1]], offset=offs, region.dim=reg, band = raster at file@band)
if (!is.vector(result)) { result <- as.vector(result) }
}
+
raster at data@values <- as.vector(result)
if (rownr < 0) {
raster at data@indices <- c(1, ncell(raster))
Modified: pkg/raster/R/values.R
===================================================================
--- pkg/raster/R/values.R 2009-01-29 15:11:08 UTC (rev 216)
+++ pkg/raster/R/values.R 2009-01-30 03:31:38 UTC (rev 217)
@@ -27,14 +27,43 @@
valuesRow <- function(raster, rownr) {
- if (!(validRows(raster, rownr))) {stop(paste(rownr,'is not a valid rownumber')) }
- if (dataContent(raster) == 'sparse') {return (.valuesRow.sparse(raster, rownr))
- } else if (dataContent(raster) != 'all') {stop('cannot do. Need all data')
- } else {
+ if (dataContent(raster) == 'nodata') {
+ stop('no values in memory. First read or set values')
+ }
+ if (!(validRows(raster, rownr))) {
+ stop(paste(rownr,'is not a valid rownumber'))
+ }
+ if (dataContent(raster) == 'sparse') {
+ return (.valuesRow.sparse(raster, rownr))
+ } else if (dataContent(raster) == 'row') {
startcell <- cellFromRowCol(raster, rownr, 1)
endcell <- startcell+ncol(raster)-1
+ if (dataIndices(raster) == c(startcell, endcell)) {
+ return(values(raster))
+ } else {
+ stop('this row is not in memory. First use readRow() or readAll')
+ }
+ } else if (dataContent(raster) == 'block') {
+ firstcol <- colFromCell(raster, dataIndices(raster)[1])
+ lastcol <- colFromCell(raster, dataIndices(raster)[2])
+ if (firstcol != 1 | lastcol != ncol(raster)) {
+ stop('the block data in this raster does not have complete rows')
+ }
+ firstrow <- rowFromCell(raster, dataIndices(raster)[1])
+ lastrow <- rowFromCell(raster, dataIndices(raster)[2])
+ if (rownr < firstrow | rownr > lastrow) {
+ stop('this row is not in memory. First use readRow() or readAll')
+ }
+ startcell <- ((rownr - firstrow) * ncol(raster) + 1)
+ endcell <- startcell + ncol(raster) - 1
return(values(raster)[startcell:endcell])
- }
+ } else if (dataContent(raster) == 'all'){
+ startcell <- cellFromRowCol(raster, rownr, 1)
+ endcell <- startcell+ncol(raster)-1
+ return(values(raster)[startcell:endcell])
+ } else {
+ stop('something is wrong with the RasterLayer dataContent')
+ }
}
Modified: pkg/raster/man/PolygonsToRaster.Rd
===================================================================
--- pkg/raster/man/PolygonsToRaster.Rd 2009-01-29 15:11:08 UTC (rev 216)
+++ pkg/raster/man/PolygonsToRaster.Rd 2009-01-30 03:31:38 UTC (rev 217)
@@ -8,7 +8,7 @@
Two implementations of a polygon to raster conversion (number 2 should be slowest in most cases)
}
\usage{
-polygonsToRaster(sppoly, raster, field=0, filename="", overwrite=FALSE)
+polygonsToRaster(sppoly, raster, field=0, filename="", overwrite=FALSE, updateRaster=FALSE, updateValue="NA")
polygonsToRaster2(sppoly, raster, field=0, filename="", overwrite=FALSE)
}
@@ -18,7 +18,8 @@
\item{field}{ The index of the column in the SpatialPolygonsDataFrame to be transfered to the RasterLayer }
\item{filename}{ }
\item{overwrite}{ }
-
+ \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 'all', 'NA', '!NA', or 'zero' }
}
\details{
@@ -39,6 +40,15 @@
r <- polygonsToRaster(polys, r)
#plot(r)
#plot(polys, add=T)
+cds4 <- rbind(c(-180,20), c(0,90), c(40,85), c(150, 45), c(180,-5), c(145,-40), c(40,-35), c(0,-25), c(-25, -15), c(-180,20))
+addpoly <- SpatialPolygons(list(Polygons(list(Polygon(cds4)), 1)))
+addpoly <- as(addpoly, "SpatialPolygonsDataFrame")
+addpoly at data[1,1] <- 10
+r2 <- polygonsToRaster(addpoly, r, field=1, updateRaster=TRUE, updateValue="NA")
+#plot(r2)
+#plot(polys, add=T)
+#plot(addpoly, add=T, border="red", lwd=2)
+
}
\keyword{ spatial }
More information about the Raster-commits
mailing list