terrain

Terrain discovery and consolidation for HEC-RAS projects.

Overview

The terrain module discovers named terrain layers from a HEC-RAS project's rasmap configuration and consolidates each surface's TIFF members independently. This is useful for:

• Inspecting terrain configuration: Enumerate all terrain layers, their CRS, resolution, and file locations
• Consolidating terrain: Merge the TIFF members of each named terrain into its own authoritative file
• Downsampling without upsampling: Select a whole native-cell multiple with a 5 ft publication floor
• Recovering relocated projects: Consolidate an explicit, priority-ordered TIFF list when stored RASMapper paths cannot be resolved on the processing host
• Publishing source construction: Export source TIFF footprints and terrain-modification vectors
• Creating HEC-RAS terrain HDFs: Generate new terrain HDF files via RasProcess.exe (required for result mapping)

How Terrain Discovery Works

1. Reads the project's .rasmap file to get terrain names in priority order
2. For each terrain name, locates the corresponding .hdf file in the Terrain/ directory
3. Discovers associated .tif files by matching the HDF stem against TIF file names
4. Optionally reads CRS and resolution from TIF files using rasterio

Terrain Name Matching

TIF files are associated with a terrain by matching the file stem against the terrain name. The matching is case-insensitive and allows suffixes separated by ., _, or -:

TIF Stem	Terrain Name	Match?
Terrain50	Terrain50	Yes (exact)
Terrain50.muncie_clip	Terrain50	Yes (dot separator)
Terrain50_tile2	Terrain50	Yes (underscore separator)
Terrain50-highres	Terrain50	Yes (dash separator)
Terrain50WithChannel	Terrain50	No (alphanumeric continuation)

How Terrain Consolidation Works

1. Discover terrain TIFs from rasmap (priority ordered)
2. Keep terrain names separate: different named RASMapper terrains are never merged implicitly
3. Choose a target grid: preserve native resolution at or above 5 ft; otherwise use the smallest whole native-cell multiple at or above 5 ft. A mixed-resolution mosaic requires an explicit target that is a whole multiple of its coarsest native grid; every source factor is retained in provenance.
4. Merge by windows: reproject each member to the target grid and let the first RASMapper source win in overlaps without allocating the full mosaic in memory
5. Optionally create HEC-RAS terrain HDF via RasTerrain.create_terrain_from_rasters() (requires RasProcess.exe)
6. Optionally register the new terrain in the project's rasmap

Steps 5-6 require RasProcess.exe (Windows or Wine). Steps 1-4 are pure Python (rasterio).

API Reference

ras2cng.terrain.TerrainInfo dataclass

Information about a single terrain layer discovered from a RAS project.

Source code in ras2cng/terrain.py

@dataclass
class TerrainInfo:
    """Information about a single terrain layer discovered from a RAS project."""
    name: str
    hdf_path: Optional[Path] = None
    hdf_exists: bool = False
    tif_files: list[Path] = field(default_factory=list)
    crs: Optional[str] = None
    resolution: Optional[str] = None       # e.g. "50.0 x 50.0 ft"
    bounds: Optional[tuple] = None         # (xmin, ymin, xmax, ymax)
    total_size_mb: float = 0.0

ras2cng.terrain.discover_terrains(project_path)

Discover terrain layers from rasmap in priority order.

Uses RasMap.get_terrain_names() + rasmap_df['terrain_hdf_path']. For each terrain HDF, discovers associated .tif files in same directory.

Parameters:

Name	Type	Description	Default
project_path	Path	Path to .prj file or project directory	required

Returns:

Type	Description
list[TerrainInfo]	List of TerrainInfo in rasmap priority order

Source code in ras2cng/terrain.py

