[GRASS-user] How to initialize raster maps.

Marcello_Gorini · December 8, 2011, 9:26pm

Dear all,

Inside a loop in a shell script, I keep incrementing the values of a raster
map using r.mapcalc with some logical expressions.

Therefore, I need to initialize the raster map before the loop itself.
Ideally, I would like to initialize it as a raster with the same dimensions
as my base map, but containing only zeros.

In order to do that, I use the following expression:

r.mapcalc "initialized_map=if(base_map>999999,1,0)"

Besides being one of the ugliest pieces of code one could ever write, it
also incurs in waste of processing time.

Which would be a more elegant and optimized way of initializing these raster
maps?

Thanks in advance for any thoughts you could share.

Best regards,
Marcello.

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/How-to-initialize-raster-maps-tp7075938p7075938.html
Sent from the Grass - Users mailing list archive at Nabble.com.

mlennert · December 8, 2011, 10:04pm

On 08/12/11 22:26, Marcello Gorini wrote:

Dear all,

Inside a loop in a shell script, I keep incrementing the values of a raster
map using r.mapcalc with some logical expressions.

Therefore, I need to initialize the raster map before the loop itself.
Ideally, I would like to initialize it as a raster with the same dimensions
as my base map, but containing only zeros.

In order to do that, I use the following expression:

r.mapcalc "initialized_map=if(base_map>999999,1,0)"

Besides being one of the ugliest pieces of code one could ever write, it
also incurs in waste of processing time.

Which would be a more elegant and optimized way of initializing these raster
maps?

Thanks in advance for any thoughts you could share.

g.region rast=YourBaseMap
r.mapcalc initialized_map=0

Moritz

Marcello_Gorini · December 8, 2011, 11:42pm

how difficult, isn’t it?

Thanks!

But there is no way to avoid r.mapcalc, is there?

On Thu, Dec 8, 2011 at 8:04 PM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 08/12/11 22:26, Marcello Gorini wrote:

Dear all,

Inside a loop in a shell script, I keep incrementing the values of a raster
map using r.mapcalc with some logical expressions.

Therefore, I need to initialize the raster map before the loop itself.
Ideally, I would like to initialize it as a raster with the same dimensions
as my base map, but containing only zeros.

In order to do that, I use the following expression:

r.mapcalc “initialized_map=if(base_map>999999,1,0)”

Besides being one of the ugliest pieces of code one could ever write, it
also incurs in waste of processing time.

Which would be a more elegant and optimized way of initializing these raster
maps?

Thanks in advance for any thoughts you could share.

g.region rast=YourBaseMap
r.mapcalc initialized_map=0

Moritz

Glynn_Clements1 · December 9, 2011, 10:14pm

Marcello Gorini wrote:

>> Inside a loop in a shell script, I keep incrementing the values of a
>> raster map using r.mapcalc with some logical expressions.

Note that using the same map as both input and output results in
undefined behaviour. Even if it happens to work, there's no guarantee
that it will continue to work in future versions.

If you need to implement an iterative algorithm, you should generate a
new map, then replace the original with g.remove+g.rename afterwards.

>> Ideally, I would like to initialize it as a raster with the same
>> dimensions as my base map, but containing only zeros.
>>
>> In order to do that, I use the following expression:
>>
>> r.mapcalc "initialized_map=if(base_map>**999999,1,0)"

There is no point in referring to your original map in the r.mapcalc
expression; the result will always be created with the bounds and
resolution of the current region. AFAIK, the only modules which don't
behave this way are r.in.*, which create the output to match the size
of the input file.

But there is no way to avoid r.mapcalc, is there?

I'm sure that there are other ways to create a map consisting of
zeros, but there's no reason to suspect that any of them will be
faster than r.mapcalc.

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

Marcello_Gorini · December 9, 2011, 11:57pm

Glynn:

Note that using the same map as both input and output results in
undefined behaviour. Even if it happens to work, there’s no guarantee
that it will continue to work in future versions.

If you need to implement an iterative algorithm, you should generate a
new map, then replace the original with g.remove+g.rename afterwards.

Now, you scared me, because my ongoing work is pretty much based on this kind of code. I use r.mapcalc to create maps that serve as counters such as the command i++ in C. I initialize them as zero maps and then based on some logic I keep incrementing them, pretty much like the following pseudo-code to implement a sum:

sum = 0
for i = 1 to 10
sum = sum + i
end

But I do more that that.

For instance, if certain conditions are met in other maps at a given iteration and the initialized map is zero, then update the initialized map with the value of “i”. If not, don’t change it. In this way, my map is updated every iteration until no more cells meet the condition or some threshold is met and processing stops. I even use d.rast to display the map as it is being updated.

How should I proceed to avoid the “undefined behaviour” and still get what I want? (I hope I was able to make myself understood )

