[Raster-commits] r229 - in pkg/raster: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jan 31 15:04:35 CET 2009
Author: rhijmans
Date: 2009-01-31 15:04:35 +0100 (Sat, 31 Jan 2009)
New Revision: 229
Modified:
pkg/raster/R/polygonToRaster.R
pkg/raster/R/raster.create.R
pkg/raster/man/LinesToRaster.Rd
pkg/raster/man/PolygonsToRaster.Rd
pkg/raster/man/create.raster.Rd
Log:
Modified: pkg/raster/R/polygonToRaster.R
===================================================================
--- pkg/raster/R/polygonToRaster.R 2009-01-31 10:28:59 UTC (rev 228)
+++ pkg/raster/R/polygonToRaster.R 2009-01-31 14:04:35 UTC (rev 229)
@@ -40,22 +40,21 @@
.intersectLinePolygon <- function(line, poly) {
-# for a simple line (connecting 2 points) and a single poly
resxy <- matrix(NA, ncol=2, nrow=0)
- if (min(poly[,2]) > max(line[,2]) | max(poly[,2]) < min(line[,2])) {
- return(resxy)
- } else if (min(poly[,1]) > max(line[,1]) | max(poly[,1]) < min(line[,1])) {
- return(resxy)
- } else {
- for (i in 2:length(poly[,1])) {
-#compute intersection, the last one is not necessary as the the first point and last point coincide.
- 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 (!is.na(xy[1])) {
- resxy <- rbind(resxy, xy)
- }
+ miny <- min(line[,2])
+ maxy <- max(line[,2])
+ xyxy <- cbind(poly, rbind(poly[-1,], poly[1,]))
+ xyxy <- subset(xyxy, !( (xyxy[,2] > maxy & xyxy[,4] > maxy ) | (xyxy[,2] < miny & xyxy[,4] < miny)) )
+ if (length(xyxy) < 1) {
+ return(resxy)
+ }
+ for (i in 1:length(xyxy[,1])) {
+ xy <- .intersectSegments(xyxy[i,1], xyxy[i,2], xyxy[i,3], xyxy[i,4], line[1,1], line[1,2], line[2,1], line[2,2] )
+ if (!is.na(xy[1])) {
+ resxy <- rbind(resxy, xy)
}
- return(resxy)
}
+ return(resxy)
}
@@ -105,7 +104,7 @@
}
raster <- setDatatype(raster, class(putvals[1]))
- adj <- 0.49 * xres(raster)
+ adj <- 0.5 * xres(raster)
v <- vector(length=0)
rxmn <- xmin(raster) + 0.1 * xres(raster)
rxmx <- xmax(raster) - 0.1 * xres(raster)
@@ -116,19 +115,18 @@
ly <- yFromRow(raster, r)
myline <- rbind(c(lxmin,ly), c(lxmax,ly))
- uly <- ly + 0.01 * yres(raster)
- lly <- ly - 0.01 * yres(raster)
for (i in 1:npol) {
- if (info[i,2] > uly | info[i,3] < lly) {
- # do nothing
+
+ if (info[i,2] > ly | info[i,3] < ly) {
+ # entire polygon above or below row. do nothing
} else {
for (j in 1:info[i,1]) {
if ( max ( spPolys at polygons[[i]]@Polygons[[j]]@coords[,2] ) < ly | min( spPolys at polygons[[i]]@Polygons[[j]]@coords[,2] ) > ly ) {
- # do nothing
+ # polygon part above or below row. do nothing
} else {
mypoly <- spPolys at polygons[[i]]@Polygons[[j]]@coords
+ intersection <- .intersectLinePolygon(myline, mypoly)
- intersection <- .intersectLinePolygon(myline, mypoly)
if (nrow(intersection) > 0) {
x <- sort(intersection[,1])
for (k in 1:round(nrow(intersection)/2)) {
@@ -140,11 +138,11 @@
# adjust to skip first cell if the center is not covered by this polygon
x1a <- x1 + adj
x2a <- x2 - adj
- x1a <- max(rxmn, x1a)
- x2a <- min(rxmx, x2a)
+ x1a <- min(rxmx, max(rxmn, x1a))
+ x2a <- min(rxmx, max(rxmn, x2a))
col1 <- colFromX(raster, x1a)
col2 <- colFromX(raster, x2a)
- if (is.na(col1) | is.na(col2) | col1 > col2) { next }
+ if (col1 > col2) { next }
if ( spPolys at polygons[[i]]@Polygons[[j]]@hole ) {
holes[col1:col2] <- TRUE
@@ -192,6 +190,8 @@
# This is based on sampling by points. Should be slower except when polygons very detailed and raster las ow resolution
# but it could be optimized further
+# this version does not deal with polygon holes
+
# check if bbox of raster and spPolys overlap
filename <- trim(filename)
raster <- setRaster(raster, filename)
Modified: pkg/raster/R/raster.create.R
===================================================================
--- pkg/raster/R/raster.create.R 2009-01-31 10:28:59 UTC (rev 228)
+++ pkg/raster/R/raster.create.R 2009-01-31 14:04:35 UTC (rev 229)
@@ -30,23 +30,28 @@
return(rs)
}
-
-
#if (!isGeneric("values")) {
# setGeneric("values", function(object, ...)
# standardGeneric("values"))
#}
#setMethod('values', signature(object='Raster'),
-rasterFromBbox <- function(bndbox, nrows=1, ncols=1) {
+
+rasterFromBbox <- function(bndbox, nrows=10, ncols=10, nudge=TRUE) {
crs <- newCRS('NA')
try(crs <- projection(bndbox, asText=F), silent = T)
- bndbox <- getBbox(bndbox)
+ bb <- getBbox(bndbox)
+ if (nudge) {
+ bb at xmin <- floor(bb at xmin)
+ bb at ymin <- floor(bb at ymin)
+ bb at xmax <- ceiling(bb at xmax)
+ bb at ymax <- ceiling(bb at ymax)
+ }
nr = as.integer(round(nrows))
nc = as.integer(round(ncols))
if (nc < 1) { stop("ncols should be > 0") }
if (nr < 1) { stop("nrows should be > 0") }
- raster <- new("RasterLayer", bbox = bndbox, crs=crs, ncols = nc, nrows = nr )
+ raster <- new("RasterLayer", bbox = bb, crs=crs, ncols = nc, nrows = nr )
return(raster)
}
Modified: pkg/raster/man/LinesToRaster.Rd
===================================================================
--- pkg/raster/man/LinesToRaster.Rd 2009-01-31 10:28:59 UTC (rev 228)
+++ pkg/raster/man/LinesToRaster.Rd 2009-01-31 14:04:35 UTC (rev 229)
@@ -13,7 +13,7 @@
\arguments{
\item{spLines}{ a SpatialLines or a SpatialLinesDataFrame object (sp package)}
\item{raster}{ a RasterLayer object}
- \item{field}{ The index of the column in the SpatialPolygonsDataFrame to be transfered to the RasterLayer }
+ \item{field}{ The index of the column in the SpatialLinesDataFrame 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{updateRaster}{logical. If \code{TRUE} the values of the input RasterLayer are updated where the polygons overlap cells }
Modified: pkg/raster/man/PolygonsToRaster.Rd
===================================================================
--- pkg/raster/man/PolygonsToRaster.Rd 2009-01-31 10:28:59 UTC (rev 228)
+++ pkg/raster/man/PolygonsToRaster.Rd 2009-01-31 14:04:35 UTC (rev 229)
@@ -16,8 +16,8 @@
\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}{ }
- \item{overwrite}{ }
+ \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 'all', 'NA', '!NA', or 'zero' }
}
@@ -28,7 +28,7 @@
\author{Robert J. Hijmans \email{r.hijmans at gmail.com}}
-\seealso{ \code{\link[rgdal]{rgdal}} }
+\seealso{ \code{\link[raster]{lines to raster}} }
\examples{
cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25), c(-50,0))
@@ -48,8 +48,8 @@
addpoly at data[1,1] <- 10
r2 <- polygonsToRaster(addpoly, r, field=1, updateRaster=TRUE, updateValue="NA")
plot(r2)
-plot(polys, border="blue", lwd=2, add=TRUE)
-plot(addpoly, add=TRUE, border="red", lwd=2)
+#plot(polys, border="blue", lwd=2, add=TRUE)
+#plot(addpoly, add=TRUE, border="red", lwd=2)
}
Modified: pkg/raster/man/create.raster.Rd
===================================================================
--- pkg/raster/man/create.raster.Rd 2009-01-31 10:28:59 UTC (rev 228)
+++ pkg/raster/man/create.raster.Rd 2009-01-31 14:04:35 UTC (rev 229)
@@ -16,7 +16,7 @@
\usage{
raster(xmn = -180, xmx = 180, ymn = -90, ymx = 90, nrows = 180, ncols = 360, projstring = "+proj=longlat +datum=WGS84")
newRaster(xmn = -180, xmx = 180, ymn = -90, ymx = 90, nrows = 180, ncols = 360, projstring = "+proj=longlat +datum=WGS84")
- rasterFromBbox(bndbox, nrows = 1, ncols = 1)
+ rasterFromBbox(bndbox, nrows = 10, ncols = 10, nudge=TRUE)
rasterFromFile(filename, values=FALSE, band=1)
closeHandle(raster)
}
@@ -30,6 +30,7 @@
\item{ymx}{ maximum y coordinate of raster (top border) }
\item{nrows}{ number of rows on raster }
\item{ncols}{ number of columns on raster }
+ \item{nudge}{logical. If \code{TRUE}, the boundingbox is nudged outwards (to the next integer) }
\item{projstring}{ PROJ4 type description of a map projection}
\item{bndbox}{bounding box used to crop a raster. Any object that is or has a bounding box (see below) }
\item{values}{logical. If \code{TRUE}, RasterLayer values will be read into memory with 'readAll()'}
More information about the Raster-commits
mailing list