Summary Areas, Thumbnails, and Charting with geeViz

• geeViz.getSummaryAreasLib provides functions to retrieve common summary and study area ee.FeatureCollection objects

• geeViz.thumbLib generates automatic thumbnails from Earth Engine data with auto-visualization detection

• geeViz.chartingLib summarizes and charts EE data with automatic thematic/continuous detection

• geeViz.reportLib combines all three into styled HTML reports with LLM-generated narratives

• Every function accepts an area parameter (ee.Geometry, ee.Feature, or ee.FeatureCollection) and returns spatially filtered results

• This notebook demonstrates each module’s capabilities

• Datasets used: LCMS, Annual NLCD, MTBS, and Sentinel-2 composites via superSimpleGetS2

Copyright 2026 Ian Housman

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

# Boilerplate imports
try:
    import geeViz.geeView as geeView
except:
    !python -m pip install geeViz
    import geeViz.geeView as geeView

import geeViz.getImagesLib as gil
from geeViz.outputLib import charts as cl
import geeViz.getSummaryAreasLib as sal
from geeViz.outputLib import thumbs as tl
from geeViz.outputLib import reports as rl

ee = geeView.ee
Map = geeView.Map
Map.clearMap()
print('Imports complete')

from IPython.display import HTML, display

Define a Study Area

• We’ll use a broad area centered on the Wasatch Front in Utah

• This area intersects national forests, urban areas, multiple counties, and protected lands

• All getSummaryAreasLib functions will filter to features intersecting this area

# Study area: point near Salt Lake City with an 80km buffer
study_area = ee.Geometry.Point([-111.9, 40.7]).buffer(80000)
print('Study area defined')

Load Datasets

• We load several datasets that will be summarized across the different summary areas throughout the notebook

• LCMS: Land cover, land use, and change (1985-2024)

• Annual NLCD: Annual National Land Cover Database (1985-2024)

• MTBS: Monitoring Trends in Burn Severity

• Global Land Cover: GLC-FCS30D annual land cover

• Sentinel-2: Cloud-free composites via superSimpleGetS2

# LCMS
lcms = ee.ImageCollection('USFS/GTAC/LCMS/v2024-10')

# Annual NLCD (preferred over official NLCD releases — annual data from 1985-2024)
nlcd = 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 = nlcd.map(lambda img: img.rename('LC').set(nlcd_viz_props))

# MTBS burn severity (CONUS mosaics) - select by index and rename for consistent band naming
mtbs = ee.ImageCollection('USFS/GTAC/MTBS/annual_burn_severity_mosaics/v1').select([0], ['Severity'])

# Sentinel-2 unfiltered collection (will be filtered to study areas later)
s2 = gil.superSimpleGetS2(study_area, '2023-06-01', '2023-09-30')

print('Datasets loaded')

---

1. Countries (getAdminBoundaries(level=0))

• Returns country boundaries from geoBoundaries v6 (default), FAO GAUL, or FieldMaps

• getAdminBoundaries supports levels 0–4: countries, states/provinces, districts/counties, sub-districts, localities

• Use getAdminNameProperty(level, source) to get the name column for any source

countries = sal.getAdminBoundaries(study_area, level=0)
name_prop = sal.getAdminNameProperty(level=0)
print('Countries intersecting study area:')
print(countries.aggregate_array(name_prop).getInfo())

# Also try GAUL source — different name property
countries_gaul = sal.getAdminBoundaries(study_area, level=0, source='gaul')
gaul_name = sal.getAdminNameProperty(level=0, source='gaul')
print(f'Countries (GAUL, prop={gaul_name}):', countries_gaul.aggregate_array(gaul_name).getInfo())

---

2. States/Provinces (getAdminBoundaries(level=1))

• Returns state or province boundaries

• Default source is geoBoundaries v6; also supports GAUL, GAUL 2024, and FieldMaps

• Use getAdminNameProperty to get the correct name column for any source

states = sal.getAdminBoundaries(study_area, level=1)
name_prop = sal.getAdminNameProperty(level=1)
print(f'States/provinces (prop={name_prop}):', states.aggregate_array(name_prop).getInfo())

---

3. US States (getUSStates)

• Returns US state boundaries from TIGER 2018

• Includes NAME, STUSPS (abbreviation), STATEFP (FIPS)

us_states = sal.getUSStates(study_area)
print('US States:', us_states.aggregate_array('NAME').getInfo())

---

4. US Counties (getUSCounties) + LCMS Land Use Charting

• Returns US county boundaries with FULL_NAME, STATEFP, STUSPS

• Optional state_fips or state_abbr filters

