Terrain and bathymetry

Terrain adds ground elevation. Bathymetry adds positive water depth and keeps turbine bases at zero height. Neither moves turbines horizontally, and both are synthetic rather than models of real ground. A SurfaceModel passed to Sampler(..., surface=...) is generated over the boundary bounds.

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

import oowl
from oowl import boundaries as B
from oowl import layouts as L
from oowl import surfaces as SF

Surface models

GaussianTerrain and GaussianBathymetry add a zero-mean Gaussian field to a mean plane with a regional gradient. The field has pointwise scale standard_deviation and squared-exponential covariance \(\exp\!\left(-d^2/2\ell^2\right)\) with \(\ell\) the correlation_length. It is drawn on an internal grid of spacing \(\ell/8\), so resolution does not change the realization, and it is never rescaled or clipped. The grid is drawn by circulant embedding with the FFT, after Dietrich and Newsam, “Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix”, SIAM Journal on Scientific Computing 18 (1997), 1088-1107. It is padded by \(8\ell\), so its covariance is exact to about \(10^{-14}\), and a domain of any size within the surface allocation budget takes a fraction of a second.

PerlinTerrain and PerlinBathymetry blend random unit gradients on a square lattice of cell size correlation_length with the fade \(f(t)=6t^5-15t^4+10t^3\) of Perlin, “Improving noise”, ACM Transactions on Graphics 21 (2002), 681-682. Corner dot products are bounded by \(\sqrt2\) and the blend is convex, so noise \(n\) maps into the range by

\[v=v_{min}+(v_{max}-v_{min})\left(\frac12+\frac{n}{2\sqrt{2}}\right).\]

The range is a bound, not the attained extrema. Equal correlation_length values give different feature sizes in the two models.

boundary = B.OffshorePolygon(area=8_000_000, seed=7)
models = {
    "GaussianTerrain": SF.GaussianTerrain(
        mean=100.0, standard_deviation=20.0, gradient=(0.001, 0.0), seed=7
    ),
    "GaussianBathymetry": SF.GaussianBathymetry(
        mean=50.0, standard_deviation=5.0, seed=8
    ),
    "PerlinTerrain": SF.PerlinTerrain(elevation_range=(20, 180), seed=7),
    "PerlinBathymetry": SF.PerlinBathymetry(
        depth_range=(15, 55), correlation_length=1_500, seed=8
    ),
}
farms = {
    name: oowl.Sampler(boundary, L.Grid(24, seed=3), surface=model).sample()[0]
    for name, model in models.items()
}
fig, axes = plt.subplots(2, 2, figsize=(11, 10))
for ax, (name, farm) in zip(axes.flat, farms.items()):
    farm.surface.plot(ax=ax, title=name)
    farm.plot(ax=ax, facecolor="none", edgecolor="C1", alpha=1.0, color="C1")
fig.tight_layout()

Terrain sets the turbine base z to the ground elevation. Bathymetry keeps it at zero, and surface.at(points) interpolates the depth bilinearly anywhere in the grid.

print("terrain [x, y, z]:", np.round(farms["GaussianTerrain"].coordinates[:2], 1))
offshore = farms["GaussianBathymetry"]
print("bathymetry base z:", offshore.coordinates[:2, 2])
print("water depth [m]:", np.round(offshore.surface.at(offshore.positions[:2]), 1))

terrain [x, y, z]: [[ -977.2 -1090.6    90.3]
 [ -361.8 -1090.6    72.3]]
bathymetry base z: [0. 0.]
water depth [m]: [45.6 41.8]

Sign rules

Terrain is signed elevation, so negative values are basins. Bathymetry requires positive depth at every node, whether generated, built by hand or loaded. PerlinBathymetry rejects a depth_range touching or crossing zero, and GaussianBathymetry raises if a generated depth is not positive.

values = np.full((3, 3), 25.0)
values[1, 1] = -2.0
attempts = {
    "Terrain with a -2 m node": lambda: SF.Terrain([0, 1, 2], [0, 1, 2], values),
    "Bathymetry with a -2 m node": lambda: SF.Bathymetry([0, 1, 2], [0, 1, 2], values),
    "PerlinTerrain((-5, 40))": lambda: SF.PerlinTerrain(elevation_range=(-5, 40)),
    "PerlinBathymetry((0, 40))": lambda: SF.PerlinBathymetry(depth_range=(0, 40)),
}
for label, attempt in attempts.items():
    try:
        attempt()
        print(f"{label:30s} accepted")
    except ValueError as error:
        print(f"{label:30s} ValueError: {error}")

Terrain with a -2 m node       accepted
Bathymetry with a -2 m node    ValueError: bathymetry depths must be positive
PerlinTerrain((-5, 40))        accepted
PerlinBathymetry((0, 40))      ValueError: depth_range must contain only positive depths

Resolution

Grids of different resolution agree exactly at shared nodes for the same boundary and seed. Coarse grids lose detail, so check convergence when slopes or heights matter.

rectangle = B.Rectangle(4, 2, area=8_000_000)
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
surfaces = []
for ax, resolution in zip(axes, [1_000, 250]):
    surface = SF.PerlinTerrain(
        elevation_range=(20, 180), resolution=resolution, seed=7
    ).generate(rectangle)
    surfaces.append(surface)
    surface.plot(ax=ax, title=f"resolution {resolution:g} m")
    ax.collections[0].set_clim(20, 180)
fig.tight_layout()
difference = np.max(np.abs(surfaces[0].values - surfaces[1].values[::4, ::4]))
print(f"largest shared-node difference: {difference:.3g} m")

largest shared-node difference: 0 m

Saving and loading

Bathymetry is saved to windIO’s site.bathymetry as points with their depth, which the windIO site schema accepts, or as a NetCDF grid with surface_format="netcdf", which windIO loads but its schema does not accept. Terrain uses site.elevation and turbine z. Elevation is only a permissive future windIO field, so other consumers may not read it.

with tempfile.TemporaryDirectory() as directory:
    loaded = oowl.load(offshore.save(Path(directory) / "offshore.yaml"))
print(type(loaded.surface).__name__, loaded.surface.values.shape)
print("depths preserved:", np.allclose(loaded.surface.values, offshore.surface.values))

Bathymetry (36, 50)
depths preserved: True