[GIS] Writing Loop Script for GRASS

grasspython

I am new to Python and scripting in general.

My overall goal is to create a script that I can cycle through all the raster maps in a mapset and then save the univariate stats to a .csv file for each raster map. I am working on the first part: cycling through all the maps and printing the stats to the terminal window. Later, I will try and save the data to a .csv file. I have seven maps in my practice mapset and I get seven errors, but my problem is that I do not know what variable to assign to the parameter MAP (hopefully, I stated that correctly).

Here is what I have thus far:

import sys
import os
import grass.script as grass
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb="/Users/kc/GrassData"
location="Data_2013"
mapset="methods_check"
gsetup.init(gisbase, gisdb, location, mapset)
grass.run_command("g.list", _type="rast") 
mymaps = grass.parse_command("g.list", _type="rast") 
print len(mymaps) 
for items in mymaps:
    grass.run_command('r.univar', map='items', separator=",")

What should I put for map = 'items' so that it cycles through all the raster maps of a mapset?

I am using GRASS 7.0, but the same thing happens in GRASS 6.5. In Grass 6.5, the parameter separator must be changed to fs. Both versions label the input file "map".

Best Answer

You need to understand the difference between:

  • grass.run_command which performs a GRASS command
grass.run_command('g.remove', raster="raster1")
  • grass.read_command which is used to receive a result (generally a string, same result as the command in GRASS GIS, of little interest to Python, you need to use the re (Regular Expression) module to process the result)
grass.read_command("g.region",flags="p")  
'projection: 99 (Lambert Conformal Conic)\nzone:       0\ndatum:      bel72\nellipsoid:  international\nnorth:      112990.07863053\nsouth:      108258.30027318\nwest:       227552.40388837\neast:       235027.97855494\nnsres:      6.35138035\newres:      6.35138035\nrows:       745\ncols:       1177\ncells:      876865\n'
  • grass.parse_command that parse the result for Python (generally a dictionary)
grass.parse_command("g.region",flags="p")  
{'ellipsoid:  international': None, 'zone:       0': None, 'nsres:      6.35138035': None, 'south:      108258.30027318': None, 'ewres:      6.35138035': None, 'west:       227552.40388837': None, 'rows:       745': None, 'east:       235027.97855494': None, 'projection: 99 (Lambert Conformal Conic)': None, 'cells:      876865': None, 'north:      112990.07863053': None, 'cols:       1177': None, 'datum:      bel72': None}

So, in your case, it is not:

grass.run_command("g.list", _type="rast")

whose result is 0

but:

rasters = grass.read_command("g.list",_type="rast")

whose result is a string with all the layers in the variable rasters

type(raster)
<type 'str'>

But it is not interesting at all here, because you only need the names of the layers:

rasters = grass.parse_command("g.list", _type="rast")

whose result is a Python dictionary and for further, see Python Scripts For GRASS GIS: Using "parse_command()"

The result is a list with the layer names:

['raster1', 'raster2', 'raster3']

so the processing is:

for items in mymaps:
     grass.parse_command('r.univar', map=items)

.

Related Question