Contents
This tutorial is an introduction to geospatial programming in GRASS GIS with Python. In this tutorial you will learn how to automatically import, process, render maps of a time series of digital surface models from aerial surveys of the LSU Hilltop Arboretum. For this tutorial install GRASS GIS and download hilltop_drone_data.zip, extract it, and move it your GRASS database directory. Optionally install a text editor like Atom. With the GRASS Python Scripting Library you can:
There are two ways to run Python scripts in GRASS. In the Interactive Python Shell in the Layer Manager you can open a simple editor where you can write and run scripts. Or you can write your scripts in a text editor and then launch the script from the GRASS file menu.
GRASS Interactive Python Shell |
---|
GRASS scripting with a text editor |
---|
Start GRASS and open the Python tab in the Layer Manager. This an Interactive Python Shell. Try some Python here.
print("Hello World!")
Open the Console in the Layer Manager
and run a GRASS command.
Import a raster map with the module
r.in.gdal.
Set your input data file path to one of the rasters
in the hilltop_drone_data
directory.
Use the full file path.
r.in.gdal input=surface_2020_03_19.tif output=surface_2020_03_19
Open the Interactive Python Shell in the Layer Manager.
In the bottom left corner of the Interactive Python Shell
open the Simple Editor.
To run the same command in the Simple Editor using Python,
first add the Python shebang,
then import the
GRASS Python Scripting Library
with import grass.script as gscript
,
and use the run_command()
function from the grass.script
package
to run a GRASS module.
Use the full file path for the input map in the example below.
#!/usr/bin/env python3
import grass.script as gscript
gscript.run_command(
'r.in.gdal',
input='surface_2020_03_19.tif',
ouput='surface_2020_03_19',
overwrite=True)
Write a script that uses a for loop
to import all of the maps in a directory.
Start GRASS and
create a new location called laspm_hilltop
with EPSG code 6478 for NAD83(2011) / Louisiana South Meters
.
In the script first import
the os
library and the grass.script
library.
Then set the GRASS environment settings.
Assign the GRASS database directory
to the variable gisdbase
.
The use the os.path.join
method to set the filepath
to hilltop_drone_data
directory
inside of the GRASS database directory.
Use a for loop to iterate through all of the files in that directory.
For each file, use an if statement to test if it is a geotiff.
If it is then import it with the module
r.in.gdal.
#!/usr/bin/env python3
import os
import grass.script as gscript
# settings
env = gscript.gisenv()
overwrite = True
env['GRASS_OVERWRITE'] = overwrite
env['GRASS_VERBOSE'] = False
env['GRASS_MESSAGE_FORMAT'] = 'standard'
gisdbase = env['GISDBASE']
location = env['LOCATION_NAME']
mapset = env['MAPSET']
# set path
data = os.path.join(
gisdbase,
"hilltop_drone_data")
# set region
gscript.run_command('g.region', res=0.1)
# iterate through files in directory
for file in os.listdir(data):
filename = os.path.splitext(file)[0]
# iterate through geotiffs
if file.endswith('.tif'):
# check if raster is already in mapset
gscript.run_command('r.in.gdal',
input=os.path.join(data, file),
output=filename,
overwrite=overwrite)
else:
pass
To automatically process maps in GRASS,
create a list of maps with the module
g.list
then use a for loop to iterate through the list of maps,
running modules to process each map.
This example iterates through the digital surface models
that you just imported from hilltop_drone_data
to generate shaded relief maps for each.
#!/usr/bin/env python3
import os
import grass.script as gscript
# settings
env = gscript.gisenv()
overwrite = True
env['GRASS_OVERWRITE'] = overwrite
env['GRASS_VERBOSE'] = False
env['GRASS_MESSAGE_FORMAT'] = 'standard'
gisdbase = env['GISDBASE']
location = env['LOCATION_NAME']
mapset = env['MAPSET']
# list rasters in mapset
raster_list = gscript.list_grouped('rast',
pattern='surface_*')[mapset]
for raster in raster_list:
# set region
gscript.run_command(
'g.region',
raster=raster,
res=0.1)
# set color tables
gscript.run_command(
'r.colors',
map=raster,
color='viridis',
flags='e')
# compute hillshade
gscript.run_command(
'r.relief',
input=raster,
output='relief'+raster[-11:],
zscale=3,
overwrite=overwrite)
# compute shaded relief
gscript.run_command(
'r.shade',
shade='relief'+raster[-11:],
color=raster,
output='shaded_relief'+raster[-11:],
brighten=30,
overwrite=overwrite)
To automatically render maps from GRASS, create a list of maps with the module g.list, iterate through the list with a for loop, start a graphics display monitor with d.mon using the Cairo driver, add a map with d.rast, and stop the monitor. The maps in the graphics monitors will written directly to files in the current GRASS location.
#!/usr/bin/env python3
import os
import grass.script as gscript
# settings
env = gscript.gisenv()
overwrite = True
env['GRASS_OVERWRITE'] = overwrite
env['GRASS_VERBOSE'] = False
env['GRASS_MESSAGE_FORMAT'] = 'standard'
gisdbase = env['GISDBASE']
location = env['LOCATION_NAME']
mapset = env['MAPSET']
# list rasters in mapset
raster_list = gscript.list_grouped('rast',
pattern='shaded_relief_*')[mapset]
for raster in raster_list:
# set region
gscript.run_command(
'g.region',
raster=raster,
res=0.1)
# write map to image file
gscript.run_command(
'd.mon',
start="cairo",
width=1000,
height=1000,
output=os.path.join(
gisdbase,
location,
raster+'.png'),
overwrite=overwrite)
gscript.run_command(
'd.rast',
map=raster)
gscript.run_command(
'd.mon',
stop="cairo")
March 19, 2020 |
---|
April 25, 2020 |
---|
June 1, 2020 |
---|
June 30, 2020 |
---|
August 3, 2020 |
---|
August 29, 2020 |
---|