• We’ll use counties as summary zones for LCMS Land Use with area charting on the map

• summarize_and_chart with feature_label + ee.ImageCollection produces per-feature time series subplots

• Also shows how to use getAdminBoundaries which is more versatile - works outside the US

# Get Utah counties in the study area
counties = sal.getUSCounties(study_area, state_abbr='UT')
print('Utah counties:', counties.aggregate_array('FULL_NAME').getInfo())

counties2 = sal.getAdminBoundaries(study_area, level = 2)
print('Utan counties using getAdminBoundaries:',counties2.aggregate_array('shapeName').getInfo())

# Map: LCMS Land Use with counties as selectable areas
Map.clearMap()

Map.addLayer(lcms.select(['Land_Use']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True, 'sankey': True}
}, 'LCMS Land Use')

Map.addLayer(lcms.select(['Land_Cover']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True}
}, 'LCMS Land Cover', False)

Map.addSelectLayer(counties, {
    'strokeColor': '00F',
    'selectLayerNameProperty': 'FULL_NAME'
}, 'Utah Counties')

Map.turnOnSelectionAreaCharting()
Map.centerObject(study_area,8)
Map.view()

# Inline: LCMS Land Use time series across counties (per-feature subplots)
# Passing an ImageCollection + feature_label produces one time series subplot per county
# Height is automatically scaled by number of features
top3_counties = counties.sort('Shape_Area', False).limit(3)

per_feature_dfs, fig = cl.summarize_and_chart(
    lcms.select(['Land_Use']).filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    top3_counties,
    feature_label='FULL_NAME',
    title='LCMS Land Use Time Series by County (2000-2024)',
    scale=60,
    width=900,
)
fig.show()

# per_feature_dfs is a dict of {county_name: DataFrame}
for name, feat_df in per_feature_dfs.items():
    print(f'\n--- {name} ---')
    display(feat_df.tail(3))


    

# Sankey: LCMS Land Use transitions for the largest county
sankey_df, sankey_html, matrix_dict = cl.summarize_and_chart(
    lcms.select(['Land_Use']).filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    top3_counties,
    chart_type='sankey',
    transition_periods=[2000, 2010, 2024],
    sankey_band_name='Land_Use',
    title='LCMS Land Use Transitions (2000-2024)',
    min_percentage=0.5,
)
display(HTML(cl.sankey_iframe(sankey_html)))

for k,i in matrix_dict.items():
    print(k)
    display(i)

# Thumbnail: Multi-feature grid in Albers projection with topo basemap
# Note: rotated north arrow and projected extent box in the inset
result = tl.generate_thumbs(
    lcms.select(['Land_Use']), top3_counties,
    dimensions=300, max_features=3, columns=3,
    crs='EPSG:5070',               # Albers Equal Area projection
    basemap='esri-topo',            # topographic basemap
    scalebar_units='imperial',      # show miles instead of km
    inset_map=True,                 # inset should show all 3 features
    inset_basemap='esri-hillshade', # different basemap for inset
    title='LCMS Land Use by County (Albers)',
)
display(HTML(result['html']))

# Same data as 2-column layout with geometry burn-in, no clip, no inset
result2 = tl.generate_thumbs(
    lcms.select(['Land_Use']), top3_counties,
    dimensions=250, max_features=3, columns=2,  # 2 columns = taller grid
    burn_in_geometry=True,
    geometry_outline_color='yellow',
    clip_to_geometry=False,         # show data beyond boundary
    inset_map=False,                # no inset map
)
display(HTML(result2['html']))

---

5. US Urban Areas (getUSUrbanAreas) + Annual NLCD Charting

• Returns TIGER 2024 urban area boundaries

• We’ll summarize Annual NLCD land cover across urban areas

urban = sal.getUSUrbanAreas(study_area)
print('Urban areas:', urban.aggregate_array('NAME20').getInfo())

# Map: NLCD over urban areas
Map.clearMap()

Map.addLayer(nlcd, {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'chartType': 'stacked-bar'}
}, 'Annual NLCD Land Cover')

Map.addSelectLayer(urban, {
    'strokeColor': 'F0F',
    'selectLayerNameProperty': 'NAME20'
}, 'Urban Areas')

Map.turnOnSelectionAreaCharting()
Map.centerObject(study_area)
Map.view()

# Inline: NLCD mode grouped bar across the 3 largest urban areas
top3_urban = urban.sort('ALAND20', False).limit(3)

