Notebook

State-level Overview of Sewersheds¶

In this Jupyter notebook for Python, we visualize sewersheds at the state level.

Data requirements¶

Approximate positions of sewersheds represented as points
- Each point needs to have an attribute with size of the population the given sewershed serves.
State boundary

Software setup¶

We will use a couple of standard Python packages and GRASS GIS.

For now, the setup here assumes Linux. Instructions for Windows are available at GRASS GIS Jupyter notebooks wiki page.

In [ ]:

# Import Python standard library and IPython packages we need.
import os
import subprocess
import sys
from pathlib import Path

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

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

Get the data ready¶

This notebooks needs sewershed locations as vector points with at attributes for population size and name. The file should be in directory data/sewers_state and should be named sewer_points with file extension appropriate for the format, e.g. sewer_points.shp. Either rename the files or modify the code below if needed.

In [ ]:

data_directory = "data/sewers_state"
point_file = "sewer_points"  # Filename without file extension
state_boundary_file = "state_boundary"

In [ ]:

grass_project = "data/state_overview"

To compute the data, we will use a GRASS project (aka location).

In [ ]:

!grass -e -c $data_directory $grass_project

In [ ]:

session = gj.init(grass_project)

From now on, we will use the following variables to refer to the data in the GRASS project.

In [ ]:

point_vector = "points"
state_boundary_vector = "state_boundary_vector"

In [ ]:

gs.run_command("v.import", input=data_directory, layer=point_file, output=point_vector)

In [ ]:

gs.run_command(
    "v.import",
    input=data_directory,
    layer=state_boundary_file,
    output=state_boundary_vector,
)

Render¶

In [ ]:

# Extract the large sewersheds.
large_sewersheds = "large_sewersheds"

gs.run_command(
    "v.extract",
    input=point_vector,
    output=large_sewersheds,
    where="pop_2020 > 100000",
)

In [ ]:

gs.run_command("g.region", vector="nc_state", res=1000, grow=50)
state_map = gj.Map(width=1024, use_region=True)
state_map.d_background(color="white")
state_map.d_vect(flags="s", map=state_boundary_vector, fill_color="none")
state_map.d_vect_thematic(
    map=point_vector,
    column="pop_2020",
    legend_title="Population per sewershed",
    colors="#ADA8A7,#6A00A8,#B12A90,#E16462,#FCA636,#F0F921",
    breaks=[100000, 200000, 300000, 400000],
    boundary_color="none",
    icon="basic/circle",
    size=10,
)
# We rendered the points colored by population. Now we take the large
# sewersheds and render their names as labels in the map (this happens in two
# steps (first, labels are extracted and laid out, then rendered in
# the second step).
gs.run_command(
    "v.label.sa",
    map=large_sewersheds,
    column="site_name",
    labels="points_site_name",
    font="DejaVuSans",
    size=30000,
    color="black",
    isize=25000,
)
state_map.d_labels(labels="points_site_name")
state_map.d_legend_vect(flags="b", at=(3, 34), columns=2)
state_map.d_barscale(flags="n", at=(2, 7), length=200, units="kilometers")
state_map.save(Path(data_directory) / "nc_sewersheds.png")
state_map.show()