Geoprocessing with GRASS GIS and Jupyter Notebooks

This tutorial demonstrates how to use GRASS GIS in Jupyter Notebooks.
geospatial
GRASS v8.5
jupyter
lidar
STAC
Author

Corey T. White

Published

August 22, 2024

GRASS v8.5 (Preview)

  • GRASS GIS is a Geospatial Processing Engine
  • Open Source (GPL v2)
  • Developed by International and Multi-institutional groups and individuals (GRASS Development Team)
  • Member of the Open Source Geospatial Foundatispace
  • Recieved Open Source Security Foundation (OpenSSF) Best Practices Badge - 2024

Modern Tooling

  • Jupyter Notebooks
  • GRASS GIS Python API
  • Actinia REST API
  • Easy intergration with other Data Science tools in Python and R ecosystem

Community

  • Active community of users and developers
  • Mailing lists, chat, and forums
  • Conferences and workshops
  • Documentation and tutorials
  • Mentoring and outreach programs

Leader in Open Science

  • Open Access
  • Community Mantaince and Support
  • Reproducible Research
  • Citations

Setup Environment

Import python packages and set up the GRASS GIS environment.

# import standard Python packages
import os
import sys
import subprocess
import json
from io import StringIO
from pathlib import Path

%reload_ext autoreload
%autoreload 2
# Visulization packages
import pandas as pd
import geopandas as gpd
import seaborn as sns
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)
# import GRASS GIS python packages
%reload_ext autoreload
%autoreload 2
import grass.script as gs
import grass.jupyter as gj
# create a temporary folder where to place our GRASS project
import tempfile
tempdir = tempfile.TemporaryDirectory()
print(tempdir.name)
/tmp/tmpnp2dcvk3

Area of Interest

Chaco Cultural National Historical Park

NPS Park Image

The full dataset is available through OpenTopography

Data Summary

  • Funding: National Science Foundation (NSF) Earth Sciences (EAR) Instrumentation and Facilities (IF) Program
  • Partner: University of New Mexico
  • Collector: National Center for Airborne Laser Mapping (NCALM)

Data Characteristics

  • Area 542.72 km^2
  • Over 13 Billion Points
  • Point Density 25.56 pts/m^2

Raster Resolution 0.5 m

Coordinate System: Horizontal: NAD83 (2011) (EPOCH:2010) / UTM Zone 13N Meters [EPSG: 6342] Vertical: NAVD88 [EPSG: 5703]

Units: Meters

Complete metadata can be found at: OpenTopography

Dorshow, W. (2019). 3D Landscape Reconstruction and Land Use Modeling, Chaco Canyon, NM 2016. National Center for Airborne Laser Mapping (NCALM). Distributed by OpenTopography. https://doi.org/10.5069/G9XG9P8D.. Accessed: 2024-08-19

Download Data

The data used in this tutorial is available through GitHub.

Lidar Data

Visualize the point cloud

from IPython.display import IFrame

# URL of the website to be embedded
url = 'https://ot-process2.sdsc.edu/potree/index.html?t=%5B233574.5,3994716,2188.5%5D&p=%5B234704.19367662142,3993579.3112938125,2825.660091521463%5D&r=%22https://ot-process2.sdsc.edu/appEntwineEPTService1724096588005642548939/pc1724096516919%22&m=9&era=%5B1858,2519%5D'
# Dimensions of the IFrame
width = 800
height = 600
# Display the IFrame in the notebook
IFrame(url, width=width, height=height)

GRASS Project Setup

Create a new project in GRASS for Chaco Culture National Historical Park

gs.create_project(path=tempdir.name, name="ChacoCanyon2016", epsg="6342", overwrite=True)
# start GRASS in the recently created project
session = gj.init(Path(tempdir.name,"ChacoCanyon2016"))

Download Add-ons

GRASS GIS Addons Over 400 add-ons available!

Let’s download the add-ons for the project.

with open("extensions.txt", "r") as f:
    lines = f.readlines()
    for line in lines:
        line = line.strip()
        print(f"Installing: {line}")
        gs.run_command("g.extension", extension=line, operation="add")

