[R-sig-Geo] grass ascii support in the raster package

[Cornel Pop]

Dear all,

I just joined this list and I thought I would start off by making a small
contribution that may be of interest to some of you, namely a patch for the
raster package (replaces raster/R/rasterFromAscii.R) which:

1. Adds support for importing GRASS ASCII raster files.
2. Fixes what I think is an issue with the offset parameter in the import
function. Currently, if the offset is greater than the number of header
lines the function silently removes data without verifying that the number
of rows matches the nrows header parameter. This leads to potentially
difficult to diagnose errors (e.g. getValues(x) != ncell(x)). My modified
version determines the number of header lines and sets the correct offset
(i.e. excluding any data lines that may have been read in).
3. Generates an easy to understand error message if the header format is
not recognized.

In implementing these changes I tried to keep modifications of the original
code to a minimum (see the attached patch.txt). I'm sure things could be
improved, but I've used this patch for some time now and I haven't had any
issues.

Cheers,

Cornel M. Pop
---
PhD Candidate,
Max Planck Institute for Evolutionary Anthropology,
Leipzig, DE
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151118/e2b1e3e7/attachment.html>
-------------- next part --------------
--- rasterFromASCII_orig.R	2015-11-18 08:37:21.203485973 +0100
+++ final/rasterFromASCII.R	2015-11-18 08:47:39.073364056 +0100
@@ -3,12 +3,13 @@
 # Version 0.9
 # Licence GPL v3
 
+# Copyright 2015 (addendum) Cornel M. Pop <cornel at popgeology.com> Changed offset and added
+# option to load GRASS ASCII files. Offset parameter now ignored.
 
 .rasterFromASCIIFile <- function(filename, offset=6, crs=NA, ...) {
-	
-	offset <- as.integer(offset)
-	stopifnot(offset > 2)
-	
+	# Maximum real offset (GRASS = 6 required + 3 opts; ESRI = 5 required + 1 opt)
+  	offset <- 9
+  
 	splitasc <- function(s) {
 		s <- trim(s)
 		spl <- unlist(strsplit(s, ''), use.names = FALSE)
@@ -25,24 +26,49 @@
 	close(con)
 	ini <- lapply(lines, splitasc) 
 	ini <- matrix(unlist(ini, use.names = FALSE), ncol=2, byrow=TRUE)
-	
-	ini[,1] = toupper(ini[,1]) 
-	
-	w <- getOption('warn')
-	on.exit(options('warn' = w))
-	options('warn'=-1) 
-	test <- sum(as.numeric(ini[,1]), na.rm=TRUE) > 0
-	options('warn' = w)
-	if (test) {
-		m <- 'The header of this file appears to be incorrect: there are numbers where there should be keywords'
-		if (offset != 6) {
-			m <- paste(m, '\n  Are you using a wrong offset?', sep='')
-		} 
-		stop(m)
-	}
-	
-	
+
+	# Data may begin with a noval, which may not be numeric (i.e. NaN, NA)
+  	# Note that NODATA is not in the ESRI specs - using here for compatibility
+  	noval_kw <- c("NULL:", "NODATA", "NODATA_VALUE")
+
+  	noval <- ini[which(toupper(ini[,1]) %in% noval_kw == T), 2]
+  	if(length(noval) == 1){
+    	  # True offset. Rest is overshot (data)
+    	  offset <- length(which(is.na(suppressWarnings(try(as.numeric(ini[which(ini[,1] != noval),1]))))))
+  	} else {
+    	  offset <- length(which(is.na(suppressWarnings(try(as.numeric(ini[,1]))))))
+  	}
+  	stopifnot(offset > 2) # Note: This is not necessary.
+  	ini <- ini[1:offset,] # Keep only real header elements
+
+	ini[,1] = toupper(ini[,1])
+  
 	nodataval <- xn <- yn <- d <- nr <- nc <- xc <- yc <- NA
+  	grass_mandatory <- c("NORTH:", "SOUTH:", "EAST:", "WEST:", "ROWS:", "COLS:")
+  	if (!FALSE %in% (grass_mandatory %in% ini[,1])){
+   	  # GRASS ASCII - note: this is a bit different from how ESRI files are
+   	  # handled, in that we check for all mandatory header keys, and don't process
+    	  # the file if they are missing because, well, they should not be missing!
+    	  nr <- as.integer(ini[which(ini[,1] == "ROWS:"),2])
+    	  nc <- as.integer(ini[which(ini[,1] == "COLS:"),2])
+    	  xn <- as.numeric(ini[which(ini[,1] == "WEST:"),2])
+    	  yn <- as.numeric(ini[which(ini[,1] == "SOUTH:"),2])
+    	  xx <- as.numeric(ini[which(ini[,1] == "EAST:"),2])
+    	  yx <- as.numeric(ini[which(ini[,1] == "NORTH:"),2])
+    
+	  # Optional fields:
+    	  if ("NULL:" %in% ini[,1]){
+      	    nodataval <- as.numeric(ini[which(ini[,1] == "NULL:"),2])
+    	  } else {
+      	    warning('"NULL" not detected. Setting it to -Inf\n  You can set it to another value with function "NAvalue"')
+      	    nodataval <- -Inf
+    	  }
+    	  if ("TYPE:" %in% ini[,1]){
+     	    warning('Optional "TYPE" parameter detected, but ignoring and using the default data type.')
+    	  }
+    	  if ("MULTIPLIER:" %in% ini[,1]){multiplier <- as.numeric(ini[which(ini[,1] == "MULTIPLIER:"),2])}
+
+  	} else if ("NCOLS" %in% ini[,1]){
 	for (i in 1:nrow(ini)) {
 		if (ini[i,1] == "NCOLS") { nc <- as.integer(ini[i,2])
 		} else if (ini[i,1] == "NROWS") { nr <- as.integer(ini[i,2])
@@ -104,6 +130,9 @@
 	
 	xx <- xn + nc * d
 	yx <- yn + nr * d
+  	} else {
+    	  stop("File type not supported. Supported file types are: ESRI ASCII, GRASS ASCII")
+  	}
 
 	x <- raster(ncols=nc, nrows=nr, xmn=xn, ymn=yn, xmx=xx, ymx=yx, crs='')	
 
@@ -118,6 +147,9 @@
 		projection(x) <- crs
 	}
 	
+  	if(exists("multiplier")){
+		x <- setValues(x, getValues(x)*multiplier)
+  	}
 	return(x)
 }
 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rasterFromASCII.R
Type: text/x-r-source
Size: 5209 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20151118/e2b1e3e7/attachment.bin>