geeViz Report Generation Examples

This notebook demonstrates the geeViz.outputLib.reports module for generating automated, publication-ready reports from Google Earth Engine data. Each example shows a realistic user scenario with different study areas, datasets, and chart types.

Key concepts:

• report = rl.Report(title, theme, tone) — create a report

• report.add_section(ee_obj, geometry, chart_types=[...], ...) — add data sections

• report.generate(format="html") — produce the output

• chart_types accepts a list of 0–3 chart types per section (e.g. ["sankey", "line+markers"])

• Study areas come from sal.* functions (counties, forests, protected areas, etc.), but can be any EE Geometry, Feature, or FeatureCollection

Requirements:

• geeViz with outputLib

• GOOGLE_API_KEY in environment for LLM narratives (optional — reports work without it)

#	Example	Study Area	Datasets	chart_types	thumb_format
1	Single-section basics	County (sal.getUSCounties)	LCMS Land Cover 2024	auto (bar)	png
2	Fire analysis	EDW fire perimeter	MTBS + Landsat + LCMS sankey	sankey, line+markers	png, gif
3	National Forest comparison	USFS Forests (sal.getUSFSForests)	LCMS Land Cover	stacked_line+markers	png
4	NLCD urbanization	Counties (sal.getUSCounties)	Annual NLCD	sankey, stacked_line+markers	gif
5	Vegetation trends	Protected area (sal.getProtectedAreas)	Landsat NDVI/NBR	line+markers	none
6	Multi-chart snapshot	State (sal.getUSStates)	LCMS Land Use 2024	bar, donut	png

# Setup — shared by all examples
import os, time
import geeViz.geeView as gv
import geeViz.getImagesLib as gil
import geeViz.getSummaryAreasLib as sal
import geeViz.edwLib as edw
from geeViz.outputLib import reports as rl
from IPython.display import display, IFrame, HTML
ee = gv.ee

output_dir = os.path.join(os.path.dirname(os.getcwd()), "examples", "outputs", "reports")
os.makedirs(output_dir, exist_ok=True)

def show_report(path, height=800):
    """Display a generated HTML/PDF report inline in the notebook."""
    if path.endswith(".pdf"):
        rel = os.path.relpath(path, os.getcwd())
        display(HTML(f'<object data="{rel}" type="application/pdf" width="100%" height="{height}px">'
                     f'<p>PDF cannot be displayed inline. <a href="{rel}">Download here</a>.</p></object>'))
    else:
        rel = os.path.relpath(path, os.getcwd())
        display(IFrame(src=rel, width="100%", height=height))

print(f"Reports will be saved to: {output_dir}")
print("Use show_report(path) after generate() to display the report inline.")

---

Example 1: Single-Section Report (Beginner)

User question: “What does the current land cover look like in Salt Lake County, Utah?”

This is the simplest possible report — one dataset, one study area, one section. The chart type is auto-detected (bar chart for a single Image).

MCP equivalent

If using the geeViz MCP server, you would call:

create_report(title="Salt Lake County Land Cover", theme="dark", tone="informative")

Then in run_code:

study_area = sal.getUSCounties(ee.Geometry.Point([-111.89, 40.76]))
lcms = ee.ImageCollection("USFS/GTAC/LCMS/v2024-10")
lcms_2024 = lcms.filter(ee.Filter.calendarRange(2024, 2024, 'year')).filterBounds(study_area).mosaic()
lcms_2024 = lcms_2024.copyProperties(lcms.first())

Then:

add_report_section(
    ee_obj_var="lcms_2024", geometry_var="study_area",
    title="LCMS Land Cover 2024",
    band_names="Land_Cover", scale=60,
    basemap="esri-satellite", burn_in_geometry=True
)
generate_report(format="html", output_filename="ex1_salt_lake_county.html")

# Example 1: Single-section report — LCMS Land Cover in Salt Lake County
study_area_1 = sal.getUSCounties(ee.Geometry.Point([-111.89, 40.76]))
print("County:", study_area_1.aggregate_array("NAMELSAD").getInfo())

lcms = ee.ImageCollection("USFS/GTAC/LCMS/v2024-10")
lcms_2024 = lcms.filter(ee.Filter.calendarRange(2024, 2024, "year")) \
    .filterBounds(study_area_1).mosaic().copyProperties(lcms.first())