Thanks for your comments.

Glynn:

There is no point in referring to your original map in the r.mapcalc
expression; the result will always be created with the bounds and
resolution of the current region. AFAIK, the only modules which don’t
behave this way are r.in.*, which create the output to match the size
of the input file.

Yes, Moritz’s suggestion showed me that. I perfectly understand it now, thanks.

Glynn:

I’m sure that there are other ways to create a map consisting of
zeros, but there’s no reason to suspect that any of them will be
faster than r.mapcalc.

There really isn’t. I just suspected it was taking too long and found out in my little experiment that was only my fault, putting a logic expression where there was no need of it. Setting the region and r.mapcalc zero_map = 0 does the job.

Thank you.

Best,
Marcello.

mlennert · December 10, 2011, 2:25pm

On 10/12/11 00:57, Marcello Gorini wrote:

Glynn:

    Note that using the same map as both input and output results in
    undefined behaviour. Even if it happens to work, there's no guarantee
    that it will continue to work in future versions.

    If you need to implement an iterative algorithm, you should generate a
    new map, then replace the original with g.remove+g.rename afterwards.

Now, you scared me, because my ongoing work is pretty much based on this
kind of code. I use r.mapcalc to create maps that serve as counters such
as the command i++ in C. I initialize them as zero maps and then based
on some logic I keep incrementing them, pretty much like the following
pseudo-code to implement a sum:

sum = 0
for i = 1 to 10
      sum = sum + i
end

But I do more that that.

For instance, if certain conditions are met in other maps at a given
iteration and the initialized map is zero, then update the initialized
map with the value of "i". If not, don't change it. In this way, my map
is updated every iteration until no more cells meet the condition or
some threshold is met and processing stops. I even use d.rast to display
the map as it is being updated.

How should I proceed to avoid the "undefined behaviour" and still get
what I want? (I hope I was able to make myself understood )

As Glynn suggests:

for i = 1 to 10
      sum_new = sum + i
      g.remove rast=sum
      g.rename rast=sum_new,sum

or

for i = 1 to 10
sum_new = sum + i
g.rename rast=sum_new,sum --overwrite

Glynn, is there any reason not to do use g.rename with --overwrite ?

Moritz

Marcello_Gorini · December 10, 2011, 5:05pm

Moritz:

As Glynn suggests:

for i = 1 to 10

sum_new = sum + i
g.remove rast=sum
g.rename rast=sum_new,sum

OK, sorry, I got so scared I overlooked Glynn’s suggestion. I did another experiment to see the differences in processing time. Since I iterate 1260 times in my script, I tried the following alternatives with “i” varying from 1 to 1260:

1 - r.mapcalc “sum=sum+$i” # my suspicious way

2 - r.mapcalc “sum_new = sum + $i” # Glynn’s suggestion after Moritz’s explanation
g.remove rast=sum --q
g.rename rast=sum_new,sum --q

3 - r.mapcalc “sum_new = sum + $i” # Moritz’s suggestion
g.rename rast=sum_new,sum --overwrite --q

Results:

1 - 125 seconds
2 - 183 seconds
3 - 154 seconds

Bottom line:

If Glynn says Moritz suggestion with the g.rename --overwrite is OK, I will give back the 30 seconds Moritz had saved me before and end up with 30 seconds more of processing time, but with no undefined behaviour.

If method number 2 must be used, I will end up with 60 seconds more, what won’t do me any harm anyway.

Thanks a lot for this interesting exercise.

Cheers,
Marcello.

Glynn_Clements1 · December 11, 2011, 6:37pm

Moritz Lennert wrote:

Glynn, is there any reason not to do use g.rename with --overwrite ?

Not that I know of. I had forgotten that g.rename supported
--overwrite.

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

Glynn_Clements1 · December 11, 2011, 6:53pm

Marcello Gorini wrote:

I use r.mapcalc to create maps that serve as counters such as
the command i++ in C.

Raster maps aren't supposed to be modified in-place; that applies to
all raster commands, not just r.mapcalc.

At the moment, you will usually get away with it, because:

1. the actual raster data is written to a temporary file, which
replaces the original file when the map is closed, and

2. the support files (cellhd, categories, colour table, etc) are
normally opened, read into (or written from) memory, then closed,
rather than being held open.

However: you won't necessarily get away with it if you're using
GDAL-linked maps (r.external, r.external.out), as those are likely to
be overwritten in place.

If we eventually change the raster format in 7.x, there's a chance
that older formats will be supported by this mechanism. Even if we
support multiple versions natively, fixing issues with concurrent
access is on the agenda, and that may end up breaking in-place
modification.

IOW, when I say that in-place modification is "not guaranteed" to
work, I mean that there are reasons to suspect that it may actually
fail in specific cases now, and in the general case in the future.

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