Skip to content

NLCD Land Cover Layer Authoring for RASMapper

Create a RASMapper land-cover sidecar HDF from an NLCD-coded raster, register it in .rasmap, write the authoritative base Manning's n table to the plain-text geometry file, and associate the compiled geometry HDF with the new layer.

This workflow uses BaldEagleCrkMulti2D for reproducibility. The project already includes a local NLCD-coded land-cover raster, so the notebook is offline-safe; production workflows can use Annual NLCD downloads from USGS/MRLC through EarthExplorer, the MRLC viewer or mosaic download site, ScienceBase, Annual NLCD web services, or USGS AWS S3 in us-west-2.

Workflow Scope

Workflow step ras-commander API used here Status
NLCD code table for sidecar authoring ManningsFromLandCover.default_landcover_classification_table() Added for notebook 212; emits RAS-safe class names.
Land-cover sidecar HDF + GeoTIFF RasMap.add_landcover_layer() Creates Raster Map, Variables, sidecar HDF, sidecar GeoTIFF, and .rasmap registration.
RASMapper layer catalog registration RasMap.list_landcover_layers() Confirms the new LandCoverLayer entry in .rasmap.
Authoritative base Manning's n values GeomLandCover.replace_base_mannings_n() Replaces or creates the plain-text LCMann Table= block in .g##.
Enable spatially varied Manning's n GeomStorage.set_2d_flow_area_settings() Sets spatially_varied_mann_on_faces=True so HEC-RAS imports the land cover table.
Geometry HDF layer association RasMap.associate_geometry_layers() Writes /Geometry association attributes for terrain and land cover.
Full HEC-RAS geometry compile RasCmdr.compute_plan(clear_geompre=True) Clears .c## files (preserves .g##.hdf associations), forces geometry recompilation. Opt-in via RUN_HECRAS_PREPROCESSOR.

Key constraint: the sidecar HDF stores pixel-to-class mapping and starter variables for RASMapper display. The plain-text geometry LCMann Table= is the authoritative source for base Manning's n values that HEC-RAS imports during geometry preprocessing. The spatially_varied_mann_on_faces flag must be enabled for HEC-RAS to use the land cover table. Use clear_geompre (not force_geompre) to preserve the .g##.hdf land cover association.

Python
from pathlib import Path
import logging

import h5py
import matplotlib.pyplot as plt
from matplotlib.patches import Patch, Rectangle
import numpy as np
import pandas as pd
import rasterio
from IPython.display import display

from ras_commander import (
    GeomLandCover,
    RasCmdr,
    RasExamples,
    RasMap,
    init_ras_project,
    ras,
)
from ras_commander.geom import GeomStorage
from ras_commander.geom.ManningsFromLandCover import ManningsFromLandCover
from ras_commander.hdf import HdfLandCover

for logger_name in [
    "ras_commander.RasExamples",
    "ras_commander.RasMap",
    "ras_commander.RasPrj",
    "ras_commander.RasUtils",
    "ras_commander.hdf.HdfBase",
    "ras_commander.hdf.HdfLandCover",
]:
    logging.getLogger(logger_name).setLevel(logging.WARNING)

plt.rcParams.update({
    "figure.dpi": 140,
    "savefig.dpi": 180,
    "axes.titlesize": 13,
    "axes.labelsize": 10,
    "legend.fontsize": 8,
})

RAS_VERSION = "7.0"
GEOMETRY_NUMBER = "12"
AUTHORED_LAYER_NAME = "NLCD Land Cover 212"
AUTHORED_HDF_NAME = "LandCover_NLCD_212.hdf"
CELL_SIZE_FT = 100.0
RUN_HECRAS_PREPROCESSOR = True

Extract the Example Project

Use the reproducible Bald Eagle multi-2D example and target geometry g12, the smallest 2D geometry in this project. The smaller mesh keeps sidecar authoring and visual review fast while still exercising the same RASMapper land-cover contracts.

Python
project_path = RasExamples.extract_project('BaldEagleCrkMulti2D')
init_ras_project(project_path, RAS_VERSION)
project_path = Path(project_path)

geom_row = ras.geom_df.loc[
    ras.geom_df["geom_number"].astype(str).str.zfill(2) == GEOMETRY_NUMBER
].iloc[0]
plan_row = ras.plan_df.loc[
    ras.plan_df["geometry_number"].astype(str).str.zfill(2) == GEOMETRY_NUMBER
].iloc[0]

geometry_path = Path(geom_row["full_path"])
geometry_hdf_path = Path(geom_row["hdf_path"])
plan_number = str(plan_row["plan_number"]).zfill(2)
source_nlcd_raster = project_path / "Land Classification" / "LandCover.tif"
authored_hdf_path = project_path / "Land Classification" / AUTHORED_HDF_NAME

