[GRASSLIST:2925] Re: Puzzling mapcalc command

[Glynn Clements]

SWlab wrote:

> > > I want to compute some maps with mapcalc (GRASS53 v2004/03.. on a SuSE
> > > box). Depending on the way I parse the info to the command I have results
> > > or none
> 
> > Actually, I'm slightly surprised that r.mapcalc actually made it to
> > the last line without reporting an error. In retrospect, I can see
> > why; however, r.mapcalc should probably be changed to prohibit certain
> > aspects of the above. In particular, it should probably prohibit using
> > an identifier on the LHS of an assignment after it has occurred
> > elsewhere.
> 
> Glynn, thanks a lot for your enlighting answer. For a bit more background:
> I'm porting a set of bash scripts to perl. The bash scripts were basically a 
> succession of calls to "r.mapcalc": maps were then evaluated at once, and 
> everything worked fine. Except of course that I kept on calling "r.mapcalc", 
> and I assume that there should be some initialization each time I call the 
> command (I don't speak C, so the sources are not really meaningful to me).

A multi-statement r.mapcalc program is different from calling
r.mapcalc multiple times. The former does:

	open all input maps
	open all output maps
	for each row
		read row from all input maps
		for each output map
			calculate and write row of output map
	close all input maps
	close all output maps

while the latter does

	for each output map
		open all input maps
		open output map
		for each row
			read row from all input maps
			calculate and write row of output map
		close output map
		close all input maps

The main semantic consequence of this is that, in a multi-statement
r.mapcalc program, you can't use the output maps as maps, only as
variables.

OTOH, if any maps occur on the RHS of more than one assignment,
performing multiple assignments in one command will be more efficient,
as the input maps are only opened/read/closed once per command, not
once per output map.

> So I thought it'd be quicker (and more elegant) to put all my calls into one 
> big string, and call mapcalc just once (with the small function below). Fatal 
> mistake as you explained it so well... It seems that a workaround would be to 
> create shortr command lines, making sure the maps are created before using 
> them...

For complex calculations, you ideally want a two-level structure, i.e. 
a sequence of multi-statement r.mapcalc programs, with each program
doing as much as possible to minimise the total number of programs.

Also, to avoid writing intermediate variables as maps, use eval(). 
E.g.

	V1 = 1
	V2 = 2
	Out = V1+V2

will create maps V1, V2 and Out, while

	Out = eval(V1 = 1, V2 = 2, V1+V2)

will only create "Out".

[eval() returns its last argument; the other arguments are only useful
insofar as they have side effects.]

BTW, you probably shouldn't use a variable on both the LHS and RHS of
an assignment (e.g. V=V+1). Whilst it will work in most cases, it
isn't guaranteed. Instead, use several variables, assigning each one
only once.

> Is there any alternative to "r.mapcalc" to create/modify rasters in scripts ? 
> I thought about using PDLs in Perl, but the interaction with GRASS has to be 
> optimized to be really useful (I started using r.in.ascii/r.out.ascii to 
> read/write GRASS maps to/from PDLs, but the process is a bit too long...). Of 
> course, not talking C is a handicap, and I'm not really keen on using  
> FORTRAN either... Any idea would be welcomed.

Note that the limitation of not being able to read a map until it has
been created is inherent in GRASS itself. When creating maps (which
includes overwriting existing maps), the data is always written to a
temporary file. The new map isn't "visible" until it has been closed
(i.e. completely written).

-- 
Glynn Clements <glynn.clements at virgin.net>