per_feature_dfs, fig = cl.summarize_and_chart(
    nlcd,
    top3_urban,
    feature_label='NAME20',
    title='Annual NLCD Land Cover Mode - Top 3 Urban Areas',
    scale=60,
    width=800,
    height=400*3,
)
fig.show()
# per_feature_dfs is a dict of {county_name: DataFrame}
for name, feat_df in per_feature_dfs.items():
    print(f'\n--- {name} ---')
    display(feat_df.tail(3))

# Thumbnail: Multi-feature grid with USGS imagery basemap
result = tl.generate_thumbs(
    nlcd, top3_urban,
    dimensions=300, max_features=3,
    basemap='usgs-imagery-topo',    # USGS imagery + topo hybrid
    overlay_opacity=0.7,            # semi-transparent EE overlay
    clip_to_geometry=False,         # show data beyond boundary
    burn_in_geometry=True,
    geometry_outline_color='lime',
)
display(HTML(result['html']))

# Thumbnail from a single ee.Feature (not a FeatureCollection)
single_urban = top3_urban.first()  # ee.Feature
result2 = tl.generate_thumbs(
    nlcd, single_urban,             # ee.Feature works too
    dimensions=400,
    basemap='esri-street',
    inset_map=True, inset_on_map=True,  # inset overlaid on map corner
    title='Largest Urban Area - NLCD',
)
display(HTML(result2['html']))

---

6. Admin Level 2 (getAdminBoundaries(level=2))

• Returns admin-level-2 boundaries (districts/counties) from any supported source

• Works globally — use getAdminNameProperty for the name column

admin2 = sal.getAdminBoundaries(study_area, level=2)
name_prop = sal.getAdminNameProperty(level=2)
print(f'Admin-2 districts: {admin2.size().getInfo()} features')
print(f'Names (prop={name_prop}):', admin2.aggregate_array(name_prop).getInfo())

---

7. USFS National Forests (getUSFSForests) + LCMS Land Cover Charting

• Returns USFS National Forest boundaries

• Properties include FORESTNAME, REGION, GIS_ACRES

• Optional region filter

forests = sal.getUSFSForests(study_area)
print('National Forests:', forests.aggregate_array('FORESTNAME').getInfo())

# Map: LCMS Land Cover over national forests + Sentinel-2 basemap
Map.clearMap()

Map.addLayer(s2.filterBounds(forests.geometry().bounds(500)).median(), gil.vizParamsFalse10k, 'Sentinel-2 Summer 2023', True)

Map.addLayer(lcms.select(['Land_Cover']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True, 'sankey': True}
}, 'LCMS Land Cover')

Map.addLayer(lcms.select(['Change']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True}
}, 'LCMS Change', False)

Map.addSelectLayer(forests, {
    'strokeColor': '0F0',
    'selectLayerNameProperty': 'FORESTNAME'
}, 'National Forests')

Map.turnOnSelectionAreaCharting()
Map.centerObject(forests)
Map.view()

# Inline: LCMS Land Cover time series over the Uinta-Wasatch-Cache NF
uwc = forests.filter(ee.Filter.stringContains('FORESTNAME', 'Uinta'))

df, fig = cl.summarize_and_chart(
    lcms.select(['Land_Cover']),
    uwc,
    title='LCMS Land Cover - Uinta-Wasatch-Cache NF',
    stacked=False,
)
fig.show()
display(df.head())

# Sankey: LCMS Land Cover transitions for the same forest
sankey_df, sankey_html, matrix_dict = cl.summarize_and_chart(
    lcms.select(['Land_Cover']),
    uwc,
    chart_type='sankey',
    transition_periods=[1985, 2000, 2010, 2024],
    sankey_band_name='Land_Cover',
    title='LCMS Land Cover Transitions - Uinta-Wasatch-Cache NF',
    min_percentage=0.5,
)
display(HTML(cl.sankey_iframe(sankey_html)))
for label, mdf in matrix_dict.items():
    print(f'\n{label}')
    display(mdf)

# Thumbnail: UTM projection with hillshade basemap and geometry fill
crs = gil.getUTMEpsg(ee.Geometry.Point([-111.5, 40.5]))
result = tl.generate_thumbs(
    lcms.select(['Land_Cover']), uwc,
    dimensions=400,
    crs=crs,                        # UTM zone for this location
    basemap='esri-hillshade',       # hillshade basemap
    burn_in_geometry=True,
    geometry_outline_color='yellow',
    geometry_fill_color='00000044', # semi-transparent black outside
    inset_map=True,
    inset_on_map=True,              # overlay inset on the map corner
    title='Uinta-Wasatch-Cache NF',
)
display(HTML(result['html']))