def discover_terrains(project_path: Path) -> list[TerrainInfo]:
    """Discover terrain layers from rasmap in priority order.

    Uses RasMap.get_terrain_names() + rasmap_df['terrain_hdf_path'].
    For each terrain HDF, discovers associated .tif files in same directory.

    Args:
        project_path: Path to .prj file or project directory

    Returns:
        List of TerrainInfo in rasmap priority order
    """
    from ras2cng.project import resolve_project_path
    from ras_commander import init_ras_project

    project_dir, prj_file = resolve_project_path(Path(project_path))
    ras = init_ras_project(project_dir, ras_object="new", load_results_summary=False)

    terrains: list[TerrainInfo] = []

    # Try to get terrain names from rasmap
    terrain_names = _get_terrain_names_safe(project_dir)

    # Try to get terrain HDF paths from rasmap_df
    terrain_hdf_paths: dict[str, Path] = {}
    if ras.rasmap_df is not None and not ras.rasmap_df.empty:
        if "terrain_hdf_path" in ras.rasmap_df.columns:
            for _, row in ras.rasmap_df.iterrows():
                hdf_p = row.get("terrain_hdf_path")
                name = row.get("terrain_name", "")
                if hdf_p and str(hdf_p).strip():
                    terrain_hdf_paths[str(name)] = Path(str(hdf_p))

    # If no rasmap terrain info, fall back to scanning Terrain/ directory
    if not terrain_names and not terrain_hdf_paths:
        terrain_dir = project_dir / "Terrain"
        if terrain_dir.exists():
            hdf_files = sorted(terrain_dir.glob("*.hdf"))
            for hdf_f in hdf_files:
                name = hdf_f.stem
                tif_files = _discover_tifs_for_hdf(hdf_f)
                info = _get_raster_info(tif_files)
                terrains.append(TerrainInfo(
                    name=name,
                    hdf_path=hdf_f,
                    hdf_exists=hdf_f.exists(),
                    tif_files=tif_files,
                    crs=info.get("crs"),
                    resolution=info.get("resolution"),
                    bounds=info.get("bounds"),
                    total_size_mb=sum(f.stat().st_size for f in tif_files if f.exists()) / (1024 * 1024),
                ))
            # Also check for standalone TIFs
            if not hdf_files:
                tif_files = _glob_tifs(terrain_dir)
                if tif_files:
                    info = _get_raster_info(tif_files)
                    terrains.append(TerrainInfo(
                        name="Terrain",
                        tif_files=tif_files,
                        crs=info.get("crs"),
                        resolution=info.get("resolution"),
                        bounds=info.get("bounds"),
                        total_size_mb=sum(f.stat().st_size for f in tif_files if f.exists()) / (1024 * 1024),
                    ))
        return terrains

    # Build terrain info from rasmap data
    seen_names = set()
    for name in terrain_names or list(terrain_hdf_paths.keys()):
        if name in seen_names:
            continue
        seen_names.add(name)

        hdf_path = terrain_hdf_paths.get(name)
        if hdf_path and not hdf_path.is_absolute():
            hdf_path = project_dir / hdf_path

        tif_files = _discover_tifs_for_hdf(hdf_path) if hdf_path else []
        # Also check Terrain/ directory for TIFs matching the name
        if not tif_files:
            terrain_dir = project_dir / "Terrain"
            if terrain_dir.exists():
                all_tifs = _glob_tifs(terrain_dir)
                # Exact stem match: "Terrain" matches "Terrain.tif" and
                # "Terrain_tile2.tif" but NOT "TerrainWithChannel.tif"
                tif_files = sorted(
                    f for f in all_tifs
                    if _stem_matches_name(f.stem, name)
                )

        info = _get_raster_info(tif_files)
        terrains.append(TerrainInfo(
            name=name,
            hdf_path=hdf_path,
            hdf_exists=hdf_path.exists() if hdf_path else False,
            tif_files=tif_files,
            crs=info.get("crs"),
            resolution=info.get("resolution"),
            bounds=info.get("bounds"),
            total_size_mb=sum(f.stat().st_size for f in tif_files if f.exists()) / (1024 * 1024),
        ))

    return terrains

ras2cng.terrain.consolidate_terrain(project_path, output_dir, *, terrain_name='Consolidated', downsample_factor=None, target_resolution=None, terrain_names=None, horizontal_units='Feet', units='Feet', ras_version='6.6', create_hdf=True, register_rasmap=True)

Consolidate project terrains and create a new HEC-RAS terrain HDF.

Full pipeline: 1. Discover terrain TIFFs from rasmap (priority ordered) 2. Merge via rasterio.merge.merge(method='first') -- first wins in overlaps 3. Optionally downsample (reduce resolution) 4. Create HEC-RAS terrain HDF via RasTerrain.create_terrain_from_rasters() 5. Register new terrain in rasmap via RasMap.add_terrain_layer()

