[Raster-commits] r461 - pkg/raster/R

[noreply at r-forge.r-project.org]

Author: rhijmans
Date: 2009-05-19 04:29:44 +0200 (Tue, 19 May 2009)
New Revision: 461

Modified:
   pkg/raster/R/projectRaster.R
   pkg/raster/R/writeAscii.R
Log:


Modified: pkg/raster/R/projectRaster.R
===================================================================
--- pkg/raster/R/projectRaster.R	2009-05-19 00:15:19 UTC (rev 460)
+++ pkg/raster/R/projectRaster.R	2009-05-19 02:29:44 UTC (rev 461)
@@ -10,19 +10,49 @@
 	validObject(projection(projs, asText=FALSE))
 	projfrom <- projection(object)
 	projto <- projection(projs)
-	b <- extent(object)
-	xy <- rbind(c(b at xmin, b at ymax), c(b at xmax, b at ymax), c(b at xmin, b at ymin), c(b at xmax, b at ymin))
+	
+	xy1 <- xyFromCell(object, cellFromCol(object, 1))
+	xy1[,1] <- xy1[,1] - 0.5 * xres(object)
+	xy1[1,2] <- xy1[1,2] + 0.5 * yres(object)
+	xy1[nrow(object),2] <- xy1[nrow(object),2] + 0.5 * yres(object)
+	
+	xy2 <- xyFromCell(object, cellFromCol(object, ncol(object)))
+	xy2[,1] <- xy2[,1] + 0.5 * xres(object)
+	xy2[1,2] <- xy2[1,2] + 0.5 * yres(object)
+	xy2[nrow(object),2] <- xy2[nrow(object),2] + 0.5 * yres(object)
+	
+	xy <- rbind(xy1, xy2)
+	
 	res <- .Call("transform", projfrom, projto, nrow(xy), xy[,1], xy[,2], PACKAGE="rgdal")
-	p <- cbind(res[[1]], res[[2]])
-	if (any(is.infinite(p[,1])) || any(is.infinite(p[,2]))) {
-		stop("non finite transformation detected. This probably means that the map projection chosen is not appropriate for this geographic area")
+	
+	x <- res[[1]]
+	y <- res[[2]]
+	xy <- cbind(x, y)
+	xy <- subset(xy, !(is.infinite(xy[,1]) | is.infinite(xy[,2])) )
+	x <- xy[,1]
+	y <- xy[,2]
+	
+	if (length(y) == 0 | length(y) ==0) { stop("cannot do this transformation") }
+	minx <- min(x)
+	maxx <- max(x)
+	if (maxx == minx) {
+		maxx <- maxx + 0.5
+		minx <- minx - 0.5
 	}
-	obj <- raster(newExtent(min(p[,1]), max(p[,1]), min(p[,2]),  max(p[,2])), nrows=nrow(object), ncols=ncol(object), proj=(projs))
+	miny <- min(y)
+	maxy <- max(y)
+	if (maxy == miny) {
+		maxy <- maxy + 0.5
+		miny <- miny - 0.5
+	}
+	
+	obj <- raster(newExtent(minx, maxx, miny,  maxy), nrows=nrow(object), ncols=ncol(object), proj=(projs))
 	return(obj)
 }
 
 
 projectRaster <- function(from, to, method="ngb", filename="", filetype='raster', datatype='FLT4S', overwrite=FALSE, track=-1)  {
+	
 	if (dataContent(from) != 'all' & dataSource(from) == 'ram') { stop('no vales for "from". Nothing to do') }
 
 	validObject(to)

Modified: pkg/raster/R/writeAscii.R
===================================================================
--- pkg/raster/R/writeAscii.R	2009-05-19 00:15:19 UTC (rev 460)
+++ pkg/raster/R/writeAscii.R	2009-05-19 02:29:44 UTC (rev 461)
@@ -14,8 +14,10 @@
 	
 	if (dataIndices(raster)[1] == 1) {
 		resdif <- abs((yres(raster) - xres(raster)) / yres(raster) )
-		if (resdif > 0.01) {
+		if (resdif > 0.001) {
 			stop(paste("raster has unequal horizontal and vertical resolutions","\n", "these data cannot be stored in arc-ascii format"))
+		} else if (resdif > 0.00001) {
+			warning("arc-ascii format ignore that this raster has slightly unequal horizontal and vertical resolutions")
 		}
 		if (!overwrite & file.exists(filename)) {
 				stop(paste(filename, "exists. Use 'overwrite=TRUE'"))