summary = pd.DataFrame([
    {"item": "project", "value": project_path.name},
    {"item": "plan", "value": plan_number},
    {"item": "geometry", "value": geometry_path.name},
    {"item": "2D cell count", "value": int(geom_row["mesh_cell_count"])},
    {"item": "source NLCD raster", "value": source_nlcd_raster.name},
])
display(summary)

Clip Extent and Classification Table

The notebook clips the source raster to the target 2D mesh extent with a 500 ft buffer and writes the sidecar at an NLCD-scale 100 ft grid. default_landcover_classification_table() maps NLCD raster codes to sidecar class IDs, RAS-safe class names, starter Manning's n values, and optional percent impervious values.

Python
def mesh_bounds_from_geometry_hdf(hdf_path: Path, buffer_ft: float = 500.0) -> tuple[float, float, float, float]:
    with h5py.File(hdf_path, "r") as hdf_file:
        coords = hdf_file["Geometry/2D Flow Areas/Cell Points"][()]
    left, bottom = coords.min(axis=0)
    right, top = coords.max(axis=0)
    return (
        float(left - buffer_ft),
        float(bottom - buffer_ft),
        float(right + buffer_ft),
        float(top + buffer_ft),
    )

mesh_extent = mesh_bounds_from_geometry_hdf(geometry_hdf_path)
percent_impervious = {
    11: 100.0,
    21: 10.0,
    22: 25.0,
    23: 55.0,
    24: 85.0,
}
classification_table = ManningsFromLandCover.default_landcover_classification_table(
    percent_impervious=percent_impervious,
)
base_mannings_table = pd.concat(
    [
        pd.DataFrame([
            {"Land Cover Name": "NoData", "Base Mannings n Value": 0.035}
        ]),
        classification_table.rename(
            columns={
                "class_name": "Land Cover Name",
                "mannings_n": "Base Mannings n Value",
            }
        )[["Land Cover Name", "Base Mannings n Value"]],
    ],
    ignore_index=True,
)

with rasterio.open(source_nlcd_raster) as src:
    source_profile = {
        "width": src.width,
        "height": src.height,
        "cell_size_ft": abs(float(src.transform.a)),
        "crs": str(src.crs),
    }

display(pd.DataFrame([source_profile]))
display(classification_table.head(12))
assert not classification_table["class_name"].str.contains(r"[\/]", regex=True).any()

Author and Register the Sidecar Layer

RasMap.add_landcover_layer() reprojects/clips/resamples the NLCD-coded raster, writes the RASMapper sidecar HDF and companion GeoTIFF, and registers a distinct LandCoverLayer in the project .rasmap.

Python
authored_hdf = RasMap.add_landcover_layer(
    project_path,
    source_nlcd_raster,
    classification_table,
    cell_size=CELL_SIZE_FT,
    output_hdf_path=authored_hdf_path,
    restrict_to_extent=mesh_extent,
    layer_name=AUTHORED_LAYER_NAME,
)
authored_tif = authored_hdf.with_suffix(".tif")

raster_map = HdfLandCover.get_landcover_raster_map(authored_hdf)
landcover_layers = RasMap.list_landcover_layers(project_path)
registered = landcover_layers.loc[
    landcover_layers["resolved_path"] == str(authored_hdf)
][["name", "filename", "selected_parameter", "classification_kind"]]

with h5py.File(authored_hdf, "r") as hdf_file:
    sidecar_attrs = {
        key: hdf_file.attrs[key].decode("utf-8")
        for key in ["File Type", "LC Type", "Version"]
    }

display(pd.DataFrame([sidecar_attrs]))
display(registered)
display(raster_map.head(12))
assert sidecar_attrs == {"File Type": "HEC Land Cover", "LC Type": "LandCover", "Version": "2.0"}
assert len(registered) == 1

Visual Check of Authored NLCD Classes

The map below is the generated sidecar GeoTIFF. It uses the clipped 2D mesh extent, project coordinates in feet, a north arrow, and a scale bar so the figure can be reviewed as an engineering map rather than a raw array dump.

Python
with rasterio.open(authored_tif) as src:
    landcover_array = src.read(1)
    bounds = src.bounds

unique_codes, counts = np.unique(landcover_array, return_counts=True)
code_to_name = dict(zip(raster_map["pixel_value"], raster_map["class_name"]))
nonzero = [(int(code), int(count)) for code, count in zip(unique_codes, counts) if int(code) != 0]
nonzero_sorted = sorted(nonzero, key=lambda item: item[1], reverse=True)
top_codes = [code for code, _ in nonzero_sorted[:10]]
plot_codes = top_codes + [code for code, _ in nonzero_sorted[10:]]
code_to_plot = {code: index + 1 for index, code in enumerate(top_codes)}
plot_array = np.zeros_like(landcover_array, dtype="int16")
for code, plot_id in code_to_plot.items():
    plot_array[landcover_array == code] = plot_id