Prepare Data

Let’s examine point cloud data first with pdal before importing into GRASS GIS.

!pdal info --summary metadata/points2.laz
{
  "file_size": 18912844,
  "filename": "metadata/points2.laz",
  "now": "2024-08-22T17:20:40-0400",
  "pdal_version": "2.3.0 (git-version: Release)",
  "reader": "readers.las",
  "summary":
  {
    "bounds":
    {
      "maxx": 229974.27,
      "maxy": 3997397.29,
      "maxz": 2030.95,
      "minx": 229377.17,
      "miny": 3996851.14,
      "minz": 1846.72
    },
    "dimensions": "X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, ScanAngleRank, UserData, PointSourceId, GpsTime, ScanChannel, ClassFlags",
    "num_points": 9017333,
    "srs":
    {
      "compoundwkt": "COMPD_CS[\"NAD83(2011) / UTM zone 13N + NAVD88 height\",PROJCS[\"NAD83(2011) / UTM zone 13N\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-105],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6342\"]],VERT_CS[\"NAVD88 height\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"5703\"]]]",
      "horizontal": "PROJCS[\"NAD83(2011) / UTM zone 13N\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-105],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6342\"]]",
      "isgeocentric": false,
      "isgeographic": false,
      "prettycompoundwkt": "COMPD_CS[\"NAD83(2011) / UTM zone 13N + NAVD88 height\",\n    PROJCS[\"NAD83(2011) / UTM zone 13N\",\n        GEOGCS[\"NAD83(2011)\",\n            DATUM[\"NAD83_National_Spatial_Reference_System_2011\",\n                SPHEROID[\"GRS 1980\",6378137,298.257222101,\n                    AUTHORITY[\"EPSG\",\"7019\"]],\n                AUTHORITY[\"EPSG\",\"1116\"]],\n            PRIMEM[\"Greenwich\",0,\n                AUTHORITY[\"EPSG\",\"8901\"]],\n            UNIT[\"degree\",0.0174532925199433,\n                AUTHORITY[\"EPSG\",\"9122\"]],\n            AUTHORITY[\"EPSG\",\"6318\"]],\n        PROJECTION[\"Transverse_Mercator\"],\n        PARAMETER[\"latitude_of_origin\",0],\n        PARAMETER[\"central_meridian\",-105],\n        PARAMETER[\"scale_factor\",0.9996],\n        PARAMETER[\"false_easting\",500000],\n        PARAMETER[\"false_northing\",0],\n        UNIT[\"metre\",1,\n            AUTHORITY[\"EPSG\",\"9001\"]],\n        AXIS[\"Easting\",EAST],\n        AXIS[\"Northing\",NORTH],\n        AUTHORITY[\"EPSG\",\"6342\"]],\n    VERT_CS[\"NAVD88 height\",\n        VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n            AUTHORITY[\"EPSG\",\"5103\"]],\n        UNIT[\"metre\",1,\n            AUTHORITY[\"EPSG\",\"9001\"]],\n        AXIS[\"Gravity-related height\",UP],\n        AUTHORITY[\"EPSG\",\"5703\"]]]",
      "prettywkt": "PROJCS[\"NAD83(2011) / UTM zone 13N\",\n    GEOGCS[\"NAD83(2011)\",\n        DATUM[\"NAD83_National_Spatial_Reference_System_2011\",\n            SPHEROID[\"GRS 1980\",6378137,298.257222101,\n                AUTHORITY[\"EPSG\",\"7019\"]],\n            AUTHORITY[\"EPSG\",\"1116\"]],\n        PRIMEM[\"Greenwich\",0,\n            AUTHORITY[\"EPSG\",\"8901\"]],\n        UNIT[\"degree\",0.0174532925199433,\n            AUTHORITY[\"EPSG\",\"9122\"]],\n        AUTHORITY[\"EPSG\",\"6318\"]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"latitude_of_origin\",0],\n    PARAMETER[\"central_meridian\",-105],\n    PARAMETER[\"scale_factor\",0.9996],\n    PARAMETER[\"false_easting\",500000],\n    PARAMETER[\"false_northing\",0],\n    UNIT[\"metre\",1,\n        AUTHORITY[\"EPSG\",\"9001\"]],\n    AXIS[\"Easting\",EAST],\n    AXIS[\"Northing\",NORTH],\n    AUTHORITY[\"EPSG\",\"6342\"]]",
      "proj4": "+proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs",
      "units":
      {
        "horizontal": "metre",
        "vertical": "metre"
      },
      "vertical": "VERT_CS[\"NAVD88 height\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"5703\"]]",
      "wkt": "PROJCS[\"NAD83(2011) / UTM zone 13N\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-105],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6342\"]]"
    }
  }
}

