[R-sig-Geo] write.asciigrid

[Roger Bivand]

On Wed, 10 May 2006, Zev Ross wrote:

> Hi Roger,
> 
> No, I didn't use the calculate statistics trick. I just opened the
> rasters and when they didn't display properly, I simply changed to
> "Stretched" and then back to "Classified" and it looked fine. I should
> mention though, that I added each band separately rather than adding both
> bands at once. Zev

Me too - it seems that on occasion with more than a single band, Arc 
things it is a 3-band image and tries to apply RGB with curious results. 
ENVI reads writeGDAL() GeoTiffs, including multiple FLOAT bands, right 
first time. On single bands, the "Compute Statistics" trick helped for me.

Roger

> 
> Roger Bivand wrote:
> 
>  On Tue, 9 May 2006, Zev Ross wrote:
> 
>   
> 
>  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.
>     
> 
>  Did you do the:
> 
> ArcToolbox | Data Management Tools | Raster | Calculate Statistics
> 
> trick too? I think that is what Arc (sometimes) does on the fly when 
> changing symbology (dreadful term!). But when NAs are present, it can get 
> stuck.
> 
>   
> 
>  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,
> 
>     
> 
>  Yes, they need to be equal for the Arc ASCII grid route to work. If they 
> are, it works (when the decimal point in FLOAT is correct too).
> 
> Roger
> 
>   
> 
>  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 at zevross.com
> 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 at zevross.com
> 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 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