This notebook introduces the grass.jupyter package which simplifies the usage of GRASS GIS in Jupyter Notebook.
The grass.jupyter package was initially written as part of Google Summer of Code in 2021 by Caitlin Haedrich and was experimentally included in version 8.0.0. Caitlin further improved it thanks to the GRASS Mini Grant 2022. The package was officially released for the first time as part of version 8.2.0. Credits for mentoring and additional development go to Vaclav Petras, Helena Mitasova, Stefan Blumentrath, and Anna Petrasova as well as to many members of the GRASS community who provided important feedback.
In addition to simplifying the launch of GRASS GIS with a dedicated init function, grass.jupyter has two main display classes, Map and InteractiveMap. Using the GRASS rendering engine in the background, Map creates maps as PNG images. InteractiveMap displays GRASS rasters and vectors with folium, a leaflet library for Python. The package includes also Map3D and TimeSeriesMap.
This interactive notebook is available online thanks to the Binder service. To run the select part (called a cell), hit Shift + Enter
.
import os
import subprocess
import sys
# Ask GRASS GIS where its Python packages are.
sys.path.append(
subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)
# Import GRASS packages
import grass.script as gs
import grass.jupyter as gj
# Start GRASS Session
session = gj.init("~/data/grassdata", "nc_basic_spm_grass7", "user1")
# Set computational region to the elevation raster.
gs.run_command("g.region", raster="elevation")
The Map
class creates and displays GRASS maps as PNG images. There are two ways to add elements to the display. First, the name of the GRASS display module can be called as an attribute by replacing the "." with "_" in the module name. For example:
m = Map()
m.d_rast(map="elevation")
Alternatively, GRASS display modules can be called with the run()
method:
m = Map()
m.run("d.rast", map="elevation")
To display the image, call show()
.
# Create Map instance
example_map = gj.Map()
# Add a raster, vector and legend to the map
example_map.d_rast(map="elevation")
example_map.d_vect(map="streams")
example_map.d_legend(raster="elevation", at=(55, 95, 80, 84), flags="b")
# Display map
example_map.show()
We can also have multiple instances of Map
. Here, we create another map then go back and modify the first map.
# Make a second instance.
# Just for variety, we'll make this one a different size
small_map = gj.Map(height=200, width=220)
# Add a some layers
# We can also add layers with the run() methods
small_map.run("d.rast", map="elevation_shade")
small_map.run("d.vect", map="roadsmajor")
# Display second map
small_map.show()
# Then, we return to the first instance and continue to modify and display it
# Notice that layers a drawn in the order they are added
example_map.run("d.vect", map = "zipcodes", color="red", fill_color="none")
example_map.show()
By default the display extent (and resolution if applicable) is derived from the first raster or vector layer:
nc_map = gj.Map()
nc_map.d_vect(map="boundary_state")
nc_map.d_rast(map="geology")
nc_map.show()
To respect computational region, set use_region=True
:
geol_map = gj.Map(use_region=True)
geol_map.d_vect(map="boundary_state")
geol_map.d_rast(map="geology")
geol_map.show()
You can also use a saved region:
gs.run_command("g.region", save="myregion", n=224000, s=222000, w=633500, e=637300)
myregion_map = gj.Map(saved_region="myregion")
myregion_map.d_rast(map="elevation")
myregion_map.d_rast(map="lakes")
myregion_map.show()
# Create Interactive Map
raleigh_map = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
raleigh_map.add_raster("elevation")
raleigh_map.add_vector("roadsmajor")
raleigh_map.add_layer_control(position = "bottomright")
raleigh_map.show()
To share or embed the map in a website, we can export it has an HTML file.
raleigh_map.save(filename="raleigh_map.html")
We can also pass GRASS rasters and vectors directly to folium with the Raster and Vector classes. This provides much more flexibility when creating maps since we can access all of folium's capabilities.
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)
# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
[35.781608,-78.675800], popup="<i>Point of Interest</i>", tooltip=tooltip
).add_to(m)
# Add the map to the figure
fig.add_child(m)
# Display figure
fig
The Map3D
class creates 3D visualizations as PNG images. The m.nviz.image module is used in the background and the function render()
accepts parameters of this module.
The Map3D
objects have overlay
attribute which can be used in the same way as Map
and 2D images on top of the 3D visualization.
To display the image, call show()
.
First, let's create the object:
elev_map = gj.Map3D()
Now, render a 3D visualization of an elevation raster as a surface colored using, again, the elevation raster:
elev_map.render(elevation_map="elevation", color_map="elevation", perspective=20)
To add a raster legend on the image as an overlay using the 2D rendering capabilities accessible with overlay.d_legend
:
elev_map.overlay.d_legend(raster="elevation", at=(60, 97, 87, 92))
Finally, we show
elev_map.show()
Now, let's color the elevation surface using a landuse raster (note that the call to render
removes the result of the previous render
as well as the current overlays):
elev_map.render(elevation_map="elevation", color_map="landuse", perspective=20)
elev_map.show()
The init
function returns a reference to a session object which can be used to manipulate the current session. The session is global, i.e., the global state of the environment is changed. The session object is a handle for accessing this global session. When the kernel for the notebooks shuts down or is restarted, the session ends automatically. The session can be explicitly ended using session.finish()
, but that's usually not needed in notebooks.
Additionally, the session object can be used to change the current mapset. Here, we will switch to mapset called PERMANENT:
session.switch_mapset("PERMANENT")
Now we could add more data to the PERMANENT mapset or modify the existing data there. We don't need to do anything there, so we switch back to the mapset we were in before:
session.switch_mapset("user1")