Remove Outliers from Point Cloud

!pdal pipeline pdal/preprocessing.json

Import Data into GRASS GIS

Let’s calculate the point density of the point cloud data at 1m resolution.

gs.run_command('r.in.pdal',
            input='metadata/points2_clean.laz', 
            output='points_n',
            method='n', # Count number of points per cell
            resolution=1, # 1 meter
            flags="ewn",
            overwrite=True)

Visulalize the Raster Data

gs.run_command("r.colors", map="points_n", color="bcyr", flags="e")
m = gj.Map()
m.d_rast(map="points_n")
m.d_legend(raster="points_n", at=(5, 10, 50, 90), flags="b")
m.d_barscale(at=(5, 6), flags="n")
m.show()

points_n_info = gs.parse_command('r.info', map='points_n', format="json")
points_n_df = pd.DataFrame(points_n_info)

points_n_df.head(1)
north south nsres east west ewres rows cols cells datatype ... creator title timestamp units vdatum semantic_label source1 source2 description comments
0 3997398 3996851 1 229975 229377 1 547 598 327106 CELL ... coreywhite Raw X,Y,Z data binned into a raster grid by ce... None None None None metadata/points2_clean.laz generated by r.in.pdal r.in.pdal --overwrite -w -e -n input="metadata...

1 rows × 30 columns

Histograms of Raster Data

hist = gj.Map()
hist.d_histogram(map="points_n", flags="c")
hist.show()

univar_json = gs.parse_command('r.univar', map='points_n', format="json")
univar_df = pd.DataFrame(univar_json)
univar_df.head()
n null_cells cells min max range mean mean_of_abs stddev variance coeff_var sum
0 327106 0 327106 0 135 135 27.567006 27.567006 8.416197 70.832372 30.529964 9017333 
gs.run_command('r.in.pdal',
            input='metadata/points2_clean.laz', 
            output='points_median',
            method='median', # median of the z values
            resolution=1, # meter
            flags="ewn",
            overwrite=True)

gs.run_command("r.colors", map="points_median", color="elevation", flags="")
m = gj.Map()

m.d_rast(map="points_median")
m.d_legend(raster="points_median", at=(60, 95, 85, 90), flags="bd")
m.d_barscale(at=(5, 6), flags="n")
m.show()

Now let’s look at the histograms of the median elevation.

hist = gj.Map()
hist.d_histogram(map="points_median")
hist.show()

univar_json = gs.parse_command('r.univar', map='points_median', format="json", flags="e")
univar_df = pd.DataFrame(univar_json)
univar_df.head()
n null_cells cells min max range mean mean_of_abs stddev variance coeff_var sum first_quartile median third_quartile percentiles
0 327090 16 327106 1846.920044 1918.099976 71.179932 1896.1288 1896.1288 17.64089 311.200993 0.930363 6.202048e+08 1889.959961 1902.719971 1908.689941 [{'percentile': 90, 'value': 1911.6400146484375}] 
def spatial_resolution_analysis():
    output_maps = []
    for i in [0.5, 1, 3, 5, 10]:
        output_map = f'points_{i}m_mean'
        gs.run_command('r.in.pdal',
                input='metadata/points2_clean.laz', 
                output=output_map,
                method='mean', # mean of the z values
                resolution=i, # meter
                flags="ewn",
                overwrite=True)
        univar_json = gs.parse_command('r.univar', map=output_map, format="json", flags="e")
        univar_json[0]['resolution'] = i
        output_maps.append(univar_json[0])

    return pd.DataFrame(output_maps)