# Same area with NatGeo basemap, no clip, outline arrow style
result2 = tl.generate_thumbs(
    lcms.select(['Land_Cover']), uwc,
    dimensions=400,
    basemap='esri-natgeo',
    overlay_opacity=0.6,
    north_arrow_style='outline',    # outline arrow style
    clip_to_geometry=False,         # bounding box, not clipped
)
display(HTML(result2['html']))

---

8. USFS Ranger Districts (getUSFSDistricts) + MTBS Charting

• Returns USFS Ranger District boundaries

• Optional forest_name and region filters

• We’ll summarize MTBS burn severity across districts

districts = sal.getUSFSDistricts(study_area)
print('Ranger Districts:', districts.aggregate_array('DISTRICTNA').getInfo())

# Filter to a specific forest
uwc_districts = sal.getUSFSDistricts(study_area, forest_name='Uinta-Wasatch-Cache National Forest')
print('\nUinta-Wasatch-Cache districts:', uwc_districts.aggregate_array('DISTRICTNA').getInfo())

# Map: MTBS burn severity with ranger districts
Map.clearMap()

Map.addLayer(mtbs, {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {
        'chartType': 'stacked-bar',
        'visible': [False, True, True, True, True, True, False]
    }
}, 'MTBS Burn Severity')

Map.addSelectLayer(uwc_districts, {
    'strokeColor': 'FF0',
    'selectLayerNameProperty': 'DISTRICTNA'
}, 'Ranger Districts')

Map.turnOnSelectionAreaCharting()
Map.centerObject(uwc_districts)
Map.view()

# Inline: MTBS time series across all ranger districts (per-feature subplots)
per_feature_dfs, fig = cl.summarize_and_chart(
    mtbs.filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    uwc_districts.limit(3),
    feature_label='DISTRICTNA',
    title='MTBS Burn Severity by Ranger District (2000-2024)',
    chart_type = 'stacked_bar',
    class_visible={'Background':False,'Non-Mapping Area':False},
    scale=60,
)
fig.show()

# Show years with fire activity for each district
for name, feat_df in per_feature_dfs.items():
    fire_years = feat_df[feat_df.sum(axis=1) > 0]
    if len(fire_years) > 0:
        print(f'\n--- {name} (years with fire) ---')
        display(fire_years)

# Thumbnail: Multi-feature with dark basemap and classic north arrow
result = tl.generate_thumbs(
    mtbs.filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    uwc_districts.limit(3),
    dimensions=300,
    feature_label='DISTRICTNA',
    basemap='esri-dark-gray',       # dark gray canvas basemap
    north_arrow_style='classic',    # classic arrow style
    burn_in_geometry=True,
    geometry_outline_color='cyan',
    inset_map=True,
    title='MTBS Burn Severity by District',
)
display(HTML(result['html']))

# Same as single-column layout
result2 = tl.generate_thumbs(
    mtbs.filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    uwc_districts.limit(2), columns=1,  # vertical stack
    feature_label='DISTRICTNA',
    basemap='esri-satellite',
    burn_in_geometry=True,
    geometry_outline_color='cyan',
    overlay_opacity=0.85,
)
display(HTML(result2['html']))

---

9. USFS Regions (getUSFSRegions)

• Returns USFS administrative region boundaries

• Properties include REGION, REGIONNAME, REGIONHEAD

regions = sal.getUSFSRegions(study_area)
for feat in regions.getInfo()['features']:
    p = feat['properties']
    print(f"Region {p['REGION']}: {p['REGIONNAME']} (HQ: {p['REGIONHEAD']})")

---

10. US Census Tracts (getUSCensusTracts) + Sentinel-2 Spectral Charting

• Returns TIGER 2020 census tracts

• We’ll use a smaller area and chart mean Sentinel-2 spectral values across tracts

# Use a smaller area for census geographies
small_area = ee.Geometry.Point([-111.9, 40.7]).buffer(5000)

tracts = sal.getUSCensusTracts(small_area)
print(f'Census tracts: {tracts.size().getInfo()}')

block_groups = sal.getUSBlockGroups(small_area)
print(f'Block groups: {block_groups.size().getInfo()}')

# Map: Sentinel-2 composite with census tracts
Map.clearMap()

s2_small = gil.superSimpleGetS2(small_area, '2023-06-01', '2023-09-30').median().clip(small_area)

Map.addLayer(s2_small, {
    'min': 0, 'max': 3000,
    'bands': 'swir1,nir,red',
    'canAreaChart': True,
    'areaChartParams': {'bandNames': 'blue,green,red,nir,swir1,swir2', 'palette': '00D,0D0,D00,D0D,0DD,DD0'}
}, 'Sentinel-2 Summer 2023')

