make.lake

make.lake will modify an initial topobathymetric surface to add new water areas based on user-provided contours. This may be of use in areas that have batymetric surveys in the form of contours. The following is typical usage for make.lake.

A code snippet that uses make.lake may be found in the file digger/examples/pre-run/synthetic/make_lake_synthetic_example.py. In this example, a synthetic lake is placed into flat topography.

"""An example of using digger.make.lake.

This example uses synthetic topography and provides an example of using
digger.make.lake.
"""

import fiona
import matplotlib.pylab as plt
import numpy as np
import rasterio
import shapely
import yaml

from digger import make


# Define the herbie function (Lee and others, 2011) as a means
# to generate contours.
# Lee, H., Gramacy, R., Linkletter, C., and Gray, G., 2011, Optimization
# subject to hidden constraints via statistical emulation. Pacific Journal
# of Optimization, 7(3), 467-478, last accessed June 5, 2025 at
# http://www.ybook.co.jp/online2/oppjo/vol7/p467.html
def herbie(X, Y):
    # standard herbie:
    Z = (
        np.exp(-((X - 1) ** 2))
        + np.exp(-0.8 * (X + 1) ** 2)
        - 0.05 * np.sin(8 * (X + 0.1))
    )
    Z += (
        np.exp(-((Y - 1) ** 2))
        + np.exp(-0.8 * (Y + 1) ** 2)
        - 0.05 * np.sin(8 * (Y + 0.1))
    )
    # additional adjustments to turn herbie into a lake with
    # edges
    Z *= 10
    Z -= 15
    Z[Z < 0] = 0
    return Z


# For illustration purposes, we need to make a raster representing the lake
# we will then write it out as a shapefile and the shapefile to generate
# a lake. In more typical application, the shapefile would be derived from
# surveys.
x = np.linspace(-5, 5, 201)
y = np.linspace(-5, 5, 201)
X, Y = np.meshgrid(x, y)
Z = herbie(X, Y)
clines = np.arange(0, 7, 1)  # contour lines to plot and save

fig, ax = plt.subplots(1, 1, dpi=300)
cs = ax.contour(X, Y, Z, clines, cmap="viridis", linewidths=0.5)
ax.axis("equal")
ax.set_title("Lake Herbie")
cbar = fig.colorbar(cs)
fig.savefig("herbie_contours.png")

# Create a topography file
transform = rasterio.transform.from_bounds(X.min(), Y.min(), X.max(), Y.max(), *Z.shape)
with rasterio.open(
    "lake_herbie_topography.tif",
    "w",
    driver="GTiff",
    height=Z.shape[0],
    width=Z.shape[1],
    count=1,
    dtype=Z.dtype,
    # crs='+proj=latlong',
    transform=transform,
) as sink:
    sink.write(Z * 0, 1)


# Write Lake Herbie to a shapefile
schema = {
    "geometry": "LineString",
    "properties": {"level": "float"},
}
with fiona.open("lake_herbie.shp", "w", driver="ESRI Shapefile", schema=schema) as sink:
    # Loop through contour levels and extract contours using
    # rasterio.features.shapes.
    # This is a slightly different method than using the matplotlib contours
    # but is a more reliable method for generating shapefiles from raster contours.
    for level in clines:
        mask = Z > level
        features = rasterio.features.shapes(
            mask.astype(rasterio.int8), transform=transform
        )

        for polygon_mapping, _ in features:
            polygon = shapely.geometry.shape(polygon_mapping)
            lines = []
            lines.append(polygon.exterior)
            lines.extend(polygon.interiors)
        for line in lines:
            sink.write(
                {
                    "geometry": shapely.geometry.mapping(shapely.LineString(line)),
                    "properties": {"level": float(level)},
                }
            )


# Run digger.make.lake
lake_info = make.lake(
    topo_path="lake_herbie_topography.tif",
    lake_contour_path="lake_herbie.shp",
    depth_field="level",
    depth_unit_convert=1.0,
    write_tif=True,
    write_tt3=False,
    fig_path="lake_herbie.png",
)

# Save lake info to a yaml file
with open("lake_info.yaml", "w") as file:
    yaml.dump(lake_info, file, sort_keys=False)

After this code runs, it will produce geotif and topotype3 files that specify the modified topobathymetry.

It will also create diagnostic figures.

In this example, we generated a lake from the following synthetic contour set:

Fig. 7 The contour set used by make.lake

The standard diagnostic figure from make.lake

Fig. 8 An example of the diagnostic output provided by digger.make.lake.

A dictionary containing some diagnostic information is returned. It contains the following elements:

x1: -2.014925373134328
x2: 1.9154228855721396
y1: -1.9154228855721396
y2: 2.014925373134328
lake_surface_height: 0.0