Using GRASS GIS in Jupyter Notebooks:
An Introduction to grass.jupyter

Caitlin Haedrich, Vaclav Petras, Anna Petrasova

grass_jupyter

NCSU GeoForAll Lab at the Center for Geospatial Analytics
North Carolina State University

Caitlin Haedrich

  • PhD student GeoForAll Lab at North Carolina State University
  • Advised by Helena Mitasova (and Anna Petrasova and Vashek Petras)
  • Interested in geovisualization and increasing accessibilty of GIS tools
  • Working with integrating GRASS GIS with Jupyter through Google Summer of Code 2021 and Mini Project Grant from GRASS GIS
  • TA for graduate-level Geospatial Computation and Simulation

geoforall

ncstate

Teaching Geospatial Analytics

oldGIS714

GIS714: geospatial computation and simulation

Graduate level course required for Geospatial Analytics PhD students

Previously taught with GRASS tutorials, instructions posted on webpage

Goal: convert to GRASS-based assignments in Jupyter Notebooks

GRASS GIS

  • Open-source GIS software with over 500 tools and 400 addons
  • Interfaces: graphical user interface, command line, Python, C
  • 3rd Party Interfaces: R, QGIS, actinia (REST API)

GUI

grassgis

Project Jupyter

“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

Jupyter

The default kernel is IPython.

In [1]:
import sys

v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")

We are using Python 3.9.5

Output can be interactive.

In [2]:
from ipywidgets import interact

def f(x):
    return x

interact(f, x=10);

grass.jupyter

...enhances the existing GRASS Python API to allow Jupyter Notebook users to easily manage GRASS data, visualize data including spatio-temporal datasets and 3D visualizations

  • Released with main GRASS GIS distribution starting in version 8.2
  • Python package that makes GRASS more accessible in Jupyter by providing:
    • Data management functions
    • Visualizations classes for viewing and interacting with GRASS through in-line figures
  • Created with syntax consistent with GRASS's existing Python API and command line interface

You can also try it out in Binder:

Binder

Teaching Geospatial Computation and Simulation with grass.jupyter

GIS714: Geospatial Computation and Simulation, Spring 2022


Course Repo on GitHub

Binder

We'll use an example from the course on water simulation...

Getting Started

In [3]:
# Import Python standard library and IPython packages we need.
import subprocess
import sys

# Ask GRASS GIS where its Python packages are.
sys.path.append(
    subprocess.check_output(["grass80", "--config", "python_path"],
                            text=True, shell=True).strip()
)

# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj

# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "HW3_water_simulation")

# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")

gj.init starts session, sets GRASS environmental variables for notebooks

And then we are ready to begin using GRASS!

Create a map of flooded area¶

We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:

In [4]:
gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5,
               lake="flood1", coordinates="638728,220278")

# See results
flood1 = gj.Map(use_region=True, width=200)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()
Out[4]:
Question 4¶

Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.

In [5]:
#### Your Answer Here

We can also make 3D maps using similar syntax

In [6]:
elev_map = gj.Map3D()
elev_map.render(elevation_map="elevation", color_map="elevation",
                perspective=20)
elev_map.overlay.d_legend(raster="elevation", at=(60, 97, 87, 92))
elev_map.show()
Out[6]:

Syntax: Integrating GRASS Display Modules with grass.jupyter

modules

We can call display modules using Python magic short-cut (__getattr__):

To add a raster, we need d.rast so we use Map.d_rast(), Map3D.overlay.d_rast

Estimating inundation extent using HAND methodology¶

We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011).

For this section, we will change our computation region to elevation which is a larger study area than we used above. We use r.watershed to compute the flow accumulation, drainage and streams.

In [7]:
gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc",
               drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k",
               type="line")

Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.

In [8]:
gs.run_command("r.stream.distance", stream_rast="streams_100k",
               direction="drainage", elevation="elevation",
               method="downstream", difference="above_stream")

With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...

In [9]:
gs.run_command("r.lake.series", elevation="above_stream",
               start_water_level=0, end_water_level=5, water_level_step=0.5, 
               output="inundation", seed_raster="streams_100k")

... or visualize the flood with TimeSeriesMap.

In [ ]:
flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.d_rast(map="elevation")
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(2,30,2,5)) #Add legend
flood_map.show()

TimeSeriesMapOutput

flood-animation

folium Integration for Interactive Maps

Interactivity is often useful for explore geospatial data.

folium...

  • is a popular Python library for creating Leaflet.js maps
  • creates interactive HTML maps that can be displayed in-line
  • has built-in tilesets from OpenStreetMap, Mapbox, and Stamen, and supports custom tilesets

foliumLogo

Moving data from GRASS GIS location projection to folium is a challenge

folium is projected in Web Mercator (EPSG:3857)

However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)


We pass data to folium by reprojecting to temporary locations

grass_folium

In [ ]:
import folium 

# Create figure
fig = folium.Figure(width=600, height=400)

# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271],
               zoom_start=13)

# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)

# make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
    [35.781608,-78.675800],
    popup="<i>Center For Geospatial Analytics</i>",
    tooltip=tooltip
).add_to(m)

# and a circle
folium.Circle(
    radius=120,
    location=[35.769781,-78.663160],
    popup="Great Picnic Area",
    color="crimson",
    fill=False,
    tooltip=tooltip
).add_to(m)

# Add the map to the figure
fig.add_child(m)
In [12]:
# Display figure
fig
Out[12]:

The InteractiveMap class allows users unfamiliar with folium to produce maps easily.

In [13]:
# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()
Out[13]: