1

I'm trying to convert my GRASS/Bash script into the Phython Shell. Now I have a raster dataset (test.01) which I want to convert in 1 or 0. In GRASS I used the command r.mapcalc

r.mapcalc test.new.01="if(test.01/test==1,1,0)"

Now in the Python Shell I figured just the first steps, but I don't get the concept how to return the calculation into a new raster dataset.

import grass.script as grass
import grass.script.setup as gsetup
for rast in grass.list_strings(type = 'rast'):
    print rast
...
...
test.01@mapset

Now I tried to use the Python function grass.raster.mapcalc() but without success, it kills my GRASS GIS process.

expr = "if(test.01@mapset/test.01@mapset==1,1,0)"
output = grass.raster.mapcalc(expr)

The question is also how to write the result into a new raster dataset.

2

You need to understand the raster map algebra: any expression needs and input raster and a resulting raster and this layer is automatically added to the mapset.

In grass.raster.mapcalc() , the example is

expr1 = '"%s" = "%s" * 10' % (output_raster, input_raster)

If I define the input and output rasters as:

input_raster = "test.01@mapset"
output_raster = "test.02@mapset"

The expression is

expr1 = '"%s" = "%s" * 10' % (output, input)
print expr1
"test.02@mapset" = "test.01@mapset" * 10

In your expression expr = "if(test.01@mapset/test.01@mapset==1,1,0)" there is no resulting raster.

 expr = "if(test.01@mapset/test.01@mapset==1,1,0)"
 print expr
 if(test.01@mapset/test.01@mapset==1,1,0)

The expression you want is:

  • if the value of test.01@mapset/test.01@mapset = 1 is true (cell by cell) then the value of the cell in test.02@mapset is 1
  • if the value of test.01@mapset/test.01@mapset = 1 is false (cell by cell) then the value of the cell in test.02@mapset is 0

The correct solution would be

"test.02@mapset" = "if(test.01@mapset/test.01@mapset==1,1,0)"
# or
"test.02@mapset = if(test.01@mapset/test.01@mapset==1,1,0)"

But the result of your expression will always give 1 because,cell value by cell value, test.01@mapset/test.01@mapset = 1

  • some how it wants me to set quotation marks around the rasterdata name. expr1 = '"%s" = if("%s"/"%s"==1,1,0)' % (output, input, input) or "test.02@mapset" = if("test.01@mapset"/"test.01@mapset"==1,1,0). – wittich Sep 17 '15 at 12:32
  • Maybe of interest: github.com/wenzeslaus/python-grass-addon - shortcourse on Python in GRASS GIS as an ipython notebook (scroll down to see how to launch the course) – markusN Sep 18 '15 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.