[R-sig-Geo] write.asciigrid

[Roger Bivand]

On Thu, 4 May 2006, Zev Ross wrote:

> Hi All,
> 
> I'd like to write the results of a kriging call to an ascii grid using
> GSTAT -- but I would like to write BOTH the predictions and the variances
> to grids. Is there a more elegant (and less dangerous) way to do it than
> what I have below?
> 

Dangerous for whom? Elegant, would be nice but life is short? The answer 
depends on the software that is going to read the output grids, and how it 
treats locales, etc., since both predictions and variances will be 
floating point.

> 
> depth_uk <- krige(DEPTH~slope.asc, depth, slope, vgm_depth_r)
> 
> #  write the predictions
> 
> write.asciigrid(depth_uk "c:/junk/rk/predictions.asc")
> 
> # replace predictions with variances and then write the variances
> 
> depth_uk$var1.pred<-depth_uk$var1.var
> write.asciigrid(depth_uk, "c:/junk/rk/.asc")
> 

If the software on the other side reads GeoTiff, my preference would be:

library(rgdal)
writeGDAL(depth_uk, "depth_uk.tif")

which I have seen work with ENVI, but not with ArcGIS 9.1; this preserves 
coordinate reference system metadata if set.

In a forthcoming release of rgdal, you should be able to pass options= to 
writeGDAL() - specifically INTERLEAVE=PIXEL, see:

http://grass.itc.it/grass61/manuals/html61_user/r.out.gdal.html

and I have seen this help with ArcGIS 9.1, although it wasn't predictable
(the legend scale showed correct values but the visualisation was wrong
sometimes - I tried on a Wednesday if that helps!). ArcGIS only accepted 
single band GeoTiff files, it thought 3-band were coloured images. ENVI 
simply read the GDAL-generated GeoTiffs (with 4 bands in the case we 
tried - point pattern kernel densities at different bandwidths) correctly 
without making any assumptions.

Depending on your locales, the ASCII grid route is being maintained in the 
maptools package and functions in the sp package will be deprecated. So

library(maptools)
writeAsciiGrid(depth_uk, "preds.txt", attr="var1.pred", dec=<your choice>)
writeAsciiGrid(depth_uk, "vars.txt", attr="var1.var", dec=<your choice>)

should get the values into ArcGIS 9.1 through the Toolbox (it is very 
sensitive to the "."/"," dec= setting). [The intention is to gather 
input/output functions in maptools and rgdal, freeing the other packages 
from having often older, duplicate copies of functions that do not get 
maintained.]

Again, how to do it does depend on what software is going to read the 
output ASCII grids, and what assumptions (often undocumented) it makes 
about the files.

Please let us know how you get on,

Roger

> PS this sample code comes from Tomislav Hengl's page
> (http://spatial-analyst.net/regkriging.php)
> 

Nice link!

> Zev
> --
> Zev Ross
> ZevRoss Spatial Analysis
> 303 Fairmount Ave
> Ithaca, NY 14850
> (607) 277-0004 (phone)
> (866) 877-3690 (fax toll-free)
> zev at zevross.com
> www.zevross.com
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no