Re: write.asciigrid

Hi Roger and Andy,

Thanks for the help. I'm using ArcGIS 9.1. writeGDAL works as promised. The one odd behavior was the initial display in ArcGIS once I first added the two bands. Using the "stretched" classification, the map initially looked all black. Under properties, I changed the symbology to "Classified" and it looked fine -- then I changed back to "Stretched" and again it looked fine. Quirky behavior like this from ArcGIS is nothing new.

I tried the suggested write ascii command from maptools on my own data set and it didn't work initially due to the fact that my grid cells apparently are not perfect squares (perhaps this is a rounding issue). No need to keep fiddling with that, though, as writeGDAL does the trick. Thanks again,

Zev

Roger Bivand wrote:

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@xxxxxxxxxxx
www.zevross.com


    

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

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@xxxxxxxxxxxxxxxxx
https://stat.ethz.ch/mailman/listinfo/r-sig-geo