if len(nonzero_sorted) > len(top_codes):
    other_mask = np.isin(landcover_array, [code for code, _ in nonzero_sorted[10:]])
    plot_array[other_mask] = len(top_codes) + 1

colors = ["#f2f2f2"] + list(plt.cm.tab20(np.linspace(0, 1, max(len(top_codes), 1))))
if len(nonzero_sorted) > len(top_codes):
    colors.append("#555555")
cmap = plt.matplotlib.colors.ListedColormap(colors)

fig, ax = plt.subplots(figsize=(12.0, 7.0))
extent = (bounds.left, bounds.right, bounds.bottom, bounds.top)
im = ax.imshow(plot_array, extent=extent, origin="upper", cmap=cmap, interpolation="nearest")
ax.set_title("NLCD Land Cover Classes for BaldEagleCr g12")
ax.set_xlabel("State Plane X (ft)")
ax.set_ylabel("State Plane Y (ft)")
ax.ticklabel_format(style="plain", useOffset=False)

legend_handles = [
    Patch(facecolor=colors[i + 1], edgecolor="black", linewidth=0.3,
          label=f"{code}: {code_to_name.get(code, 'Unknown')}")
    for i, code in enumerate(top_codes)
]
if len(nonzero_sorted) > len(top_codes):
    legend_handles.append(Patch(facecolor="#555555", label="Other mapped NLCD classes"))
ax.legend(
    handles=legend_handles,
    title="Dominant classes",
    loc="upper left",
    bbox_to_anchor=(1.01, 1.0),
    borderaxespad=0.0,
    frameon=True,
)

width_ft = bounds.right - bounds.left
height_ft = bounds.top - bounds.bottom
scale_ft = 10000.0
scale_x0 = bounds.left + width_ft * 0.06
scale_y0 = bounds.bottom + height_ft * 0.06
scale_box = Rectangle(
    (scale_x0 - width_ft * 0.018, scale_y0 - height_ft * 0.018),
    scale_ft + width_ft * 0.036,
    height_ft * 0.075,
    facecolor="white",
    edgecolor="black",
    linewidth=0.6,
    alpha=0.88,
    zorder=5,
)
ax.add_patch(scale_box)
ax.plot([scale_x0, scale_x0 + scale_ft], [scale_y0, scale_y0], color="black", linewidth=3, zorder=6)
ax.text(
    scale_x0 + scale_ft / 2,
    scale_y0 + height_ft * 0.025,
    "10,000 ft",
    ha="center",
    va="bottom",
    fontsize=9,
    zorder=7,
)
ax.annotate(
    "N",
    xy=(0.07, 0.90),
    xytext=(0.07, 0.76),
    xycoords="axes fraction",
    ha="center",
    va="center",
    fontsize=11,
    fontweight="bold",
    bbox={"boxstyle": "round,pad=0.2", "facecolor": "white", "edgecolor": "black", "alpha": 0.85},
    arrowprops={"arrowstyle": "-|>", "linewidth": 1.6, "color": "black"},
)

plt.tight_layout(rect=[0.0, 0.0, 0.78, 1.0])
plt.show()

area_summary = pd.DataFrame(
    [
        {
            "NLCD code": code,
            "class": code_to_name.get(code, "Unknown"),
            "cells": count,
            "area_acres": count * CELL_SIZE_FT * CELL_SIZE_FT / 43560.0,
        }
        for code, count in nonzero_sorted
    ]
)
display(area_summary.head(12))

Write the Authoritative Geometry Manning's n Table

The sidecar Variables/ManningsN values are useful for RASMapper display, but the plain-text geometry file is the authoritative source HEC-RAS imports during preprocessing. replace_base_mannings_n() writes the complete class table, including classes that may not appear in this clipped raster but can appear in a full NLCD mosaic.

Python
GeomLandCover.replace_base_mannings_n(
    geometry_path,
    base_mannings_table,
    table_number="16",
)
written_table = GeomLandCover.get_base_mannings_n(geometry_path)

display(written_table.head(12))
assert set(base_mannings_table["Land Cover Name"]) == set(written_table["Land Cover Name"])
assert (geometry_path.with_suffix(geometry_path.suffix + ".bak")).exists()

Enable Spatially Varied Manning's n on 2D Faces

HEC-RAS only imports the LCMann Table into per-cell values when spatially_varied_mann_on_faces is enabled in the 2D flow area settings. Without this flag, HEC-RAS assigns the uniform base Manning's n to every cell regardless of the land cover table.

Python
settings_before = GeomStorage.get_2d_flow_area_settings(geometry_path)
flow_area_name = settings_before.iloc[0]["name"]
print(f"2D flow area: {flow_area_name}")
print(f"  spatially_varied_mann_on_faces BEFORE: {settings_before.iloc[0]['spatially_varied_mann_on_faces']}")