report1 = rl.Report(
    title="Salt Lake County Land Cover",
    theme="dark",
    tone="informative",
    header_text="Current land cover from USFS LCMS v2024-10.",
)
report1.add_section(
    ee_obj=lcms_2024,
    geometry=study_area_1,
    title="LCMS Land Cover (2024)",
    band_names=["Land_Cover"],
    scale=60,
    basemap="esri-satellite",
    burn_in_geometry=True,
    geometry_outline_color="white",
)

out_path = os.path.join(output_dir, "ex1_salt_lake_county.html")
report1.generate(format="html", output_path=out_path)
print(f"\nReport saved to: {out_path}")
print(report1.metadata().to_markdown(index=False))
show_report(out_path)

---

Example 2: Multi-Section Fire Analysis

User question: “Analyze the 2020 Cameron Peak Fire in Colorado — show the burn severity, before/after imagery, and how land cover transitioned.”

This example uses:

• edw.query_features() to get the actual fire perimeter from MTBS

• MTBS burn severity mosaic

• Landsat composites (pre-fire 2019, post-fire 2021, recent 2024)

• LCMS Land Cover sankey transitions (2019 → 2021 → 2024)

MCP equivalent

# Step 1: Find the fire perimeter
search_edw(query="MTBS")
get_edw_service_info(service_name="EDW_MTBS_01")
query_edw_features(
    service_name="EDW_MTBS_01", layer_id=75,
    where="FIRE_NAME LIKE '%CAMERON PEAK%'",
    out_fields="FIRE_NAME,ACRES,YEAR"
)

# Step 2: In run_code, load the GeoJSON into EE
geojson = edw.query_features("EDW_MTBS_01", 75,
    where="FIRE_NAME LIKE '%CAMERON PEAK%'",
    out_fields="FIRE_NAME,ACRES,YEAR")
fire_fc = ee.FeatureCollection(geojson)

# Step 3: Build the report
create_report(title="Cameron Peak Fire Analysis", theme="dark")
add_report_section(ee_obj_var="mtbs", geometry_var="fire_geom", title="Burn Severity", ...)
add_report_section(
    ee_obj_var="lcms_lc", geometry_var="fire_geom",
    title="Land Cover Transitions",
    chart_types="sankey,line+markers",
    transition_periods="[[2019,2019],[2021,2021],[2024,2024]]",
    sankey_band_name="Land_Cover"
)
generate_report(format="html")

# Example 2: Get the Cameron Peak Fire perimeter from EDW using edwLib
# Layer 75 = "2020 Burned Area Boundaries" in EDW_MTBS_01

geojson = edw.query_features(
    "EDW_MTBS_01", 75,
    where="FIRE_NAME LIKE '%CAMERON PEAK%'",
    out_fields="FIRE_NAME,ACRES,YEAR",
)

for f in geojson["features"]:
    p = f["properties"]
    print(f"  {p['FIRE_NAME']} ({p['YEAR']}) — {p['ACRES']:,.0f} acres")

fire_fc = ee.FeatureCollection(geojson)
fire_geom = fire_fc.geometry()
study_area_2 = fire_geom.bounds().buffer(2000)
print(f"\nFire perimeter loaded into EE.")

# Set up all datasets for Example 2

# MTBS burn severity
mtbs_2020 = ee.ImageCollection("USFS/GTAC/MTBS/annual_burn_severity_mosaics/v1") \
    .filterDate("2020-01-01", "2020-12-31") \
    .filterBounds(study_area_2).first().select([0], ["Severity"])

# Landsat composites: pre-fire, post-fire, recent
def get_landsat_composite(area, year):
    return gil.getLandsatWrapper(
        area, year, year, startJulian=152, endJulian=273
    )["processedComposites"].first().select(["swir2", "nir", "red"])

pre_img = get_landsat_composite(study_area_2, 2019)
post_img = get_landsat_composite(study_area_2, 2021)
recent_img = get_landsat_composite(study_area_2, 2024)

# LCMS for transitions
lcms_lc = ee.ImageCollection("USFS/GTAC/LCMS/v2024-10") \
    .filterBounds(study_area_2).select(["Land_Cover"])

print("All datasets ready.")

# Build the fire report
fire_thumb = dict(
    basemap="esri-satellite",
    burn_in_geometry=True,
    geometry_outline_color="red",
    geometry_fill_color="FF000020",
)
landsat_viz = {"min": 0.05, "max": 0.45, "gamma": 1.4}

