[GRASS-dev] s.cellstats replacement for GRASS 6

[Helena Mitasova]

On May 27, 2006, at 2:09 AM, Hamish wrote:

> Hi,
>
>
> LIDAR & swath bathymetry people [3million+ input points]: I have an
> idea for a s.cellstats replacement. Will call it r.in.xyz (name is  
> from
> a 2001 wish by Roy Sanderson). Output map options will be the same as
> r.univar.

Hamish, s.cellstats inputs and outputs a site file, what you are  
proposing is supposed
to read an ascii x,y,z and output a raster file? If that is true, you  
may want to look
also at s.to.rast2  in GRASS Add-ons: http://grass.itc.it/outgoing/
(it includes an option to include data from surrounding cells, so if  
there are
small holes in data, they can be filled-in).
I have used s.to.rast a lot with lidar data or the binning option  
that provides
the on-line LDART tool (http://ekman.csc.noaa.gov/TCM/) and I miss it  
a lot in GRASS6.
Your suggestion to go directly from ascii points to raster
would get around the problems that the current vector support has  
with large data sets
(otherwise I would have suggested to do v.in.ascii and add this to  
v.to.rast so that attributes
can be used instead of z as well).

> The module will be quite memory intensive - for most stats binning  
> into
> a 10000x10000 region will take ~800mb RAM for a CELL,FCELL map & 1.6gb
> for a DCELL. Simple stats (min,max,n,sum) would take half that amount.
>
> A low memory option could cut complex stats memory requirements in  
> half
> (e.g. for parameter=mean do in two passes, first write "n" .tmp map to
> disk before populating "sum" in a second pass). Another idea to save
> memory is to cut the region up into segments and process the data in
> multiple passes. I blindly guess that the second method is more
> efficient.

I agree that processing it by segments is preferable.
>
> There'll be a test for z-range so the module can be used  
> iteratively to
> build 3D raster slices (for r.to.rast3). Same memory/segment stratagem
> could be employed to populate a 3D raster in multiple passes of a  
> single
> call to the module.
>
> The gist is to read in a line, test if the point is in the region
> bounds, then use column = int((input_x - west) / ew_res), row=same; to
> figure which cell to increment. Finally write out raster as f(a[,b]),
> row by row.
>
> Maybe the existing s.cellstats code is more efficient than my brute
> force way of keeping whole output map components in memory? The old
> sites format isn't too far from ascii x,y,z, so if s.cellstats is
> clearly a better method, maybe it is better to adapt that instead?

I have never used (or looked at the code of) s.cellstats or s.to.rast2,
so I cannot tell - let me know if it would be useful to try it out -
I have GRASS53 still running.
I have cc-d this to researchers working on the GEON project - we have  
discussed
the binning idea using quadtrees recently.
Andy may have some suggestions too,

Helena

> No idea.
>
>
> Hamish
>
> _______________________________________________
> grass-dev mailing list
> grass-dev at grass.itc.it
> http://grass.itc.it/mailman/listinfo/grass-dev