[Raster-commits] r227 - pkg/raster/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jan 31 10:13:37 CET 2009
Author: rhijmans
Date: 2009-01-31 10:13:36 +0100 (Sat, 31 Jan 2009)
New Revision: 227
Added:
pkg/raster/R/linesToRaster.R
Modified:
pkg/raster/R/conversion.R
pkg/raster/R/polygonToRaster.R
pkg/raster/R/set.R
Log:
Modified: pkg/raster/R/conversion.R
===================================================================
--- pkg/raster/R/conversion.R 2009-01-31 08:25:31 UTC (rev 226)
+++ pkg/raster/R/conversion.R 2009-01-31 09:13:36 UTC (rev 227)
@@ -97,11 +97,11 @@
setMethod('asRasterLayer', signature(object='SpatialPixels', index='numeric'),
function(object, index){
- raster <- raster()
- raster <- setBbox(raster, getBbox(object))
- raster <- setProjection(raster, object at proj4string)
- raster <- setRowCol(raster, object at grid@cells.dim[2], object at grid@cells.dim[1])
- return(raster)
+ r <- raster()
+ r <- setBbox(r, getBbox(object))
+ r <- setProjection(r, object at proj4string)
+ r <- setRowCol(r, object at grid@cells.dim[2], object at grid@cells.dim[1])
+ return(r)
}
)
Added: pkg/raster/R/linesToRaster.R
===================================================================
--- pkg/raster/R/linesToRaster.R (rev 0)
+++ pkg/raster/R/linesToRaster.R 2009-01-31 09:13:36 UTC (rev 227)
@@ -0,0 +1,169 @@
+# Author: Robert J. Hijmans, r.hijmans at gmail.com
+# International Rice Research Institute
+# Date : June 2008
+# Version 0,1
+# Licence GPL v3
+
+.specialRowFromY <- function(object, y) {
+ rownr <- 1 + (trunc((ymax(object) - y)/yres(object)))
+ rownr[y == ymin(object)] <- nrow(object)
+ rownr[y > ymax(object)] <- -1
+ rownr[y < ymin(object)] <- nrow(object) + 1
+ return(rownr)
+}
+
+.specialColFromX <- function(object, x) {
+ colnr <- (trunc((x - xmin(object))/xres(object))) + 1
+ colnr[x == xmax(object)] <- ncol(object)
+ colnr[x < xmin(object)] <- -1
+ colnr[x > xmax(object)] <- ncol(object) + 1
+ return(colnr)
+}
+
+
+getCols <- function(rs, rownr, segment, line1, line2) {
+# for a simple line (connecting 2 points) and a single poly
+ rows <- .specialRowFromY(rs, segment[,2])
+ if (rows[1] > rownr & rows[2] > rownr | rows[1] < rownr & rows[2] < rownr) { return(NA) }
+ cols <- .specialColFromX(rs,segment[,1])
+ rowcol <- cbind(rows, cols)[order(cols),]
+ if (rows[1] == rows[2]) {
+ return(rowcol[1,2]:rowcol[2,2])
+ } else {
+ if (rowcol[1,1] == rownr) {
+ col1 <- rowcol[1,2]
+ if (rowcol[2,1] < rownr) {
+ xy <- .intersectSegments(line1[1,1], line1[1,2], line1[2,1], line1[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ } else {
+ xy <- .intersectSegments(line2[1,1], line2[1,2], line2[2,1], line2[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ }
+ xy <- t(as.matrix(xy))
+ col2 <- max(colFromX(rs, xy[,1]))
+ } else if (rowcol[2,1] == rownr) {
+ col1 <- rowcol[2,2]
+ if (rowcol[1,1] < rownr) {
+ xy <- .intersectSegments(line1[1,1], line1[1,2], line1[2,1], line1[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ } else {
+ xy <- .intersectSegments(line2[1,1], line2[1,2], line2[2,1], line2[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ }
+ xy <- t(as.matrix(xy))
+ col2 <- max(colFromX(rs, xy[,1]))
+ } else {
+ xy1 <- .intersectSegments(line1[1,1], line1[1,2], line1[2,1], line1[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ xy2 <- .intersectSegments(line2[1,1], line2[1,2], line2[2,1], line2[2,2], segment[1,1], segment[1,2], segment[2,1], segment[2,2] )
+ xy <- rbind(xy1, xy2)
+ cols <-colFromX(rs, xy[,1])
+ col1 <- min(cols)
+ col2 <- max(cols)
+ }
+ return(col1:col2)
+ }
+}
+
+
+
+
+
+linesToRaster <- function(spLines, raster, field=0, filename="", overwrite=FALSE, updateRaster=FALSE, updateValue="NA") {
+ 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)
+
+# check if bbox of raster and spLines overlap
+ spbb <- bbox(spLines)
+ rsbb <- bbox(raster)
+ if (spbb[1,1] > rsbb[1,2] | spbb[2,1] > rsbb[2,2]) {
+ stop('lines and raster have no overlapping areas')
+ }
+ nline <- length(spLines at lines)
+ info <- matrix(NA, nrow=nline, ncol=3)
+ for (i in 1:nline) {
+# holes <- sapply(rings, function(y) slot(y, "hole"))
+# areas <- sapply(rings, function(x) slot(x, "area"))
+
+ info[i,1] <- length(spLines at lines[[i]]@Lines)
+ miny <- NULL
+ maxy <- NULL
+ for (j in 1:info[i,1]) {
+ miny <- min(miny, min(spLines at lines[[i]]@Lines[[j]]@coords[,2]))
+ maxy <- max(maxy, max(spLines at lines[[i]]@Lines[[j]]@coords[,2]))
+ }
+ info[i,2] <- miny
+ info[i,3] <- maxy
+ }
+ lxmin <- min(spbb[1,1], rsbb[1,1]) - xres(raster)
+ lxmax <- max(spbb[1,2], rsbb[1,2]) + xres(raster)
+
+ if (class(spLines) == 'SpatialPolygons' | field == 0) {
+ putvals <- as.integer(1:nline)
+ } else {
+ putvals <- as.vector(spLines at data[,field])
+ if (class(putvals) == 'character') {
+ stop('selected field is charater type')
+ }
+ }
+ raster <- setDatatype(raster, class(putvals[1]))
+
+ v <- vector(length=0)
+ rxmn <- xmin(raster) + 0.1 * xres(raster)
+ rxmx <- xmax(raster) - 0.1 * xres(raster)
+ for (r in 1:nrow(raster)) {
+ rv <- rep(NA, ncol(raster))
+
+ ly <- yFromRow(raster, r)
+ line1 <- rbind(c(lxmin, ly + 0.5*yres(raster)), c(lxmax,ly + 0.5*yres(raster)))
+ line2 <- rbind(c(lxmin, ly - 0.5*yres(raster)), c(lxmax,ly - 0.5*yres(raster)))
+
+ uly <- ly + 0.01 * yres(raster)
+ lly <- ly - 0.01 * yres(raster)
+ for (i in 1:nline) {
+ if (info[i,2] > uly | info[i,3] < lly) {
+ # do nothing
+ } else {
+ for (j in 1:info[i,1]) {
+ if ( max ( spLines at lines[[i]]@Lines[[j]]@coords[,2] ) < ly | min( spLines at lines[[i]]@Lines[[j]]@coords[,2] ) > ly ) {
+ # do nothing
+ } else {
+ segment <- spLines at lines[[i]]@Lines[[j]]@coords
+ colnrs <- getCols(raster, r, segment, line1, line2)
+ if ( length(colnrs) > 0 ) {
+ rv[colnrs] <- putvals[i]
+ }
+ }
+ }
+ }
+ }
+
+ 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 {
+ raster <- setValues(raster, values=rv, rownr=r)
+ raster <- writeRaster(raster)
+ }
+ }
+ if (filename == "") {
+ raster <- setValues(raster, v)
+ }
+ return(raster)
+}
+
Modified: pkg/raster/R/polygonToRaster.R
===================================================================
--- pkg/raster/R/polygonToRaster.R 2009-01-31 08:25:31 UTC (rev 226)
+++ pkg/raster/R/polygonToRaster.R 2009-01-31 09:13:36 UTC (rev 227)
@@ -39,7 +39,7 @@
}
-.overlayLinePolygon <- function(line, poly) {
+.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])) {
@@ -60,7 +60,7 @@
-polygonsToRaster <- function(sppoly, raster, field=0, filename="", overwrite=FALSE, updateRaster=FALSE, updateValue="NA") {
+polygonsToRaster <- function(spPolys, raster, field=0, filename="", overwrite=FALSE, updateRaster=FALSE, updateValue="NA") {
filename <- trim(filename)
if (updateRaster) {
oldraster <- raster
@@ -70,24 +70,24 @@
}
raster <- setRaster(raster, filename)
-# check if bbox of raster and sppoly overlap
- spbb <- bbox(sppoly)
+# check if bbox of raster and spPolys overlap
+ spbb <- bbox(spPolys)
rsbb <- bbox(raster)
if (spbb[1,1] > rsbb[1,2] | spbb[2,1] > rsbb[2,2]) {
stop('polygon and raster have no overlapping areas')
}
- npol <- length(sppoly at polygons)
+ npol <- length(spPolys at polygons)
info <- matrix(NA, nrow=npol, ncol=3)
for (i in 1:npol) {
# holes <- sapply(rings, function(y) slot(y, "hole"))
# areas <- sapply(rings, function(x) slot(x, "area"))
- info[i,1] <- length(sppoly at polygons[[i]]@Polygons)
+ info[i,1] <- length(spPolys at polygons[[i]]@Polygons)
miny <- NULL
maxy <- NULL
for (j in 1:info[i,1]) {
- miny <- min(miny, min(sppoly at polygons[[i]]@Polygons[[j]]@coords[,2]))
- maxy <- max(maxy, max(sppoly at polygons[[i]]@Polygons[[j]]@coords[,2]))
+ miny <- min(miny, min(spPolys at polygons[[i]]@Polygons[[j]]@coords[,2]))
+ maxy <- max(maxy, max(spPolys at polygons[[i]]@Polygons[[j]]@coords[,2]))
}
info[i,2] <- miny
info[i,3] <- maxy
@@ -95,10 +95,10 @@
lxmin <- min(spbb[1,1], rsbb[1,1]) - xres(raster)
lxmax <- max(spbb[1,2], rsbb[1,2]) + xres(raster)
- if (class(sppoly) == 'SpatialPolygons' | field == 0) {
+ if (class(spPolys) == 'SpatialPolygons' | field == 0) {
putvals <- as.integer(1:npol)
} else {
- putvals <- as.vector(sppoly at data[,field])
+ putvals <- as.vector(spPolys at data[,field])
if (class(putvals) == 'character') {
stop('selected field is charater type')
}
@@ -123,12 +123,12 @@
# do nothing
} else {
for (j in 1:info[i,1]) {
- if ( max ( sppoly at polygons[[i]]@Polygons[[j]]@coords[,2] ) < ly | min( sppoly at polygons[[i]]@Polygons[[j]]@coords[,2] ) > ly ) {
+ if ( max ( spPolys at polygons[[i]]@Polygons[[j]]@coords[,2] ) < ly | min( spPolys at polygons[[i]]@Polygons[[j]]@coords[,2] ) > ly ) {
# do nothing
} else {
- mypoly <- sppoly at polygons[[i]]@Polygons[[j]]@coords
+ mypoly <- spPolys at polygons[[i]]@Polygons[[j]]@coords
- intersection <- .overlayLinePolygon(myline, mypoly)
+ intersection <- .intersectLinePolygon(myline, mypoly)
if (nrow(intersection) > 0) {
x <- sort(intersection[,1])
for (k in 1:round(nrow(intersection)/2)) {
@@ -146,7 +146,7 @@
col2 <- colFromX(raster, x2a)
if (is.na(col1) | is.na(col2) | col1 > col2) { next }
- if ( sppoly at polygons[[i]]@Polygons[[j]]@hole ) {
+ if ( spPolys at polygons[[i]]@Polygons[[j]]@hole ) {
holes[col1:col2] <- TRUE
} else {
rv[col1:col2] <- putvals[i]
@@ -188,24 +188,24 @@
}
-polygonsToRaster2 <- function(sppoly, raster, field=0, filename="", overwrite=FALSE) {
+polygonsToRaster2 <- function(spPolys, raster, field=0, filename="", overwrite=FALSE) {
# 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
-# check if bbox of raster and sppoly overlap
+# check if bbox of raster and spPolys overlap
filename <- trim(filename)
raster <- setRaster(raster, filename)
- spbb <- bbox(sppoly)
+ spbb <- bbox(spPolys)
rsbb <- bbox(raster)
if (spbb[1,1] > rsbb[1,2] | spbb[2,1] > rsbb[2,2]) {
stop('polygon and raster have no overlapping areas')
}
- if (class(sppoly) == 'SpatialPolygons' | field == 0) {
- putvals <- as.integer(1:length(sppoly at polygons))
+ if (class(spPolys) == 'SpatialPolygons' | field == 0) {
+ putvals <- as.integer(1:length(spPolys at polygons))
} else {
- putvals <- as.vector(sppoly at data[,field])
+ putvals <- as.vector(spPolys at data[,field])
if (class(putvals) == 'character') {
stop('selected field is charater type')
}
@@ -225,7 +225,7 @@
} else {
rowcol[,1] <- r
sppoints <- xyFromCell(raster, cellFromRowCol(raster, rowcol[,1], rowcol[,2]), TRUE)
- over <- overlay(sppoints, sppoly)
+ over <- overlay(sppoints, spPolys)
vals <- putvals[over]
}
if (filename == "") {
Modified: pkg/raster/R/set.R
===================================================================
--- pkg/raster/R/set.R 2009-01-31 08:25:31 UTC (rev 226)
+++ pkg/raster/R/set.R 2009-01-31 09:13:36 UTC (rev 227)
@@ -41,13 +41,13 @@
stop("it is not allowed to set the filename of the output RasterLayer to that of the input RasterLayer")
}
- raster <- raster(xmn = xmin(object), xmx = xmax(object), ymn = ymin(object), ymx = ymax(object), nrows=nrow(object), ncols=ncol(object), projstring=projection(object))
- raster <- setFilename(raster, filename)
+ r <- raster(xmn = xmin(object), xmx = xmax(object), ymn = ymin(object), ymx = ymax(object), nrows=nrow(object), ncols=ncol(object), projstring=projection(object))
+ r <- setFilename(r, filename)
- if ( length(values) != 1 | ( length(values) == 1 & ncell(raster) == 1) ) {
- raster <- setValues(raster, values)
+ if ( length(values) != 1 | ( length(values) == 1 & ncell(r) == 1) ) {
+ r <- setValues(r, values)
}
- return(raster)
+ return(r)
}
setFilename <- function(object, filename) {
More information about the Raster-commits
mailing list