Map.addLayer(lcms.select(['Land_Cover']), {
    'autoViz': True,
    'canAreaChart': True
}, 'LCMS Land Cover', False)

Map.addSelectLayer(tracts, {
    'strokeColor': '0FF',
    'selectLayerNameProperty': 'NAMELSAD'
}, 'Census Tracts')

Map.turnOnSelectionAreaCharting()
Map.centerObject(small_area)
Map.view()

# Inline: Mean spectral values per tract (continuous, reduceRegions)
s2_bands = s2_small.select(['blue', 'green', 'red', 'nir', 'swir1', 'swir2'])

df, fig = cl.summarize_and_chart(
    s2_bands,
    tracts,
    feature_label='NAMELSAD',
    title='Mean S2 Spectral Values by Census Tract',
    palette=['00D', '0D0', 'D00', 'D0D', '0DD', 'DD0'],
    width=800,
    height=400,
    chart_type='stacked_bar'
)
fig.show()
display(df)

# Thumbnail: Sentinel-2 false color with custom viz params on Google Hybrid
result = tl.generate_thumbs(
    s2_small, tracts,
    viz_params={'bands': ['swir1', 'nir', 'red']},
    dimensions=250, max_features=4,
    basemap='google-hybrid',        # Google hybrid basemap
    burn_in_legend=False,           # hide legend for continuous data
    north_arrow=False,              # no north arrow
    scalebar_units='imperial',
)
display(HTML(result['html']))

# True color version with geometry burn-in
result2 = tl.generate_thumbs(
    s2_small, tracts,
    viz_params={'bands': ['red', 'green', 'blue']},
    dimensions=250, max_features=4,
    burn_in_geometry=True,
    geometry_outline_color='red',
    geometry_outline_weight=3,      # thicker boundary
    title='Census Tracts - True Color',
)
display(HTML(result2['html']))

---

11. US Census Blocks (getUSCensusBlocks) and Block Groups (getUSBlockGroups)

• Census blocks are extremely numerous - use small areas only

• Block groups are a coarser aggregation

tiny_area = ee.Geometry.Point([-111.9, 40.7]).buffer(1000)
blocks = sal.getUSCensusBlocks(tiny_area)
print(f'Census blocks in 1km radius: {blocks.size().getInfo()}')

---

12. Protected Areas (getProtectedAreas) + LCMS Change Charting

• Returns WDPA protected area polygons

• Optional iucn_cat and desig_type filters

• We’ll chart LCMS change within protected areas

protected = sal.getProtectedAreas(study_area).sort('REP_AREA', False).limit(5)
print(f'Protected areas: {protected.size().getInfo()}')

# Show a sample of names and designations
for feat in protected.getInfo()['features']:
    p = feat['properties']
    print(f"  {p.get('NAME', 'N/A')} - {p.get('DESIG_ENG', 'N/A')} (IUCN: {p.get('IUCN_CAT', 'N/A')})")

# Map: LCMS Change with protected areas as selectable zones
Map.clearMap()

Map.addLayer(lcms.select(['Change']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True, 'sankey': True}
}, 'LCMS Change')

Map.addLayer(lcms.select(['Land_Cover']), {
    'autoViz': True,
    'canAreaChart': True
}, 'LCMS Land Cover', False)

# Filter to larger protected areas for cleaner visualization
large_protected = protected.filter(ee.Filter.gt('REP_AREA', 10))

Map.addSelectLayer(large_protected, {
    'strokeColor': '0F0',
    'selectLayerNameProperty': 'NAME'
}, 'Protected Areas')

Map.turnOnSelectionAreaCharting()
Map.centerObject(study_area)
Map.view()

# Inline: LCMS Change time series over all protected areas combined
df, fig = cl.summarize_and_chart(
    lcms.select(['Change']),
    protected.limit(2),
    class_visible={'Stable':False,'Non-Processing Area Mask':False},
    title='LCMS Change - Protected Areas near Salt Lake City',
)
fig.show()
for title, d in df.items():
    print(title)
    display(d.head())

# Sankey: LCMS Change transitions over the protected areas
sankey_df, sankey_html, matrix_dict = cl.summarize_and_chart(
    lcms.select(['Change']),
    protected,
    chart_type='sankey',
    transition_periods=[1990, 2000, 2010, 2024],
    sankey_band_name='Change',
    title='LCMS Change Transitions - Protected Areas',
    min_percentage=0.2,
    scale=120 # Increase scale to avoid memory erros - 120m will include every 4th pixel in the stats
)
display(HTML(cl.sankey_iframe(sankey_html)))

