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@virgin.net>
