[GRASS-dev] CELL/FCELL/DCELL [was: Re: r.mapcal rand() strangeness]

[Glynn Clements]

Maciej Sieczka wrote:

> >> I was going to add this information, but I'm not sure if I understand 
> >> correctly that GRID3D is always floating point. Raster intro seems to 
> >> suggest so, but r3.mapcalc suggests the opposite, eg.: "Note: If you 
> >> calculate with integer numbers, the resulting map will be integer". 
> >> Please tell me.
> 
> > I suspect that's probably an artifact of verbatim copying from the
> > r.mapcalc documentation.
> > 
> > I'm fairly sure that G3D only supports FP.
> 
> Can you confirm if double precision only?

G3D has both single and double precision formats.

> I guess so by experiments with 
> r3.mapcalc - eg.:
> 
> $ r3.mapcalc 'map=float(1)'
> $ r3.info map -t datatype="double"

Right. r3.mapcalc sets the in-memory representation appropriately, but
it doesn't set the file representation. AFAICT, it requires the
following change:

--- raster/r.mapcalc/map3.c	(revision 30389)
+++ raster/r.mapcalc/map3.c	(working copy)
@@ -588,6 +588,8 @@
 	void *handle;
 	int fd;
 
+	G3d_setFileType(res_type == FCELL_TYPE ? FCELL_TYPE : DCELL_TYPE);
+
 	handle = G3d_openCellNew(
 		(char *) name,
 		res_type == FCELL_TYPE ? FCELL_TYPE : DCELL_TYPE,

> > Looking at the r3.mapcalc code, the only way that you can get CELL
>  > values from a 3D grid is if you use the various # operators to perform
>  > colour lookups.
> 
> It still yields double for me:
> 
> $ r3.mapcalc 'map_g=g#map'
> $ r3.info map_g -t
> datatype="double"

You can't store CELL values in the output; they only exist internally. 
E.g.:

	r3.mapcalc 'map_g=g#map/2'

will only ever produce integers; you won't see e.g. 1.5 in the output,
as dividing two integers always produces an integer result.

Even with the above patch, writing an expression of type CELL to a G3D
map will result in a DCELL map (FCELL cannot accurately represent the
entire 32-bit integer range, as it only has a 24-bit mantissa).

-- 
Glynn Clements <glynn at gclements.plus.com>