for k,t in matrix_dict.items():
    print(k)
    display(t)

lcms_change_composite = lcms.select(['Change']).min().copyProperties(lcms.first())
# Thumbnail: LCMS Change in Albers with USA Topo basemap
result = tl.generate_thumbs(
    lcms_change_composite, protected.limit(2),
    dimensions=300, max_features=2,
    crs='EPSG:5070',               # Albers Equal Area
    basemap='esri-usa-topo',        # USA Topo basemap
    north_arrow_style='outline',
    inset_map=True,
    inset_basemap='esri-natgeo',    # different basemap for inset
    title='LCMS Change Composite - Protected Areas',
)
display(HTML(result['html']))

# Single protected area as ee.Feature
result2 = tl.generate_thumbs(
    lcms_change_composite, protected.first(),  # ee.Feature
    dimensions=400,
    basemap='esri-hillshade-dark',  # dark hillshade
    burn_in_geometry=True,
    geometry_outline_color='lime',
    geometry_fill_color='00000066',
    inset_map=True, inset_on_map=True,
)
display(HTML(result2['html']))

---

13. Roads (getRoads)

• Returns road features from TIGER (US, 2012-2025) or GRIP4 (global)

• TIGER properties include FULLNAME, MTFCC (road type code), RTTYP

• GRIP4 properties include GP_RTP (1=Highway, 2=Primary, 3=Secondary, 4=Tertiary, 5=Local)

• Roads are line features — useful for visualization and proximity analysis

# Use a small area - roads are very numerous
road_area = ee.Geometry.Point([-111.9, 40.7]).buffer(2000)

# TIGER 2024 roads (US, default)
roads = sal.getRoads(road_area)
print(f'TIGER 2024 roads in 2km radius: {roads.size().getInfo()}')

# Filter to primary/secondary roads only (MTFCC codes)
major_roads = roads.filter(ee.Filter.inList('MTFCC', ['S1100', 'S1200']))
print(f'Major roads: {major_roads.size().getInfo()}')
print('Road names:', major_roads.aggregate_array('FULLNAME').distinct().getInfo())

# GRIP4 global roads for the same area
grip_roads = sal.getRoads(road_area, source='grip')
print(f'\nGRIP4 roads: {grip_roads.size().getInfo()}')
highways = grip_roads.filter(ee.Filter.lte('GP_RTP', 2))
print(f'GRIP4 highways+primary: {highways.size().getInfo()}')

# Map: Roads overlaid on Sentinel-2
Map.clearMap()

s2_road = gil.superSimpleGetS2(road_area, '2023-06-01', '2023-09-30').median().clip(road_area)
Map.addLayer(s2_road, {'min': 0, 'max': 3000, 'bands': 'red,green,blue'}, 'Sentinel-2 True Color')

# All roads
Map.addLayer(roads, {'strokeColor': 'AAA'}, 'All Roads', False)

# Major roads highlighted
Map.addLayer(major_roads, {'strokeColor': 'F00'}, 'Major Roads')

Map.centerObject(road_area)
Map.view()

# Thumbnail: Satellite basemap with study area outline and title
result = tl.generate_thumbs(
    s2_road, road_area,
    viz_params=gil.vizParamsFalse10k,
    dimensions=400,
    basemap='esri-satellite',
    burn_in_geometry=True,
    geometry_outline_color='lime',
    inset_map=True,
    title='Roads Study Area - False Color',
)
display(HTML(result['html']))

# True color version on OSM basemap
result2 = tl.generate_thumbs(
    s2_road, road_area,
    viz_params=gil.vizParamsTrue10k,
    dimensions=400,
    basemap='osm',                  # OpenStreetMap
    overlay_opacity=0.5,            # very transparent to see OSM detail
    scalebar=False,
)
display(HTML(result2['html']))

---

14. Buildings (getBuildings)

• Returns building footprints from VIDA Combined (default), Microsoft, or Google Open Buildings

• Automatically determines intersecting countries and merges per-country collections

• Use small areas only - building collections are very large

bldg_area = ee.Geometry.Point([-111.9, 40.7]).buffer(500)

# VIDA Combined buildings
vida_bldgs = sal.getBuildings(bldg_area, source='vida')
print(f'VIDA buildings in 500m radius: {vida_bldgs.size().getInfo()}')

# Microsoft buildings
ms_bldgs = sal.getBuildings(bldg_area, source='ms')
print(f'MS buildings in 500m radius: {ms_bldgs.size().getInfo()}')

# Map: Buildings on Sentinel-2
Map.clearMap()