mean_univar_df = spatial_resolution_analysis()
mean_univar_df.set_index('resolution', inplace=True)
mean_univar_df.head(5)
n null_cells cells min max range mean mean_of_abs stddev variance coeff_var sum first_quartile median third_quartile percentiles
resolution
0.5 1298570 7565 1306135 1846.883301 1962.815796 115.932495 1896.145970 1896.145970 17.622542 310.554003 0.929387 2.462278e+09 1889.972046 1902.729980 1908.694946 [{'percentile': 90, 'value': 1911.6434326171875}]
1.0 327090 16 327106 1846.941772 1926.239624 79.297852 1896.130754 1896.130754 17.641135 311.209630 0.930375 6.202054e+08 1889.958130 1902.725830 1908.695435 [{'percentile': 90, 'value': 1911.6436767578125}]
3.0 36599 1 36600 1847.095459 1917.331787 70.236328 1896.056648 1896.056648 17.732032 314.424974 0.935206 6.939378e+07 1889.824585 1902.695557 1908.697388 [{'percentile': 90, 'value': 1911.6749267578125}]
5.0 13200 0 13200 1847.404785 1917.161377 69.756592 1896.171552 1896.171552 17.628175 310.752564 0.929672 2.502946e+07 1890.059448 1902.774170 1908.708374 [{'percentile': 90, 'value': 1911.650634765625}]
10.0 3355 0 3355 1847.787354 1916.982178 69.194824 1895.952356 1895.952356 17.820593 317.573521 0.939928 6.360920e+06 1889.712280 1902.729980 1908.641113 [{'percentile': 90, 'value': 1911.5826416015625}] 
import seaborn as sns
import matplotlib.pyplot as plt
sns.set_theme(style="darkgrid")

# Example plot: scatter plot of mean vs. resolution
# Adjust the column names according to your DataFrame structure
plt.figure(figsize=(10, 6))
sns.scatterplot(data=mean_univar_df, x='resolution', y='mean')

# Add titles and labels
plt.title('Mean Value vs. Resolution')
plt.xlabel('Resolution (m)')
plt.ylabel('Mean Value')

# Show the plot
plt.show()

Now let’s reset our computational region to run our analysis.

gs.run_command("g.region", raster="points_median", flags="pa")
projection: 1 (UTM)
zone:       13
datum:      nad83_2011
ellipsoid:  grs80
north:      3997398
south:      3996851
west:       229377
east:       229975
nsres:      1
ewres:      1
rows:       547
cols:       598
cells:      327106

Create a Digital Terrian Model (DTM)

Import lidar data as a vector in GRASS.

gs.run_command('v.in.pdal',
            input='metadata/points2_clean.laz', 
            output='lidar_points_be',
            class_filter=2, # Bare earth points
            flags="w",
            overwrite=True)

How many bare earth ppints did we just import?

# Assuming gs.parse_command is already defined and imported
lidar_be_info = gs.read_command('v.info', map='lidar_points_be', format="json")
lidar_be_dict = json.loads(lidar_be_info)
num_be_points = lidar_be_dict["points"]
print(f"We just imported {num_be_points:,} bare earth points")
We just imported 6,631,980 bare earth points

We will now interpolate the lidar points into our digital terrain model (DTM).

gs.run_command("v.surf.rst",
               input="lidar_points_be",
               elevation="lidar_be",
               slope="lidar_be_slope",
               aspect="lidar_be_aspect",
               pcurvature="lidar_be_pcurvature",
               tcurvature="lidar_be_tcurvature",
               smooth=0.5,
               tension=40,
               overwrite=True,
               nprocs=24
            )

Visualize the DTM