report2 = rl.Report(
    title="Cameron Peak Fire Analysis (2020)",
    theme="dark",
    tone="informative",
    header_text="The Cameron Peak Fire burned over 208,000 acres in Larimer County, Colorado in 2020. "
               "Fire perimeter from USFS EDW/MTBS.",
)

# Section 1: Burn severity
report2.add_section(
    ee_obj=mtbs_2020, geometry=fire_geom,
    title="MTBS Burn Severity (2020)",
    band_names=["Severity"],
    prompt="Describe the burn severity distribution.",
    **fire_thumb,
)

# Sections 2-4: Landsat before / after / recent (thumbnail only)
for img, yr, label in [
    (pre_img, 2019, "Pre-Fire"),
    (post_img, 2021, "Post-Fire"),
    (recent_img, 2024, "Recovery"),
]:
    report2.add_section(
        ee_obj=img, geometry=fire_geom,
        title=f"{label} Landsat (Summer {yr})",
        thumb_viz_params=landsat_viz,
        prompt=f"Describe the {label.lower()} landscape.",
        generate_chart=False, generate_table=False,
        **fire_thumb,
    )

# Section 5: Land Cover sankey + time series (both charts)
report2.add_section(
    ee_obj=lcms_lc, geometry=fire_geom,
    title="Land Cover Transitions (2019 → 2021 → 2024)",
    chart_types=["sankey", "line+markers"],
    transition_periods=[[2019, 2019], [2021, 2021], [2024, 2024]],
    sankey_band_name="Land_Cover",
    scale=30,
    thumb_format=None,
    prompt="Describe land cover transitions showing fire impact and early recovery.",
)

out_path = os.path.join(output_dir, "ex2_cameron_peak_fire.html")
t0 = time.time()
report2.generate(format="html", output_path=out_path)
print(f"\nGenerated in {time.time()-t0:.0f}s: {out_path}")
print(report2.metadata().to_markdown(index=False))
show_report(out_path)

---

Example 3: Multi-Feature Comparison Across National Forests

User question: “How does vegetation compare across the National Forests of the Northern Rockies?”

This uses sal.getUSFSForests() with a region filter to get multiple forest boundaries, then produces per-feature time series subplots using feature_label.

MCP equivalent

# In run_code:
forests = sal.getUSFSForests(ee.Geometry.Point([-114, 46.5]).buffer(300000), region='01')
lcms_lc = ee.ImageCollection('USFS/GTAC/LCMS/v2024-10').select(['Land_Cover'])

add_report_section(
    ee_obj_var="lcms_lc", geometry_var="forests",
    title="Land Cover by National Forest",
    feature_label="FORESTNAME", scale=120,
    chart_types="stacked_line+markers",
    basemap="esri-satellite", burn_in_geometry=True
)

# Example 3: Multi-feature comparison — National Forests in the Northern Rockies
forests = sal.getUSFSForests(
    ee.Geometry.Point([-114, 46.5]).buffer(300000),
    region="01",
)
forest_names = forests.aggregate_array("FORESTNAME").distinct().getInfo()
print(f"Forests found ({len(forest_names)}): {forest_names}")

lcms_lc_3 = ee.ImageCollection("USFS/GTAC/LCMS/v2024-10") \
    .select(["Land_Cover"]).filterBounds(forests)

report3 = rl.Report(
    title="Northern Rockies Forest Comparison",
    theme="dark",
    tone="technical",
    header_text="Comparing LCMS land cover across National Forests in USFS Region 1.",
)

# Single section with per-feature subplots
report3.add_section(
    ee_obj=lcms_lc_3,
    geometry=forests,
    title="Land Cover by National Forest (1985–2024)",
    feature_label="FORESTNAME",
    chart_types=["stacked_line+markers"],
    scale=120,
    basemap="esri-satellite",
    burn_in_geometry=True,
    geometry_outline_color="white",
    prompt="Compare the long-term land cover trends across these National Forests.",
)

out_path = os.path.join(output_dir, "ex3_northern_rockies_forests.html")
t0 = time.time()
report3.generate(format="html", output_path=out_path)
print(f"\nGenerated in {time.time()-t0:.0f}s: {out_path}")
print(report3.metadata().to_markdown(index=False))
show_report(out_path)

---

Example 4: NLCD Urbanization Analysis with Sankey

User question: “How has urbanization changed land cover around Denver from 1985 to 2024?”

This demonstrates:

• Annual NLCD (not official NLCD releases) — requires manual band rename and class properties

• chart_types=["sankey", "stacked_line+markers"] — both charts from one section

• thumb_format="gif" — animated time-lapse