Steps 4-5 require RasProcess.exe. If create_hdf=False, only produces the merged TIFF (useful for exporting to cloud-native COG pipeline).

Parameters:

Name	Type	Description	Default
project_path	Path	Path to .prj file or project directory	required
output_dir	Path	Directory for output terrain files	required
terrain_name	str	Name for the consolidated terrain (default: "Consolidated")	'Consolidated'
downsample_factor	Optional[float]	Factor to reduce resolution (2.0 = half resolution)	None
target_resolution	Optional[float]	Target cell size in project units (overrides downsample_factor)	None
terrain_names	Optional[list[str]]	One named terrain to consolidate. Required when the project contains multiple named terrains.	None
horizontal_units	str	Horizontal raster units (Feet or Meters), used to enforce the five-foot publication floor.	'Feet'
units	str	Vertical units "Feet" or "Meters"	'Feet'
ras_version	str	HEC-RAS version for RasProcess.exe	'6.6'
create_hdf	bool	If True, create HEC-RAS terrain HDF (requires RasProcess.exe)	True
register_rasmap	bool	If True, register new terrain in project rasmap	True

Returns:

Type	Description
Path	Path to consolidated terrain HDF (if create_hdf) or TIFF (if not)

Source code in ras2cng/terrain.py

def consolidate_terrain(
    project_path: Path,
    output_dir: Path,
    *,
    terrain_name: str = "Consolidated",
    downsample_factor: Optional[float] = None,
    target_resolution: Optional[float] = None,
    terrain_names: Optional[list[str]] = None,
    horizontal_units: str = "Feet",
    units: str = "Feet",
    ras_version: str = "6.6",
    create_hdf: bool = True,
    register_rasmap: bool = True,
) -> Path:
    """Consolidate project terrains and create a new HEC-RAS terrain HDF.

    Full pipeline:
    1. Discover terrain TIFFs from rasmap (priority ordered)
    2. Merge via rasterio.merge.merge(method='first') -- first wins in overlaps
    3. Optionally downsample (reduce resolution)
    4. Create HEC-RAS terrain HDF via RasTerrain.create_terrain_from_rasters()
    5. Register new terrain in rasmap via RasMap.add_terrain_layer()

    Steps 4-5 require RasProcess.exe. If create_hdf=False, only produces
    the merged TIFF (useful for exporting to cloud-native COG pipeline).

    Args:
        project_path: Path to .prj file or project directory
        output_dir: Directory for output terrain files
        terrain_name: Name for the consolidated terrain (default: "Consolidated")
        downsample_factor: Factor to reduce resolution (2.0 = half resolution)
        target_resolution: Target cell size in project units (overrides downsample_factor)
        terrain_names: One named terrain to consolidate. Required when the
            project contains multiple named terrains.
        horizontal_units: Horizontal raster units (Feet or Meters), used to
            enforce the five-foot publication floor.
        units: Vertical units "Feet" or "Meters"
        ras_version: HEC-RAS version for RasProcess.exe
        create_hdf: If True, create HEC-RAS terrain HDF (requires RasProcess.exe)
        register_rasmap: If True, register new terrain in project rasmap

    Returns:
        Path to consolidated terrain HDF (if create_hdf) or TIFF (if not)
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    # Step 1: Discover terrains
    terrains = discover_terrains(project_path)
    if not terrains:
        raise ValueError("No terrain data found in project")

    # Select exactly one named terrain. TIFF members of that terrain are a
    # mosaic; separate named terrain surfaces are distinct model inputs.
    if terrain_names:
        if len(terrain_names) != 1:
            raise ValueError(
                "Consolidate one named terrain at a time; distinct named terrains "
                "must never be merged"
            )
        name_set = set(terrain_names)
        terrains = [t for t in terrains if t.name in name_set]
        if not terrains:
            raise ValueError(f"No terrains matching names: {terrain_names}")
    elif len(terrains) > 1:
        raise ValueError(
            "Multiple named terrains were discovered. Pass terrain_names=[<name>] "
            "and publish each named surface separately."
        )

    # Collect all TIF files in priority order
    all_tifs: list[Path] = []
    for t in terrains:
        all_tifs.extend(t.tif_files)

    if not all_tifs:
        raise ValueError("No TIFF files found for terrain consolidation")

    # Step 2: Merge TIFFs and record the publication-resolution decision.
    final_tif = consolidate_terrain_files(
        all_tifs,
        output_dir,
        terrain_name=terrain_name,
        source_terrain_name=terrains[0].name,
        downsample_factor=downsample_factor,
        target_resolution=target_resolution,
        horizontal_units=horizontal_units,
    )

    # Step 4: Create HEC-RAS terrain HDF (requires RasProcess.exe)
    if not create_hdf:
        console.print(f"[green]OK[/green] TIFF-only mode: {final_tif}")
        return final_tif

    try:
        from ras_commander import RasTerrain

        from ras2cng.project import resolve_project_path
        project_dir, prj_file = resolve_project_path(Path(project_path))

        terrain_hdf = RasTerrain.create_terrain_from_rasters(
            input_rasters=[final_tif],
            output_folder=output_dir,
            terrain_name=terrain_name,
            units=units,
            hecras_version=ras_version,
        )
        terrain_hdf = Path(terrain_hdf)
        console.print(f"  HEC-RAS terrain HDF -> {terrain_hdf.name}")
    except ImportError:
        console.print("[yellow]Warning:[/yellow] RasTerrain not available; returning TIFF only")
        return final_tif
    except Exception as e:
        console.print(f"[yellow]Warning:[/yellow] Terrain HDF creation failed: {e}")
        console.print("  Returning merged TIFF instead")
        return final_tif

    # Step 5: Register in rasmap
    if register_rasmap:
        try:
            from ras_commander import RasMap

            rasmap_path = project_dir / f"{prj_file.stem}.rasmap"
            if rasmap_path.exists():
                RasMap.add_terrain_layer(
                    terrain_hdf=terrain_hdf,
                    rasmap_path=rasmap_path,
                    layer_name=terrain_name,
                )
                console.print(f"  Registered in rasmap: {rasmap_path.name}")
        except Exception as e:
            console.print(f"[yellow]Warning:[/yellow] Could not register terrain in rasmap: {e}")

    console.print(f"[green]OK[/green] Terrain consolidation complete: {terrain_hdf}")
    return terrain_hdf

ras2cng.terrain.consolidate_terrain_files(tif_files, output_dir, *, terrain_name='Consolidated', source_terrain_name=None, downsample_factor=None, target_resolution=None, horizontal_units='Feet', source_paths=None)

Consolidate an explicit, priority-ordered terrain TIFF mosaic.

This entry point supports projects whose RAS Mapper paths cannot be resolved on the processing host. The first source wins where valid pixels overlap, and the normal no-upsample publication policy is always enforced.

Source code in ras2cng/terrain.py

def consolidate_terrain_files(
    tif_files: list[Path],
    output_dir: Path,
    *,
    terrain_name: str = "Consolidated",
    source_terrain_name: Optional[str] = None,
    downsample_factor: Optional[float] = None,
    target_resolution: Optional[float] = None,
    horizontal_units: str = "Feet",
    source_paths: Optional[list[str | Path]] = None,
) -> Path:
    """Consolidate an explicit, priority-ordered terrain TIFF mosaic.

    This entry point supports projects whose RAS Mapper paths cannot be resolved
    on the processing host. The first source wins where valid pixels overlap,
    and the normal no-upsample publication policy is always enforced.
    """

    sources = [Path(path) for path in tif_files]
    if not sources:
        raise ValueError("No TIFF files provided for terrain consolidation")

    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    source_inventory = inspect_terrain_sources(sources)
    if source_paths is not None:
        if len(source_paths) != len(source_inventory):
            raise ValueError("source_paths must contain one display path per TIFF source")
        for item, display_path in zip(source_inventory, source_paths):
            item["path"] = Path(display_path).as_posix()

    native_resolutions = [
        max(item["resolution_x"], item["resolution_y"])
        for item in source_inventory
    ]
    if downsample_factor is not None and target_resolution is not None:
        raise ValueError("Use either downsample_factor or target_resolution, not both")

    requested_resolution = target_resolution
    if downsample_factor is not None:
        if downsample_factor < 1:
            raise ValueError("downsample_factor must be at least 1; upsampling is prohibited")
        requested_resolution = max(native_resolutions) * float(downsample_factor)

    decision = select_terrain_resolution(
        native_resolutions,
        requested=requested_resolution,
        horizontal_units=horizontal_units,
    )

    console.print(f"[bold]Terrain consolidation:[/bold] {len(sources)} TIFF(s)")
    merged_tif = output_dir / f"{terrain_name}_merged.tif"
    _merge_tifs(
        sources,
        merged_tif,
        target_resolution=decision.target_resolution,
    )
    console.print(f"  Merged -> {merged_tif.name}")

    provenance = {
        "schema": "ras2cng.terrain-consolidation/v1",
        "terrain_name": source_terrain_name or terrain_name,
        "output_name": terrain_name,
        "source_priority": "first-valid-value-wins",
        "resampling": "bilinear",
        "resolution": asdict(decision),
        "sources": source_inventory,
        "output": merged_tif.name,
    }
    provenance_path = output_dir / f"{terrain_name}_terrain-provenance.json"
    provenance_path.write_text(json.dumps(provenance, indent=2) + "\n", encoding="utf-8")
    console.print(f"  Provenance -> {provenance_path.name}")
    return merged_tif

ras2cng.terrain.consolidate_project_terrains(project_path, output_dir, *, target_resolutions=None, horizontal_units='Feet')

Consolidate TIFF members separately for every named project terrain.

Source code in ras2cng/terrain.py

def consolidate_project_terrains(
    project_path: Path,
    output_dir: Path,
    *,
    target_resolutions: Optional[dict[str, float]] = None,
    horizontal_units: str = "Feet",
) -> dict[str, Path]:
    """Consolidate TIFF members separately for every named project terrain."""

    terrains = discover_terrains(project_path)
    if not terrains:
        raise ValueError("No terrain data found in project")
    targets = target_resolutions or {}
    unknown = sorted(set(targets) - {terrain.name for terrain in terrains})
    if unknown:
        raise ValueError(f"Target resolutions reference unknown terrains: {unknown}")

    outputs: dict[str, Path] = {}
    used_names: set[str] = set()
    for terrain in terrains:
        slug = re.sub(r"[^A-Za-z0-9]+", "_", terrain.name).strip("_") or "Terrain"
        output_name = slug
        suffix = 2
        while output_name.lower() in used_names:
            output_name = f"{slug}_{suffix}"
            suffix += 1
        used_names.add(output_name.lower())
        outputs[terrain.name] = consolidate_terrain(
            project_path,
            output_dir,
            terrain_name=output_name,
            terrain_names=[terrain.name],
            target_resolution=targets.get(terrain.name),
            horizontal_units=horizontal_units,
            create_hdf=False,
            register_rasmap=False,
        )
    return outputs

ras2cng.terrain.extract_terrain_source_footprints(tif_files, *, out_crs=None)

Build one queryable footprint polygon per native terrain TIFF member.

Source code in ras2cng/terrain.py

def extract_terrain_source_footprints(
    tif_files: list[Path],
    *,
    out_crs: Optional[str] = None,
):
    """Build one queryable footprint polygon per native terrain TIFF member."""

    import geopandas as gpd
    import pandas as pd
    import rasterio
    from shapely.geometry import box

    frames = []
    for priority, path in enumerate(tif_files):
        path = Path(path)
        with rasterio.open(path) as source:
            if source.crs is None:
                raise ValueError(f"Terrain source has no CRS: {path}")
            frame = gpd.GeoDataFrame(
                {
                    "priority": [priority],
                    "source_file": [path.name],
                    "size_bytes": [path.stat().st_size],
                    "source_crs": [source.crs.to_string()],
                    "resolution_x": [abs(float(source.res[0]))],
                    "resolution_y": [abs(float(source.res[1]))],
                    "width": [int(source.width)],
                    "height": [int(source.height)],
                    "dtype": [str(source.dtypes[0])],
                    "nodata": [None if source.nodata is None else float(source.nodata)],
                },
                geometry=[box(*source.bounds)],
                crs=source.crs,
            )
            if out_crs:
                frame = frame.to_crs(out_crs)
            frames.append(frame)
    if not frames:
        return gpd.GeoDataFrame(geometry=[], crs=out_crs)
    target_crs = frames[0].crs
    normalized = [
        frame if frame.crs == target_crs else frame.to_crs(target_crs)
        for frame in frames
    ]
    return gpd.GeoDataFrame(
        pd.concat(normalized, ignore_index=True),
        geometry="geometry",
        crs=target_crs,
    )

ras2cng.terrain.extract_terrain_modification_layers(terrain_hdf_path, *, crs=None)

Read RASMapper terrain modification vectors from a terrain HDF.

Source code in ras2cng/terrain.py

def extract_terrain_modification_layers(
    terrain_hdf_path: Path,
    *,
    crs: Optional[str] = None,
) -> dict[str, object]:
    """Read RASMapper terrain modification vectors from a terrain HDF."""

    import geopandas as gpd
    import h5py
    import numpy as np
    from shapely.geometry import LineString, Point, Polygon
    from ras_commander.hdf import HdfBase
    from ras_commander.terrain import RasTerrainModWriter

    terrain_hdf_path = Path(terrain_hdf_path)
    if not terrain_hdf_path.is_file():
        raise FileNotFoundError(f"Terrain HDF does not exist: {terrain_hdf_path}")
    source_crs = crs
    if source_crs is None:
        try:
            source_crs = HdfBase.get_projection(terrain_hdf_path)
        except Exception:
            source_crs = None
    if not source_crs:
        raise ValueError(
            f"Terrain modifications require a validated CRS: {terrain_hdf_path}"
        )

    metadata = RasTerrainModWriter.list_modifications(terrain_hdf_path)
    metadata_by_name = {
        str(row["name"]): row.to_dict()
        for _, row in metadata.iterrows()
    }
    lines: list[dict] = []
    polygons: list[dict] = []
    control_points: list[dict] = []
    with h5py.File(terrain_hdf_path, "r") as hdf:
        modifications = hdf.get("Modifications")
        if modifications is None:
            return {
                "terrain_modification_lines": gpd.GeoDataFrame(geometry=[], crs=source_crs),
                "terrain_modification_polygons": gpd.GeoDataFrame(geometry=[], crs=source_crs),
                "terrain_modification_control_points": gpd.GeoDataFrame(geometry=[], crs=source_crs),
            }
        for name, group in modifications.items():
            properties = _terrain_modification_properties(name, group, metadata_by_name.get(name, {}))
            if "Polyline Points" in group:
                points = np.asarray(group["Polyline Points"][:], dtype="float64")
                if points.ndim == 2 and points.shape[1] >= 2 and len(points) >= 2:
                    lines.append({**properties, "geometry": LineString(points[:, :2])})
            elif "Polygon Points" in group and "Polygon Parts" in group:
                points = np.asarray(group["Polygon Points"][:], dtype="float64")
                parts = np.asarray(group["Polygon Parts"][:], dtype="int64")
                rings = [
                    points[start:start + count, :2]
                    for start, count in parts[:, :2]
                    if count >= 3
                ]
                if rings:
                    polygons.append(
                        {
                            **properties,
                            "boundary_elevation_min": _finite_dataset_stat(group, "Boundary Elevations", "min"),
                            "boundary_elevation_max": _finite_dataset_stat(group, "Boundary Elevations", "max"),
                            "geometry": Polygon(rings[0], holes=rings[1:]),
                        }
                    )

            controls = group.get("Control Points")
            if controls is not None and "Points" in controls:
                points = np.asarray(controls["Points"][:], dtype="float64")
                elevations = (
                    np.asarray(controls["Elevations"][:], dtype="float64")
                    if "Elevations" in controls
                    else np.full(len(points), np.nan)
                )
                names = _terrain_control_names(controls, len(points))
                for index, point in enumerate(points):
                    if len(point) < 2:
                        continue
                    control_points.append(
                        {
                            "parent_modification": name,
                            "control_name": names[index],
                            "elevation": float(elevations[index]) if np.isfinite(elevations[index]) else None,
                            "modification_mode": properties.get("modification_mode"),
                            "geometry": Point(float(point[0]), float(point[1])),
                        }
                    )

    return {
        "terrain_modification_lines": gpd.GeoDataFrame(
            lines if lines else {"geometry": []}, geometry="geometry", crs=source_crs
        ),
        "terrain_modification_polygons": gpd.GeoDataFrame(
            polygons if polygons else {"geometry": []}, geometry="geometry", crs=source_crs
        ),
        "terrain_modification_control_points": gpd.GeoDataFrame(
            control_points if control_points else {"geometry": []},
            geometry="geometry",
            crs=source_crs,
        ),
    }