Interactive Map (Folium)

# Create the shaded relief map
gs.run_command(
    "r.relief",
    input="lidar_be",
    output="hillshade",
    zscale=1,
    overwrite=True,
)

m = gj.InteractiveMap(width="500", tiles="OpenStreetMap", map_backend="folium")
m.add_raster("hillshade", opacity=0.75)
m.add_raster("lidar_be", opacity=0.5)
m.show()

Interactive Map (ipyleaflet)

m = gj.InteractiveMap(width="500", map_backend="ipyleaflet")
m.query_mode = "true"
m.add_raster("hillshade", opacity=0.85)
m.add_raster("points_median", opacity=0.5)
m.add_layer_control()
m.show()

DTM, Slope, and Aspect Figures

DTM 1m Resolution
m = gj.Map()
m.d_shade(color="lidar_be", shade="hillshade")
m.d_legend(raster="lidar_be", at=(5, 9, 50, 90), flags="b", unit="m")
m.d_barscale(at=(5, 6), flags="n")
m.show()

Aspect

gs.run_command("r.colors", map="lidar_be_aspect", color="aspect")
m = gj.Map()
m.d_rast(map="lidar_be_aspect")
m.d_legend(raster="lidar_be_aspect", at=(5, 9, 50, 90), flags="b")
m.d_barscale(at=(5, 7), flags="n")
m.show()

Slope

gs.run_command("r.colors", map="lidar_be_slope", color="sepia", flags="e")
m = gj.Map()
m.d_shade(color="lidar_be_slope", shade="hillshade")
m.d_legend(raster="lidar_be_slope", at=(5, 9, 50, 90), flags="bd")
m.d_barscale(at=(5, 7), flags="n")
m.show()

Calculate Geomorphons

gs.run_command(
        "r.geomorphon",
        elevation="lidar_be",
        forms="geomorphon",
        search=21,
        overwrite=True,
    )

m = gj.Map()
m.d_shade(color="geomorphon", shade="hillshade")
m.d_legend(raster="geomorphon", at=(60, 95, 85, 90), flags="bd")
m.d_barscale(at=(5, 7), flags="n")
m.show()

3D Visualization

elevation_3dmap = gj.Map3D()
# Full list of options m.nviz.image
# https://grass.osgeo.org/grass84/manuals/m.nviz.image.html
elevation_3dmap.render(
    elevation_map="lidar_be",
    color_map="geomorphon",
    zexag=1, 
    perspective=20, 
    height=4000, 
    resolution_fine=1,
    fringe=['nw','ne','sw','se'],
    fringe_elevation=1000,
    arrow_position=[100,50],
)
elevation_3dmap.show()

Stream Extraction

gs.run_command("r.stream.extract", elevation="lidar_be", threshold=500,
                mexp=0.5, stream_length=500, memory=100000, stream_raster="stream_r",
                direction="direction_r", stream_vector="stream_vect")

Let’s view the stream network on the map.

m = gj.InteractiveMap(width="500", map_backend="folium")
m.query_mode = "true"
m.add_raster("hillshade", opacity=0.85)
m.add_raster("points_median", opacity=0.5)
m.add_vector("stream_vect", color="blue", weight=3, type="line")
m.add_layer_control()
m.show()

Overland Flow Simulation

gs.run_command(
    "r.slope.aspect",
    elevation="lidar_be",
    dx="dx",
    dy="dy",
    nprocs=4,
    overwrite=True,
)

OUTPUT_STEP = 2

# Run the simulation
gs.run_command(
        "r.sim.water",
        elevation="lidar_be",
        dx="dx",
        dy="dy",
        rain_value=50,  # mm/hr
        infil_value=0.0,  # mm/hr
        man_value=0.1,
        niterations=10,  # event duration (minutes)
        output_step=OUTPUT_STEP,  # minutes
        depth="depth",  # m
        discharge="disch",  # m3/s
        random_seed=3,
        nprocs=4,
        flags="t",
        overwrite=True,
    )

