[Raster-commits] r202 - in pkg/raster: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 28 16:53:54 CET 2009
Author: rhijmans
Date: 2009-01-28 16:53:54 +0100 (Wed, 28 Jan 2009)
New Revision: 202
Added:
pkg/raster/R/polygonToRaster.R
pkg/raster/man/PolygonsToRaster.Rd
Modified:
pkg/raster/R/conversion.R
pkg/raster/R/get.R
Log:
Modified: pkg/raster/R/conversion.R
===================================================================
--- pkg/raster/R/conversion.R 2009-01-28 09:58:22 UTC (rev 201)
+++ pkg/raster/R/conversion.R 2009-01-28 15:53:54 UTC (rev 202)
@@ -64,11 +64,10 @@
dindex <- max(1, min(nlayers(object), index))
if (dindex != index) { warning(paste("index was changed to", dindex))}
rs <- object at layers[[dindex]]
-# rs <- newRaster(xmn = xmin(object), xmx = xmax(object), ymn = ymin(object), ymx = ymax(object), nrows=nrow(object), ncols=ncol(object), projstring=projection(object))
-# if (dataContent(object) == 'all') {
-# rs <- setValues(rs, as.matrix(values(object))[,dindex])
-# }
-# return(rs)
+ if (dataContent(object) == 'all') {
+ rs <- setValues(rs, values(object)[,dindex])
+ }
+ return(rs)
}
)
Modified: pkg/raster/R/get.R
===================================================================
--- pkg/raster/R/get.R 2009-01-28 09:58:22 UTC (rev 201)
+++ pkg/raster/R/get.R 2009-01-28 15:53:54 UTC (rev 202)
@@ -66,7 +66,7 @@
colnr <- colFromX(object, x[i]) - 1
rownr <- rowFromY(object, y[i]) - 1
if ((!is.na(colnr)) & (!is.na(rownr))) {
- cell[i] <- as.integer((rownr * ncol(object) + colnr) + 1)
+ cell[i] <- rownr * ncol(object) + colnr + 1
}
}
return(cell)
@@ -87,7 +87,7 @@
colnr <- (trunc((x - xmin(object)) / xres(object))) + 1
colnr[x == xmax(object)] <- ncol(object)
colnr[x < xmin(object) | x > xmax(object) ] <- NA
- return(colnr)
+ return(as.vector(colnr))
}
Added: pkg/raster/R/polygonToRaster.R
===================================================================
--- pkg/raster/R/polygonToRaster.R (rev 0)
+++ pkg/raster/R/polygonToRaster.R 2009-01-28 15:53:54 UTC (rev 202)
@@ -0,0 +1,194 @@
+# Author: Robert J. Hijmans, r.hijmans at gmail.com
+# International Rice Research Institute
+# Date : June 2008
+# Version 0,1
+# Licence GPL v3
+
+
+.intersectSegments <- function(x1, y1, x2, y2, x3, y3, x4, y4) {
+#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/
+ 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)))
+ }
+# 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)
+ }
+ 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)
+ outx <- x1 + (ua * (x2 - x1))
+ outy <- y1 + (ua * (y2 - y1))
+# print("C")
+ return(c(outx, outy))
+ } else {
+# print("D")
+ return(NA)
+ }
+}
+
+
+.overlayLinePolygon <- function(line, poly) {
+# for a simple line (connecting 2 points)
+ 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
+ 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) {
+ 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) {
+ resxy <- rbind(resxy, xy)
+ }
+ return(resxy)
+ }
+}
+
+
+
+polygonsToRaster <- function(sppoly, raster, field=0, filename="", overwrite=FALSE) {
+# check if bbox of raster and sppoly overlap
+ filename <- trim(filename)
+ raster <- setRaster(raster, filename)
+
+ spbb <- bbox(sppoly)
+ 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)
+ info <- matrix(NA, nrow=npol, ncol=3)
+ for (i in 1:npol) {
+ info[i,1] <- length(sppoly 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]))
+ }
+ 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(sppoly) == 'SpatialPolygons' | field == 0) {
+ putvals <- as.integer(1:npol)
+ } else {
+ putvals <- as.vector(sppoly at data[,field])
+ if (class(putvals) == 'character') {
+ stop('selected field is charater type')
+ }
+ }
+ raster <- setDatatype(raster, class(putvals[1]))
+
+
+ v <- vector(length=0)
+ for (r in 1:nrow(raster)) {
+ rv <- rep(NA, ncol(raster))
+ 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
+ } 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 ) {
+ # do nothing
+ } else {
+ mypoly <- sppoly at polygons[[i]]@Polygons[[j]]@coords
+ intersection <- .overlayLinePolygon(myline, mypoly)
+ if (nrow(intersection) > 0) {
+ cols <- sort(colFromX(raster, intersection[,1]))
+ for (k in 1:round(nrow(intersection)/2)) {
+ l <- (k * 2) - 1
+ rv[cols[l]:cols[l+1]] <- putvals[i]
+ }
+ }
+ }
+ }
+ }
+ }
+ if (filename == "") {
+ v <- c(v, rv)
+ } else {
+ raster <- setValues(raster, rv, r)
+ raster <- writeRaster(raster)
+ }
+ }
+ if (filename == "") {
+ raster <- setValues(raster, v)
+ }
+ return(raster)
+}
+
+
+polygonsToRaster2 <- function(sppoly, raster, field=0, filename="", overwrite=FALSE) {
+# should be slower except when polygons very detailed and raster las ow resolution
+
+# check if bbox of raster and sppoly overlap
+ filename <- trim(filename)
+ raster <- setRaster(raster, filename)
+
+ spbb <- bbox(sppoly)
+ 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))
+ } else {
+ putvals <- as.vector(sppoly at data[,field])
+ if (class(putvals) == 'character') {
+ stop('selected field is charater type')
+ }
+ }
+ raster <- setDatatype(raster, class(putvals[1]))
+
+
+ v <- vector(length=0)
+ rowcol <- cbind(0, 1:ncol(raster))
+
+ for (r in 1:nrow(raster)) {
+ rowcol[,1] <- r
+ sppoints <- xyFromCell(raster, cellFromRowCol(raster, rowcol[,1], rowcol[,2]), TRUE)
+ over <- overlay(sppoints, sppoly)
+ if (filename == "") {
+ v <- c(v, over)
+ } else {
+ raster <- setValues(raster, putvals[over], r)
+ raster <- writeRaster(raster)
+ }
+ }
+ if (filename == "") {
+ raster <- setValues(raster, v)
+ }
+ return(raster)
+}
+
Added: pkg/raster/man/PolygonsToRaster.Rd
===================================================================
--- pkg/raster/man/PolygonsToRaster.Rd (rev 0)
+++ pkg/raster/man/PolygonsToRaster.Rd 2009-01-28 15:53:54 UTC (rev 202)
@@ -0,0 +1,43 @@
+\name{polygons to raster}
+
+\alias{polygonsToRaster}
+
+\title{ Transform polygons to a raster }
+\description{
+ xxx
+}
+\usage{
+polygonsToRaster(sppoly, raster, field=0, filename="", overwrite=FALSE)
+}
+
+\arguments{
+ \item{sppoly}{ 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}{ }
+
+}
+
+\details{
+ For SpatialPolygons, the 'field' argument is ignored. The polygon index is tranferred.
+}
+
+\author{Robert J. Hijmans \email{r.hijmans at gmail.com}}
+
+\seealso{ \code{\link[rgdal]{rgdal}} }
+\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,-60), c(-60,-20), c(-140,55), c(-180,-20))
+
+polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), Polygons(list(Polygon(cds2)), 2), Polygons(list(Polygon(cds3)), 3) ))
+r <- newRaster()
+r <- polygonsToRaster(polys, r)
+#plot(r)
+#plot(polys, add=T)
+
+}
+
+\keyword{ spatial }
More information about the Raster-commits
mailing list