Important: Annual NLCD setup

Annual NLCD bands are named b1 with no class properties. You must rename and set properties:

nlcd_lc = ee.ImageCollection('projects/sat-io/open-datasets/USGS/ANNUAL_NLCD/LANDCOVER')
nlcd_lc = nlcd_lc.map(lambda img: img.rename('LC').set(nlcd_viz_props))

Do NOT remap values to sequential integers (0, 1, 2…) — the *_class_values property maps originals (11, 12, 21…) to colors.

# Example 4: NLCD urbanization — Denver metro area
denver_counties = sal.getUSCounties(
    ee.Geometry.Point([-104.99, 39.74]).buffer(50000),
    state_abbr="CO",
)
print("Counties:", denver_counties.aggregate_array("NAMELSAD").getInfo())

# Set up Annual NLCD with proper class properties
nlcd_lc = ee.ImageCollection("projects/sat-io/open-datasets/USGS/ANNUAL_NLCD/LANDCOVER")
nlcd_viz_props = {
    "LC_class_values": [11,12,21,22,23,24,31,41,42,43,52,71,81,82,90,95],
    "LC_class_palette": [
        "466b9f","d1def8","dec5c5","d99282","eb0000","ab0000",
        "b3ac9f","68ab5f","1c5f2c","b5c58f","ccb879","dfdfc2",
        "dcd939","ab6c28","b8d9eb","6c9fb8",
    ],
    "LC_class_names": [
        "Open Water","Perennial Ice/Snow","Developed, Open Space",
        "Developed, Low Intensity","Developed, Medium Intensity",
        "Developed, High Intensity","Barren Land","Deciduous Forest",
        "Evergreen Forest","Mixed Forest","Shrub/Scrub",
        "Grassland/Herbaceous","Pasture/Hay","Cultivated Crops",
        "Woody Wetlands","Emergent Herbaceous Wetlands",
    ],
}
nlcd_lc = nlcd_lc.map(lambda img: img.rename("LC").set(nlcd_viz_props))

report4 = rl.Report(
    title="Denver Metro Urbanization (1985–2024)",
    theme="dark",
    tone="informative",
    header_text="Tracking urban growth using Annual NLCD land cover data.",
)

report4.add_section(
    ee_obj=nlcd_lc,
    geometry=denver_counties,
    title="NLCD Land Cover Change (1985–2024)",
    chart_types=["sankey", "stacked_line+markers"],
    transition_periods=[[1985, 1985], [2000, 2000], [2010, 2010], [2024, 2024]],
    sankey_band_name="LC",
    scale=30,
    thumb_format="gif",
    basemap="esri-satellite",
    burn_in_geometry=True,
    geometry_outline_color="cyan",
    prompt="Describe the urbanization trends in the Denver metro area over four decades.",
)

out_path = os.path.join(output_dir, "ex4_denver_nlcd.html")
t0 = time.time()
report4.generate(format="html", output_path=out_path)
print(f"\nGenerated in {time.time()-t0:.0f}s: {out_path}")
print(report4.metadata().to_markdown(index=False))
show_report(out_path)

---

Example 5: Continuous Vegetation Trends (NDVI/NBR)

User question: “How has vegetation greenness trended in Glacier National Park over the last 25 years?”

This demonstrates:

• Protected areas via sal.getProtectedAreas(iucn_cat="II") (IUCN Category II = National Park)

• Continuous (non-thematic) data — Landsat NDVI and NBR indices

• getLandsatWrapper with startJulian/endJulian for summer composites

• chart_types=["line+markers"] — simple time series

MCP equivalent

# In run_code:
pa = sal.getProtectedAreas(ee.Geometry.Point([-113.8, 48.7]), iucn_cat='II')
study_area = pa.filter(ee.Filter.stringContains('NAME', 'Glacier')).first().geometry()
composites = gil.getLandsatWrapper(study_area, 2000, 2024, startJulian=152, endJulian=273)['processedComposites']

add_report_section(
    ee_obj_var="composites", geometry_var="study_area",
    title="NDVI & NBR Trends", band_names="NDVI,NBR",
    chart_types="line+markers", thumb_format="none"
)

# Example 5: Continuous trends — Glacier National Park NDVI/NBR
pa = sal.getProtectedAreas(
    ee.Geometry.Point([-113.8, 48.7]),
    iucn_cat="II",
)
glacier = pa.filter(ee.Filter.stringContains("NAME", "Glacier")).first().geometry()
print("Study area: Glacier National Park")