s2_bldg = gil.superSimpleGetS2(bldg_area, '2023-06-01', '2023-09-30').median().clip(bldg_area)
Map.addLayer(s2_bldg, {'min': 0, 'max': 3000, 'bands': 'red,green,blue'}, 'Sentinel-2 True Color')
Map.addLayer(vida_bldgs, {'strokeColor': '0FF', 'fillColor': '0FF4'}, 'VIDA Buildings')
Map.addLayer(ms_bldgs, {'strokeColor': 'F0F', 'fillColor': 'F0F4'}, 'MS Buildings', False)

Map.centerObject(bldg_area)
Map.view()

# Thumbnail: OSM basemap with transparent overlay (small area)
result = tl.generate_thumbs(
    s2_bldg, bldg_area,
    viz_params=gil.vizParamsFalse10k,
    dimensions=400,
    basemap='osm',
    overlay_opacity=0.6,            # transparent to see OSM roads
    scalebar=False,                 # no scalebar at this zoom
    title='Buildings Study Area',
)
display(HTML(result['html']))

# Light gray basemap version
result2 = tl.generate_thumbs(
    s2_bldg, bldg_area,
    viz_params=gil.vizParamsTrue10k,
    dimensions=400,
    basemap='esri-light-gray',
    overlay_opacity=0.7,
    burn_in_geometry=True,
    geometry_outline_color='red',
)
display(HTML(result2['html']))

---

15. Full Summary: All Layers on One Map

• This final example puts multiple summary area layers together on a single map

• LCMS Land Cover, Land Use, and Change are all available with area charting

• Multiple feature layers are available for selection-based charting

• The Sentinel-2 composite provides a basemap

Map.clearMap()

# Raster layers with area charting
Map.addLayer(s2.filterBounds(counties).median(), {
    'min': 0, 'max': 3000,
    'bands': 'swir1,nir,red',
    'canAreaChart': True,
    'areaChartParams': {'bandNames': 'nir,swir1,swir2', 'palette': 'D0D,0DD,DD0'}
}, 'Sentinel-2 False Color', True)

Map.addLayer(lcms.select(['Land_Cover']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True, 'sankey': True}
}, 'LCMS Land Cover')

Map.addLayer(lcms.select(['Land_Use']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True, 'sankey': True}
}, 'LCMS Land Use', False)

Map.addLayer(lcms.select(['Change']), {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'line': True}
}, 'LCMS Change', False)

Map.addLayer(nlcd, {
    'autoViz': True,
    'canAreaChart': True
}, 'Annual NLCD Land Cover', False)

Map.addLayer(mtbs, {
    'autoViz': True,
    'canAreaChart': True,
    'areaChartParams': {'chartType': 'stacked-bar', 'visible': [False, True, True, True, True, True, False]}
}, 'MTBS Burn Severity', False)

# Selectable summary area layers
Map.addSelectLayer(counties, {'strokeColor': '00F', 'selectLayerNameProperty': 'FULL_NAME'}, 'Counties')
Map.addSelectLayer(forests, {'strokeColor': '0F0', 'selectLayerNameProperty': 'FORESTNAME'}, 'National Forests')
Map.addSelectLayer(uwc_districts, {'strokeColor': 'FF0', 'selectLayerNameProperty': 'DISTRICTNA'}, 'Ranger Districts')
Map.addSelectLayer(urban, {'strokeColor': 'F0F', 'selectLayerNameProperty': 'NAME20'}, 'Urban Areas')

Map.turnOnSelectionAreaCharting()
Map.centerObject(study_area)
Map.view()

---

Inline Summary Across Multiple Summary Areas

• summarize_and_chart with feature_label + ee.ImageCollection produces per-feature time series subplots

• Each feature gets its own subplot, with shared x-axis and consistent legend

• For ee.Image + feature_label, a grouped bar chart is produced instead

• This section demonstrates both chart types across different summary areas

# LCMS Land Cover time series across ranger districts (per-feature subplots)
per_feature_dfs, fig = cl.summarize_and_chart(
    lcms.select(['Land_Cover']).filter(ee.Filter.calendarRange(2000, 2024, 'year')),
    uwc_districts.limit(3),
    feature_label='DISTRICTNA',
    title='LCMS Land Cover by Ranger District (2000-2024)',
    scale=60,
    width=900,
)
fig.show()
for name, feat_df in per_feature_dfs.items():
    print(f'\n--- {name} ---')
    display(feat_df.tail(3))