# Register the output maps into a space time dataset
gs.run_command(
    "t.create",
    output="depth_sum",
    type="strds",
    temporaltype="absolute",
    title="Runoff Depth",
    description="Runoff Depth in [m]",
    overwrite=True,
)

# Get the list of naip maps (bands)
depth_list = gs.read_command(
    "g.list", type="raster", pattern="depth.*", separator="comma"
).strip()

# Register the maps
gs.run_command(
    "t.register",
    input="depth_sum",
    type="raster",
    start="2024-01-01",
    increment=f"{OUTPUT_STEP} minutes",
    maps=depth_list,
    flags="i",
    overwrite=True,
)

Visualize the results of the overland flow simulation.

m = gj.Map()
m.d_shade(color="depth.10", shade="hillshade")
m.d_legend(raster="depth.10", at=(5, 9, 50, 90), flags="bd", unit="m")
m.d_barscale(at=(5, 7), flags="n")
m.show()

Create Animation from Time Series Data

depth_sum_ts_map = gj.TimeSeriesMap(height=600, width=600, use_region=True)
depth_sum_ts_map.add_raster_series("depth_sum")
depth_sum_ts_map.d_legend()
depth_sum_ts_map.render()
depth_sum_ts_map.save(f"outputs/depth.gif")
depth_sum_ts_map.show()

Export GIF

Animation

Select a sample point for analysis

m = gj.InteractiveMap(width=500, map_backend="ipyleaflet")
m.query_mode = "true"
m.add_raster("hillshade", opacity=0.85)
m.add_raster("depth.08", opacity=0.5)
m.add_raster("depth.10", opacity=0.5)
m.add_layer_control()
m.show()

Using our sample point, we can extract the time series data for the point.

gs.parse_command('v.what.strds', input='sample', strds='depth_sum', output="depth_sum_sample")
sample_json = gs.parse_command('v.db.select', map="depth_sum_sample", format="json")

Now we can visualize the time series data for our sample point.

# This is ugly and need to be improved
records = sample_json['records'][0]
del records['cat']
sample_df = pd.DataFrame([records])
transposed = sample_df.T
transposed.columns = ['depth']
sequence = [2, 4, 6, 8, 10]
transposed['time'] = sequence
transposed.set_index('time', inplace=True)
transposed.head()
depth
time
2 0.023435
4 0.027074
6 0.029835
8 0.049045
10 0.317990

Create a Time Series Line Plot

# Create a seaborn line plot
import matplotlib.pyplot as plt
transposed.dropna(inplace=True)

plt.figure(figsize=(10, 6))
transposed.plot(kind='line', y='depth', color='blue', marker='o', linestyle='-', linewidth=2)

# Add titles and labels
plt.title('Sample Depth vs. Time')
plt.xlabel('Time (minutes)')
plt.ylabel('Depth (m)')


# Save the plot as an image file
plt.savefig('outputs/line_plot.png')

# Show the plot
plt.show()
<Figure size 1000x600 with 0 Axes>

Add the sample point to the map

depth_json = gs.parse_command("t.rast.list", input="depth_sum", format="json")

for i in depth_json['data']:
    map_name = i['name']
    file_name = f"outputs/depth_{map_name.split('.')[-1]}.png"
    print(f"Creating {file_name}")
    i['thumbnail'] = file_name
    depth_map = gj.Map(filename=file_name, use_region=True)
    depth_map.d_shade(color=map_name, shade="hillshade")
    depth_map.d_vect(map="depth_sum_sample", color="red", fill_color="red", type="point", size=12, icon="basic/point")
    depth_map.d_legend(raster=map_name, at=(5, 9, 50, 90), flags="bd", unit="m")
    depth_map.d_barscale(at=(5, 7), flags="n")
    depth_map.show()
Creating outputs/depth_02.png
Creating outputs/depth_04.png
Creating outputs/depth_06.png
Creating outputs/depth_08.png
Creating outputs/depth_10.png

Final figure

Create our final figure

from PIL import Image