composites_5 = gil.getLandsatWrapper(
    glacier, 2000, 2024,
    startJulian=152, endJulian=273,
)["processedComposites"]

report5 = rl.Report(
    title="Glacier National Park Vegetation Trends",
    theme="dark",
    tone="technical",
    header_text="Summer Landsat spectral index trends from 2000–2024.",
)

report5.add_section(
    ee_obj=composites_5,
    geometry=glacier,
    title="Summer NDVI & NBR (2000–2024)",
    band_names=["NDVI", "NBR"],
    chart_types=["line+markers"],
    scale=60,
    thumb_format=None,
    prompt="Describe the vegetation greenness and burn recovery trends.",
)

out_path = os.path.join(output_dir, "ex5_glacier_ndvi.html")
t0 = time.time()
report5.generate(format="html", output_path=out_path)
print(f"\nGenerated in {time.time()-t0:.0f}s: {out_path}")
print(report5.metadata().to_markdown(index=False))
show_report(out_path)

---

Example 6: Bar + Donut Multi-Chart Snapshot

User question: “What is the current land use breakdown for Vermont?”

This demonstrates:

• US States via sal.getUSStates()

• LCMS Land Use (not Land Cover) — a different thematic product

• chart_types=["bar", "donut"] — two chart types from the same data

• Single ee.Image (mosaic of 2024 tiles)

MCP equivalent

create_report(title="Vermont Land Use", theme="dark")
add_report_section(
    ee_obj_var="lcms_lu_2024", geometry_var="vermont",
    title="Land Use 2024", band_names="Land_Use",
    chart_types="bar,donut",
    basemap="esri-topo", burn_in_geometry=True
)
generate_report(format="html")

# Example 6: Bar + Donut — Vermont Land Use
vermont = sal.getUSStates(ee.Geometry.Point([-72.58, 44.26]))
print("State:", vermont.aggregate_array("NAME").getInfo())

lcms_lu_2024 = ee.ImageCollection("USFS/GTAC/LCMS/v2024-10") \
    .filter(ee.Filter.calendarRange(2024, 2024, "year")) \
    .filterBounds(vermont).select(["Land_Use"]) \
    .mosaic().copyProperties(
        ee.ImageCollection("USFS/GTAC/LCMS/v2024-10").first()
    )

report6 = rl.Report(
    title="Vermont Land Use (2024)",
    theme="dark",
    tone="informative",
    header_text="Current land use classification from USFS LCMS.",
)

report6.add_section(
    ee_obj=lcms_lu_2024,
    geometry=vermont,
    title="LCMS Land Use (2024)",
    band_names=["Land_Use"],
    chart_types=["bar", "donut"],
    scale=60,
    basemap="esri-topo",
    burn_in_geometry=True,
    geometry_outline_color="white",
    prompt="Describe the land use distribution across Vermont.",
)

out_path = os.path.join(output_dir, "ex6_vermont_land_use.html")
t0 = time.time()
report6.generate(format="html", output_path=out_path)
print(f"\nGenerated in {time.time()-t0:.0f}s: {out_path}")
print(report6.metadata().to_markdown(index=False))
show_report(out_path)

---

Summary

Example	Pattern	Key Features
1	Single section, auto-detect	sal.getUSCounties, single Image, bar chart auto-detected
2	Multi-section fire report	EDW fire perimeter, MTBS + Landsat + LCMS, chart_types=["sankey", "line+markers"]
3	Multi-feature comparison	sal.getUSFSForests, feature_label, per-feature stacked subplots
4	NLCD with sankey	Annual NLCD (manual setup), chart_types=["sankey", "stacked_line+markers"], GIF
5	Continuous trends	sal.getProtectedAreas, Landsat NDVI/NBR, line+markers
6	Multi-chart snapshot	sal.getUSStates, LCMS Land Use, chart_types=["bar", "donut"]

Tips

• chart_types list controls which charts per section (0–3). Auto-detected when empty.

• thumb_format: "png" (static), "gif" (animated), "filmstrip" (grid), None (no map)

• feature_label + ee.FeatureCollection → per-feature subplots (IC) or grouped bar (Image)

• MTBS: Always .select([0], ["Severity"]) for consistent band naming

• Annual NLCD: Must rename band and set class properties manually

• Landsat: getLandsatWrapper requires startJulian/endJulian (1–365)

• Tiled collections (LCMS, NLCD): Never use .first() — pass full ImageCollection