GeomStorage.set_2d_flow_area_settings(
    geometry_path,
    flow_area_name,
    spatially_varied_mann_on_faces=True,
    composite_classification=True,
)

settings_after = GeomStorage.get_2d_flow_area_settings(geometry_path)
print(f"  spatially_varied_mann_on_faces AFTER:  {settings_after.iloc[0]['spatially_varied_mann_on_faces']}")
print(f"  composite_classification AFTER:        {settings_after.iloc[0]['composite_classification']}")
assert settings_after.iloc[0]["spatially_varied_mann_on_faces"] == True

Associate the Geometry HDF with the New Layer

The plain-text .g## file records the Manning's n table. The compiled geometry HDF records the sidecar file association that RASMapper and preprocessing use to find the raster classifications.

Python
terrain_layers = RasMap.list_terrain_layers(project_path)
terrain_hdf = Path(terrain_layers.iloc[0]["resolved_path"])
RasMap.associate_geometry_layers(
    project_path,
    geometry_hdf_path,
    landcover_hdf_path=authored_hdf,
    terrain_hdf_path=terrain_hdf,
)
association = RasMap.get_hdf_geometry_association(geometry_hdf_path)
association_table = pd.DataFrame(
    [
        {"field": "landcover_raw_filename", "value": association.get("landcover_raw_filename")},
        {"field": "landcover_layer_name", "value": association.get("landcover_layer_name")},
        {"field": "terrain_raw_filename", "value": association.get("terrain_raw_filename")},
        {"field": "terrain_layer_name", "value": association.get("terrain_layer_name")},
    ]
)
display(association_table)
assert association["landcover_layer_name"] == AUTHORED_LAYER_NAME
assert Path(association["landcover_hdf_path"]) == authored_hdf

Full HEC-RAS Geometry Compile

Set RUN_HECRAS_PREPROCESSOR = True in the setup cell when running on a workstation with HEC-RAS installed. This step uses RasCmdr.compute_plan() with clear_geompre=True to delete .c## binary preprocessor files while preserving the .g##.hdf (and its land cover association attributes). HEC-RAS then re-reads the modified .g## text during preprocessing.

Why clear_geompre not force_geompre: force_geompre deletes both .g##.hdf and .c##. This destroys the land cover filename attribute that associate_geometry_layers() wrote to the geometry HDF. HEC-RAS needs that attribute to locate the raster sidecar for per-cell Manning's n sampling.

Python
if RUN_HECRAS_PREPROCESSOR:
    # Use clear_geompre=True (NOT force_geompre) so the .g##.hdf is preserved
    # with its land cover association. force_geompre deletes the .g##.hdf,
    # destroying the land cover filename attribute that HEC-RAS needs to find
    # the raster sidecar. clear_geompre deletes only .c## binary files, forcing
    # recompilation while keeping the geometry HDF associations intact.
    RasCmdr.compute_plan(
        plan_number,
        clear_geompre=True,
        num_cores=2,
    )

    # Re-read the geometry HDF to check per-cell Manning's n
    geometry_hdf_path = Path(
        ras.geom_df.loc[
            ras.geom_df["geom_number"].astype(str).str.zfill(2) == GEOMETRY_NUMBER,
            "hdf_path",
        ].iloc[0]
    )
    final_n = HdfLandCover.get_preprocessed_mannings_n(geometry_hdf_path)
    final_summary = final_n.groupby("mannings_n", as_index=False).size()
    display(final_summary)
    assert not final_n.empty
    assert final_n["mannings_n"].nunique() > 1, (
        f"Expected spatially varied Manning's n but got {final_n['mannings_n'].unique()}"
    )
    print(f"Spatially varied Manning's n confirmed: {final_n['mannings_n'].nunique()} distinct values")
else:
    display(pd.DataFrame([
        {
            "preprocessor_step": "skipped",
            "plan_number": plan_number,
            "reason": "Set RUN_HECRAS_PREPROCESSOR=True to run local HEC-RAS text-to-HDF preprocessing.",
        }
    ]))

Notebook Outputs

The authored artifacts are:

  • Land Classification/LandCover_NLCD_212.hdf
  • Land Classification/LandCover_NLCD_212.tif
  • a .rasmap LandCoverLayer named NLCD Land Cover 212
  • an updated geometry text LCMann Table=16
  • geometry HDF /Geometry association attributes pointing at the authored sidecar
CLB Engineering Corporation  ·  LLM Forward Engineering
RAS Commander is a free and open-source project maintained by CLB Engineering Corporation. For agencies and firms seeking to modernize H&H workflows with LLM Forward approaches, contact CLB to partner with the engineers who wrote the automation.