fig = plt.figure(figsize=(25, 30))
ax = fig.add_subplot(2, 2, 1)
ax.set_axis_off()

fig.subplots_adjust(hspace=0, wspace=0.1)

img = Image.open(depth_map.filename) 
plt.imshow(img)
ax.set_title("Sampe Point - 10 mintues", {"fontsize": 24, "fontweight": "bold"})

ax = fig.add_subplot(2, 2, 2)
ax.set_axis_off()
img = Image.open("outputs/line_plot.png")
plt.imshow(img)

plt.tight_layout()
plt.savefig("outputs/figure.png", bbox_inches="tight", dpi=300)

plt.show()

STAC Integration

NAIP Data

!t.stac.item url="https://planetarycomputer.microsoft.com/api/stac/v1" collection="naip" datetime="2022" format="plain" method="nearest" extent=region nprocs=2 -d
Setting bbox to current region:
BBOX: [-108.00541249, 36.07862566, -107.9985924, 36.08371681]
/home/coreywhite/Documents/GitHub/ncsu-geoforall-lab/tutorials/venv/lib/python3.10/site-packages/pystac_client/item_search.py:678: UserWarning: numberMatched or context.matched not in response
  warnings.warn("numberMatched or context.matched not in response")
Search Matched: None items
2 Assets Ready for download...
Downloading assets:   0%|                                 | 0/2 [00:00<?, ?it/s]Downloading Asset: {'href':
'https://naipeuwest.blob.core.windows.net/naip/v002/nm/2022/nm_060cm_2022/36108/m_3610864_ne_12_060_20220525.tif',
'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
'title': 'RGBIR COG tile', 'eo:bands': [{'name': 'Red', 'common_name':
'red'}, {'name': 'Green', 'common_name': 'green'}, {'name': 'Blue',
'common_name': 'blue'}, {'name': 'NIR', 'common_name': 'nir',
'description': 'near-infrared'}], 'roles': ['data'], 'collection_id':
'naip', 'item_id': 'nm_m_3610864_ne_12_060_20220525', 'file_name':
'naip.nm_m_3610864_ne_12_060_20220525.image', 'datetime':
'2022-05-25T16:00:00Z'}
Downloading Asset: {'href':
'https://naipeuwest.blob.core.windows.net/naip/v002/nm/2022/nm_060cm_2022/36107/m_3610757_nw_13_060_20220525.tif',
'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
'title': 'RGBIR COG tile', 'eo:bands': [{'name': 'Red', 'common_name':
'red'}, {'name': 'Green', 'common_name': 'green'}, {'name': 'Blue',
'common_name': 'blue'}, {'name': 'NIR', 'common_name': 'nir',
'description': 'near-infrared'}], 'roles': ['data'], 'collection_id':
'naip', 'item_id': 'nm_m_3610757_nw_13_060_20220525', 'file_name':
'naip.nm_m_3610757_nw_13_060_20220525.image', 'datetime':
'2022-05-25T16:00:00Z'}
Import Url:
/vsicurl/https://naipeuwest.blob.core.windows.net/naip/v002/nm/2022/nm_060cm_2022/36108/m_3610864_ne_12_060_20220525.tif
Import Url:
/vsicurl/https://naipeuwest.blob.core.windows.net/naip/v002/nm/2022/nm_060cm_2022/36107/m_3610757_nw_13_060_20220525.tif
Importing: naip.nm_m_3610864_ne_12_060_20220525.image
Importing: naip.nm_m_3610757_nw_13_060_20220525.image
Downloading assets: 100%|█████████████████████████| 2/2 [00:25<00:00, 12.90s/it]

We can view raw data with GRASS using the d_rgb tool.

m = gj.Map()
m.d_rgb(red="naip.nm_m_3610757_nw_13_060_20220525.image.1",
        blue="naip.nm_m_3610757_nw_13_060_20220525.image.3",
        green="naip.nm_m_3610757_nw_13_060_20220525.image.2")