# LCMS Land Use Sankey over the full study area
sankey_df, sankey_html, matrix_dict = cl.summarize_and_chart(
    lcms.select(['Land_Use']),
    study_area,
    chart_type='sankey',
    transition_periods=[1985, 1990, 2000, 2005, 2015, 2020, 2022, 2024],
    title='LCMS Land Use Transitions - Wasatch Front Region',
    min_percentage=0.5,
)
display(HTML(cl.sankey_iframe(sankey_html)))
for label, mdf in matrix_dict.items():
    print(f'\n{label}')
    display(mdf)

# NLCD grouped bar across counties (ee.Image + feature_label = grouped bar chart)
# Using .mode() to reduce the ImageCollection to a single Image first
nlcd_mode = nlcd.mode().set(nlcd_viz_props)

df, fig = cl.summarize_and_chart(
    nlcd_mode,
    counties.limit(5),
    feature_label='FULL_NAME',
    title='Annual NLCD Land Cover Mode by County - Grouped Bar',
    width=900,
    height=450,
)
fig.show()
display(df)

# Stacked bar chart of LCMS Land Cover (2022-2024) for the first 3 counties

lcms_subset = lcms.select(['Land_Cover']).filter(ee.Filter.calendarRange(2022, 2024, 'year'))
county_subset = counties.limit(3)

stacked_bar_df, stacked_bar_fig = cl.summarize_and_chart(
    nlcd_mode,
    counties.limit(5),
    feature_label='FULL_NAME',
    chart_type='stacked_bar',
    title='Annual NLCD Land Cover Mode by County - Stacked Bar',
    width=900,
    height=450
)
stacked_bar_fig.show()
display(stacked_bar_df)

---

Reports with reportLib

• geeViz.outputLib.reports builds styled HTML reports combining:

  • Thumbnails (via outputLib.thumbs) – auto-visualized satellite imagery

  • Charts (via outputLib.charts) – time series, bar charts, Sankey diagrams

  • Data tables – summary statistics from summarize_and_chart

  • LLM narratives (via Gemini) – AI-generated interpretation of each section

• All EE data requests (charts + thumbnails) run in parallel using thread pools

• Each section has toggle params: generate_table, generate_chart (default True), thumb_format (default "png", also "gif", "filmstrip", or None)

• Key API:

  • Report(title, model, header_text) – create a report

  • report.add_section(ee_obj, geometry, title, ...) – add data sections

  • report.generate(format="html", output_path=...) – generate the report

# Build a report with multiple sections
report = rl.Report(
    title='Wasatch Front Land Assessment',
    header_text='Automated analysis of land cover, land use, and burn severity trends.',
)

# Section 1: LCMS Land Cover in national forests
forests = sal.getUSFSForests(study_area)
report.add_section(
    ee_obj=lcms.select(['Land_Cover']).filter(ee.Filter.calendarRange(2020, 2024, 'year')),
    geometry=forests,
    title='LCMS Land Cover in National Forests',
    scale=60,
    feature_label='FORESTNAME',
)

# Section 2: LCMS Land Use (no thumbnail, just chart + table)
report.add_section(
    ee_obj=lcms.select(['Land_Use']).filter(ee.Filter.calendarRange(2020, 2024, 'year')),
    geometry=study_area,
    title='LCMS Land Use Trends',
    scale=60,
    thumb_format=None,
)

# Section 3: MTBS Burn Severity (no table, just chart + thumb)
report.add_section(
    ee_obj=mtbs.filter(ee.Filter.calendarRange(2010, 2024, 'year')),
    geometry=study_area,
    title='MTBS Burn Severity',
    scale=30,
    generate_table=False,
)

print(f'Report: {report.title}')
print(f'Sections: {len(report._sections)}')
for i, sec in enumerate(report._sections):
    print(f'  {i+1}. {sec.title} (table={sec.generate_table}, chart={sec.generate_chart}, thumb_format={sec.thumb_format})')

# Generate the report (parallel EE requests + Gemini narratives)
import time
start = time.time()

html_path = report.generate(format='html', output_path='wasatch_report.html')

elapsed = time.time() - start
print(f'\nReport generated in {elapsed:.1f}s')
print(f'Saved to: {html_path}')

# Display the report inline in the notebook
from IPython.display import IFrame,HTML
with open('wasatch_report.html', 'r', encoding='utf-8') as f:
    html = f.read()
display(HTML(html))

---

Available Summary Area Types

• getSummaryAreasLib provides a dictionary listing all available summary area functions

import pandas as pd

rows = [{'Type': k, 'Function': v['function'], 'Description': v['description']}
        for k, v in sal.AVAILABLE_SUMMARY_AREAS.items()]
display(pd.DataFrame(rows))