Scripting map production in GRASS?

I'm trying to automate my production of maps with GRASS.

What I want to do can be described in fairly generic terms: I have two layers (one vector, one raster), which will be identical in every map. I have one layer (a raster heat-map) that will come from a different file for every map. I can write a shell or Python script (Mac OS X) to loop some GRASS command over all the heat-map files in a directory, but I am not sufficiently competent with GRASS to have figured out how to handle that end of the process.

Could someone give me an example of GRASS command-line code that will take a few files as inputs and produce a single image file as output?

Answer

For that, you need;

1. to import the files (shapefiles, raster) in GRASS GIS



2. use the adequate modules

You can use Bash or Python from:

• the Command Console or the Python Shell of the Layer manager

• the GRASS shell (Python here, but you can also use R)

• the Mac shell (Python, here, but you can also use R)

To use bash scripts, see Shell Scripting or On scripting GRASS GIS: Building location-independent command line tools
To use Python, see GRASS and Python or Python Scripts For GRASS GIS

I will develop Python that I know best.

Python with grass module

In Python, you simply access the GRASS functions as:

debligne =grass.read_command("v.to.db", flags="p", map="testgrass", type="line", option="start", units="meters" , quiet=True)

If you want to import all the tif files in a directory, see , for example, Python: script to import multiple LANDSAT images to Grass GIS

for dirpath, dirname, filenames in os.walk(dirpath):
    # Iterate through the files 
    for raster in filenames:
          # If the suffix is '.TIF', process
          if raster.upper().endswith('.tif'):
                # full path to your file
                full_path = os.path.join(dirpath, tif_file)

                # GRASS commands
                grass.message('Importing %s -> %s@%s...' % (full_path, tif_file, dirpath))
                grass.run_command('r.in.gdal',flags = 'o',input = full_path, output = tif_file,quiet = True,overwrite = True)

See other example in the above-mentioned references or in Automatic 3D geological boreholes representation (automate v.extrude from a table ?): my solution in Python or GRASS and the Python geospatial modules (Shapely, PySAL,...)

Python with osgeo module

You can also use the osgeo (GDAL/OGR) Python module:

from osgeo import ogr
# open grass vector layer
ds = ogr.Open('/Users/username/grassdata/geol/mymapset/vector/testgrass/head')

# vector layer  
layer = ds.GetLayer(0)
layer.GetName()
'testgrass'
 feat = layer.GetFeature(0) 
 # geometry
 geom = feat.GetGeometryRef() 
 geom.ExportToWkt() 
 'LINESTRING (186139.123704173340229 53082.654193976894021,188199.122798504744424 53467.758558732457459)'