m.d_rgb(red="naip.nm_m_3610864_ne_12_060_20220525.image.1",
        blue="naip.nm_m_3610864_ne_12_060_20220525.image.3",
        green="naip.nm_m_3610864_ne_12_060_20220525.image.2")
m.d_vect(map="stream_vect", color="blue", type="line")

m.show()

Mosaic and Compiste

Let create a mosaic of the NAIP data, and give the bands some names.

def patch_and_composite_naip(year=2022):
    gs.run_command("g.region", res=1)
    naip_bands = [(1, "red"), (2, "green"), (3, "blue"), (4, "nir")]
    for band in naip_bands:
        i, band_name = band
        # Get the list of depth maps
        image_list = gs.read_command(
            "g.list", type="raster", pattern=f"*.image.{i}", separator="comma"
        ).strip()

        gs.run_command(
            "r.patch",
            input=image_list,
            output=f"naip_{year}.{band_name}",
            nprocs=4,
            memory=2100,
            overwrite=True,
        )

    gs.run_command(
        "r.composite",
        red=f"naip_{year}.red",
        green=f"naip_{year}.green",
        blue=f"naip_{year}.blue",
        output=f"naip_{year}_rgb",
        overwrite=True,
    )

patch_and_composite_naip()

Now we can view the naip compsite.

image_map = gj.Map(filename="outputs/naip_2022_rgb.png", use_region=True)
image_map.d_shade(color="naip_2022_rgb", shade="hillshade")
image_map.d_grid(siz=100, color="black", flags="a")
image_map.d_vect(map="depth_sum_sample", color="red", fill_color="red", type="point", size=12, icon="basic/point")
image_map.d_barscale(at=(5, 7), flags="n")
image_map.show()

NDVI

We can use the NAIP data to calculate NDVI.

gs.run_command("i.vi", red="naip_2022.red", nir="naip_2022.nir", output="ndvi", overwrite=True)
ndvi_map = gj.Map(filename="outputs/ndvi.png", height=600, width=600)
ndvi_map.d_shade(color="ndvi", shade="hillshade")
ndvi_map.d_grid(siz=100, color="grey", flags="a")
ndvi_map.d_legend(raster="ndvi", at=(10, 12, 50, 90), flags="bd")
ndvi_map.d_barscale(at=(5, 11), flags="n")
ndvi_map.show()

Let’s create our final figure.

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from PIL import Image

# Create a figure
fig = plt.figure(figsize=(15, 10))
fig.suptitle('Chaco Canyon 2016 - Overland Flow Simulation', fontsize=24, fontweight='bold')
# Create a GridSpec with 2 rows and 6 columns, making the bottom row larger
grid = gridspec.GridSpec(2, 6, height_ratios=[2, 1])

# Add an image to the left of the plot on the top row
ax_left = fig.add_subplot(grid[0, 1:3])
ax_left.set_title('NAIP 2022', {"fontsize": 18, "fontweight": "bold"})
img_left = Image.open(image_map.filename)
ax_left.imshow(img_left)
ax_left.axis("off")

# Add a larger line graph on the top row spanning the remaining columns
ax_large = fig.add_subplot(grid[0, 3:])
img_large = Image.open("outputs/line_plot.png")
ax_large.imshow(img_large)
# ax_large.set_title('Depth vs. Time', {"fontsize": 18, "fontweight": "bold"})
ax_large.axis("off")

ax = fig.add_subplot(grid[1, 3:4])
ax.set_title(f"Depth Maps", {"fontsize": 16, "fontweight": "bold"})
ax.axis("off")
# Add 5 subplots on the bottom row
for idx, item in enumerate(depth_json['data']):
    ax = fig.add_subplot(grid[1, idx+1])
    map_name = item['name']
    time_stamp = map_name.split('.')[-1]
    img = Image.open(item['thumbnail'])  # Open the image file
    ax.imshow(img)
    ax.set_title(f"{time_stamp} minutes", {"fontsize": 16, "fontweight": "bold"})
    ax.axis("off")

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

Funding

Pathways to Enable Open-Source Ecosystems (POSE) Phase II

NSF Award # 2303651