#!/usr/bin/env python3
"""Forecast shoaling in a channel from repeat surveys: how much material is above the depth you
need, where, how fast it's coming back, and when to act. Adds the results to an AO project.

  python3 ao_shoaling.py <project> <folder> [--shoal NAME=FROM:TO] [--as-of DATE] [--backtest]

<project> is the AO project's name, or its id if several projects share the name. <folder> holds
the surveys: a surveys.csv with the columns survey, type, date and grid, the grids it names
(depths below one vertical datum, positive down, all on the same cells) and channel.tif (1 on
the channel's cells), such as the folder ehydro_fetch.py writes.

How it works:
  1. Stations: each cell's distance along the channel, from its --start end, scaled to
     --stations (such as 505+00:950+16, from the survey names) when you know them.
  2. Volume: for every survey that covers a stretch, the material above --required, after
     counting every depth --allowance shallower, as USACE's planning quantities do.
  3. Dredging: a stretch whose volume drops by more than half (and more than --dredge-volume)
     between surveys was dredged. The trend starts again from that survey.
  4. Trend: a straight line through the volumes since the last dredging (at least --min-surveys
     surveys over --min-span days), with a range from the scatter around it (--confidence).
  5. When: the date the line reaches the shoal's trigger: --trigger-volume, or the volume it
     held before it was last dredged. Act by that date minus --lead-time.
This runs for every --segment of channel and for each --shoal you name, or, without --shoal,
for each run of segments that will hold more than --min-volume by the --horizon.

With --as-of, only the surveys up to that date are used, and later surveys are drawn on the
charts. With --backtest, it also makes a forecast from every survey date, checks each against
the surveys 9 to 18 months later, and reports how close they were.

Writes forecast-<as-of>/ in <folder> with segments.csv, shoals.csv, a chart per shoal,
forecast-<date>.tif (each cell's own trend at the horizon, the median over 5 x 5 cells), above-*.tif (the forecast where
the bottom will be above --required), shallower-*.geojson (where the bottom is above --required
now and by when) and backtest.csv. Unless --no-upload, it adds the grids, the latest
--upload-surveys surveys and the shapes to the project as visualizations, and one object per
shoal with its numbers in Notes; running again adds the objects again.

Set AO_URL (such as https://ao.example.com) and AO_EMAIL; the script asks for your password.
On MCP-enabled AO, you can set AO_TOKEN to an API token instead.

This is a sample to adapt, not part of AO and not engineering advice. A straight line can't
foresee storms, floods or other dredging: confirm with a survey before you act.
Needs: pip install requests numpy rasterio scipy matplotlib
"""

import argparse
import csv
import getpass
import hashlib
import json
import os
import re
import sys
from datetime import date, timedelta
from pathlib import Path

import numpy as np
import rasterio
import requests
from rasterio.features import shapes
from rasterio.windows import Window
from rasterio.warp import transform, transform_geom
from scipy import ndimage, sparse, stats
from scipy.sparse.csgraph import dijkstra

NODATA = -9999.0
US_FOOT = 0.3048006096012192  # meters
MIN_PIECE = 10  # cells: smaller pieces of the map layers, and smaller holes in them, are left out
SHALLOWER = [  # where the bottom is above the required depth: (file name, Vector Color)
    ("now", "#d03b3b"),
    ("in-12-months", "#ec835a"),
    ("by-the-horizon", "#fab219"),
]
session = requests.Session()
session.headers["User-Agent"] = "ao-shoaling-script"
refresh_token = None
API = None


def use_tokens(response):
    """Keeps the tokens from a sign-in or renewal's user-context header."""
    global refresh_token
    response.raise_for_status()
    tokens = json.loads(response.headers["user-context"])
    session.headers["Authorization"] = "Bearer " + tokens["auth_token"]
    refresh_token = tokens["refresh_token"]


def sign_in():
    global API
    API = os.environ["AO_URL"].rstrip("/") + "/api/v1"
    if os.environ.get("AO_TOKEN"):
        session.headers["Authorization"] = "Bearer " + os.environ["AO_TOKEN"]
        return
    email = os.environ["AO_EMAIL"]
    password = getpass.getpass(f"AO password for {email}: ")
    use_tokens(
        session.post(
            f"{API}/users/login",
            headers={"X-Email": email, "X-Password": password},
            timeout=60,
        )
    )


def sign_out():
    """Ends a password sign-in's session, so it doesn't count toward an account's 10."""
    if refresh_token:
        try:
            call("POST", "/users/logout", timeout=60)
        except requests.RequestException:
            pass


def call(method, path, upload=None, **options):
    """Calls AO. After a 401, renews the sign-in once and tries again."""
    options.setdefault("timeout", (30, 3600))
    for attempt in (1, 2):
        if upload:
            with open(upload, "rb") as f:
                response = session.request(
                    method, API + path, files={"file": (upload.name, f)}, **options
                )
        else:
            response = session.request(method, API + path, **options)
        if response.status_code == 401 and refresh_token and attempt == 1:
            use_tokens(
                session.post(
                    f"{API}/users/refresh-token",
                    headers={"X-Refresh-Token": refresh_token},
                    timeout=60,
                )
            )
            continue
        response.raise_for_status()
        return response


def find_project(name_or_id):
    """The project's record, which also holds its organization's id (org_id)."""
    me = call("GET", "/users/me").json()["items"]
    projects = call("GET", f"/surveys/users/{me['userid']}/preview").json()["items"][0][
        "surveys"
    ]
    if re.fullmatch(
        r"[0-9a-f]{8}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{12}", name_or_id
    ):
        found = [p for p in projects if p["uuid"] == name_or_id]
    else:
        found = [p for p in projects if p["name"].lower() == name_or_id.lower()]
    if len(found) != 1:
        listed = "".join(
            f"\n  {p['uuid']}  {p['name']} ({p['org_name']})" for p in found
        )
        sys.exit(
            f"{len(found)} projects match {name_or_id!r}. Give the project's id instead.{listed}"
        )
    return found[0]


def sha256(path):
    digest = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def upload(project_id, path):
    """Adds a file to the project. Returns its id, and whether it's new to the project."""
    record = {
        "name": path.name,
        "filePath": "",
        "hash": sha256(path),
        "on_conflict": "SKIP",
    }
    answer = call(
        "POST", f"/surveys/{project_id}/files", json={"files": {"files": [record]}}
    ).json()
    status = answer["items"][0]["fileCreationStatus"]
    if not status:
        sys.exit(
            f"A different file named {path.name} is already in the project: rename or delete it, then run again."
        )
    if status[0]["status"] == "dup":
        print(path.name, "is already in the project")
        return status[0]["fid"], False
    call("POST", f"/file-management/geo/{status[0]['fid']}", upload=path)
    print("uploaded", path.name)
    return status[0]["fid"], True


def visualize(project, file_id, params):
    call(
        "POST",
        f"/surveys/{project['uuid']}/files/{file_id}/render",
        json={"params": {"org_uuid": project["org_id"], **params}},
    )
    print(f"visualizing {params['name']!r}: see the project's Visualizations tab")


def grid_params(name, epsg, vertical_epsg, colors="jet"):
    """A grid of depths in meters on a vertical datum, as AO's visualize window sends it, in one
    of AO's color maps (jet, gold or greyscale). MLLW depth (EPSG 5866) is in meters; AO shows
    the depths in each viewer's display units."""
    params = {
        "name": name,
        "epsg": str(epsg),
        "epsg_authority": "EPSG",
        "data_product_type": "model",
        "data_product_subtype": colors,
        "tileset_type": "3DTILES",
        "elev_units": "METERS",
        "z_positive_down": True,
        "nodata_value": str(NODATA).removesuffix(".0"),
    }
    if vertical_epsg:
        params.update(
            custom_vertical_offset=None,
            vertical_epsg=str(vertical_epsg),
            vertical_authority="EPSG",
            vertical_kind="vertical_crs",
            vertical_value_semantic="depth",
        )
    else:
        params.update(custom_vertical_offset=0, vertical_epsg=None)
    return params


def shapes_params(name, color):
    """Shapes in a GeoJSON file, which is always in WGS 84, drawn in one color."""
    return {
        "name": name,
        "epsg": "4326",
        "epsg_authority": "EPSG",
        "data_product_type": "vector",
        "data_product_subtype": "other",
        "vector_color": color,
        "nodata_value": None,
        "custom_vertical_offset": None,
        "vertical_epsg": None,
    }


def display_copy(path, unit, out, finer=4, strip=256):
    """A copy of a depth grid for AO, in meters, made to look smooth up close: cropped to its
    data, smoothed over about a cell (which also takes out the fine ripples gridding can leave),
    on cells `finer` times smaller, with edges that follow a softened outline of the data rather
    than each cell's steps. It works in strips of rows to keep memory small."""
    with rasterio.open(path) as grid:
        profile, values = grid.profile, grid.read(1)
    rows, cols = np.nonzero(values != NODATA)
    r0, c0 = max(rows.min() - 3, 0), max(cols.min() - 3, 0)
    values = values[r0 : rows.max() + 4, c0 : cols.max() + 4]
    if unit == "ft":
        values = np.where(values == NODATA, NODATA, values * US_FOOT)
    height, width = values.shape
    t = profile["transform"]
    x0, y0 = rasterio.transform.xy(t, r0, c0, offset="ul")
    profile.update(
        height=height * finer,
        width=width * finer,
        transform=rasterio.Affine(t.a / finer, t.b, x0, t.d, t.e / finer, y0),
        dtype="float32",
        nodata=NODATA,
        compress="deflate",
        tiled=True,
        blockxsize=256,
        blockysize=256,
    )
    target = out / f"{path.stem}-m.tif"
    with rasterio.open(target, "w", **profile) as copy:
        for top in range(0, height, strip):
            lo, hi = max(top - 8, 0), min(top + strip + 8, height)
            part, n = values[lo:hi], min(strip, height - top)
            valid = part != NODATA
            block = np.full((n * finer, width * finer), NODATA, "float32")
            if valid.any():
                nearest = ndimage.distance_transform_edt(
                    ~valid, return_distances=False, return_indices=True
                )
                level = ndimage.gaussian_filter(part[tuple(nearest)], 1.0)
                edge = ndimage.gaussian_filter(valid.astype("float32"), 0.8)
                fine, keep = (
                    ndimage.zoom(a, finer, order=1, grid_mode=True, mode="nearest")
                    for a in (level, edge)
                )
                inside = slice((top - lo) * finer, (top - lo + n) * finer)
                block = np.where(keep > 0.5, fine, NODATA)[inside].astype("float32")
            copy.write(
                block, 1, window=Window(0, top * finer, width * finer, n * finer)
            )
    return target


def station_text(feet):
    return f"{int(feet // 100)}+{round(feet % 100):02d}"


def parse_station(text):
    """880+00 or 88000 as feet."""
    hundreds, _, rest = text.partition("+")
    return float(hundreds) * 100 + float(rest or 0) if rest else float(hundreds)


def month(d, horizon=None):
    """A date as a month, or past the horizon when there's none."""
    if d:
        return d.strftime("%b %Y")
    return f"after {horizon:%b %Y}" if horizon else "not within the horizon"


def load(folder):
    """The surveys in date order; each one's depths at the channel's cells; the grids' profile."""
    with open(folder / "surveys.csv", newline="") as f:
        surveys = sorted(csv.DictReader(f), key=lambda s: s["date"])
    for s in surveys:
        s["day"] = date.fromisoformat(s["date"])
    with rasterio.open(folder / "channel.tif") as channel:
        inside, profile = channel.read(1) > 0, channel.profile
    cells = np.flatnonzero(inside)
    depths = np.full((len(surveys), len(cells)), np.nan, dtype="float32")
    for i, s in enumerate(surveys):
        with rasterio.open(folder / s["grid"]) as grid:
            if grid.shape != inside.shape or grid.transform != profile["transform"]:
                sys.exit(f"{s['grid']} isn't on the same cells as channel.tif.")
            depths[i] = grid.read(1, masked=True).filled(np.nan).ravel()[cells]
    return surveys, depths, inside, cells, profile


def stations(inside, cells, cell, start, known):
    """Each channel cell's distance along the channel, in feet from the --start end, measured on
    a coarser grid and scaled to the known stations of both ends, if given."""
    f = max(1, int(round(50 / cell)))
    rows, cols = inside.shape
    pad = np.zeros((-(-rows // f) * f, -(-cols // f) * f), bool)
    pad[:rows, :cols] = inside
    coarse = pad.reshape(pad.shape[0] // f, f, pad.shape[1] // f, f).any((1, 3))
    nodes = np.flatnonzero(coarse)
    index = -np.ones(coarse.size, int)
    index[nodes] = np.arange(len(nodes))
    r, c = np.divmod(nodes, coarse.shape[1])
    edges = []
    for dr, dc in ((0, 1), (1, 0), (1, 1), (1, -1)):
        r2, c2 = r + dr, c + dc
        ok = (r2 < coarse.shape[0]) & (c2 >= 0) & (c2 < coarse.shape[1])
        ok[ok] = coarse[r2[ok], c2[ok]]
        edges.append(
            (
                np.arange(len(nodes))[ok],
                index[r2[ok] * coarse.shape[1] + c2[ok]],
                np.full(ok.sum(), f * cell * np.hypot(dr, dc)),
            )
        )
    i, j, w = (np.concatenate(x) for x in zip(*edges))
    graph = sparse.coo_matrix((w, (i, j)), shape=(len(nodes), len(nodes))).tocsr()
    pick = {
        "east": np.argmax(c),
        "west": np.argmin(c),
        "north": np.argmin(r),
        "south": np.argmax(r),
    }[start]
    from_start = dijkstra(graph, directed=False, indices=pick)
    other = int(np.argmax(np.where(np.isfinite(from_start), from_start, -1)))
    from_other = dijkstra(graph, directed=False, indices=other)
    share = from_start / (from_start + from_other)  # 0 at the start end, 1 at the other
    length = from_start[other]
    first, last = known if known else (0.0, length)
    fine_r, fine_c = np.divmod(cells, cols)
    coarse_index = index[(fine_r // f) * coarse.shape[1] + fine_c // f]
    return first + share[coarse_index] * (last - first), length


def volumes(depths, members, required, cell_volume, min_cover):
    """The material above the required depth in a set of cells, survey by survey (NaN when a
    survey covers less than min_cover of them), and the share of the cells each survey covers."""
    d = depths[:, members]
    cover = (~np.isnan(d)).mean(axis=1)
    total = np.nansum(np.maximum(required - d, 0), axis=1) * cell_volume
    return np.where(cover >= min_cover, total, np.nan), cover


def dredgings(series, cover, types, dredge_volume):
    """Indexes of the surveys that show a stretch was dredged: after-dredging (AD) surveys that
    cover at least 30% of it, and surveys after which its volume dropped by more than half."""
    seen = [i for i in range(len(series)) if np.isfinite(series[i])]
    drops = {
        b
        for a, b in zip(seen, seen[1:])
        if series[a] - series[b] > max(dredge_volume, series[a] / 2)
    }
    after = {i for i in range(len(series)) if types[i] == "AD" and cover[i] >= 0.3}
    return sorted(drops | after)


class Trend:
    """A straight line through (year, volume) points, with a range for a new survey."""

    def __init__(self, x, y, confidence):
        self.n = len(x)
        fit = stats.linregress(x, y)
        self.slope, self.intercept = fit.slope, fit.intercept
        self.mean = x.mean()
        self.sxx = np.sum((x - self.mean) ** 2)
        self.scatter = np.sqrt(
            np.sum((y - self.intercept - self.slope * x) ** 2) / (self.n - 2)
        )
        self.z = stats.t.ppf((1 + confidence) / 2, self.n - 2)

    def at(self, x):
        x = np.asarray(x, float)
        center = self.intercept + self.slope * x
        half = (
            self.z
            * self.scatter
            * np.sqrt(1 + 1 / self.n + (x - self.mean) ** 2 / self.sxx)
        )
        return center, center - half, center + half


def trend_since_dredging(series, cover, types, years, use, options):
    """The trend of a stretch since its last dredging, from the surveys in use."""
    known = np.where(use, series, np.nan)
    events = [
        e for e in dredgings(known, cover, types, options.dredge_volume) if use[e]
    ]
    start = events[-1] if events else 0
    points = [i for i in range(start, len(series)) if use[i] and np.isfinite(series[i])]
    if (
        len(points) < options.min_surveys
        or (years[points[-1]] - years[points[0]]) * 365.25 < options.min_span
    ):
        return None, events
    return Trend(years[points], series[points], options.confidence), events


def reaches(trend, steps, level):
    """The first step at which the line (and the edges of its range) reaches a volume."""
    center, low, high = trend.at(steps)

    def first(curve):
        hit = np.flatnonzero(curve >= level)
        return steps[hit[0]] if len(hit) else None

    return first(center), first(high), first(low)


def cell_trends(depths, years, use, dredge_drop, min_surveys, min_span):
    """Each cell's own straight line through its depths since its last dredging (a deepening of
    more than dredge_drop), as the depth at the as-of date (a) and the change per year (b)."""
    n_cells = depths.shape[1]
    start, last = np.zeros(n_cells, int), np.full(n_cells, np.nan, dtype="float32")
    for i in np.flatnonzero(use):
        d = depths[i]
        seen = ~np.isnan(d)
        drop = seen & ~np.isnan(last) & (d - last > dredge_drop)
        start[drop] = i
        last[seen] = d[seen]
    sums = np.zeros((5, n_cells))
    first, final = np.full(n_cells, np.inf), np.full(n_cells, -np.inf)
    for i in np.flatnonzero(use):
        d = depths[i]
        ok = ~np.isnan(d) & (i >= start)
        t, v = years[i], np.where(ok, d, 0.0)
        sums += [ok, ok * t, v, ok * t * t, v * t]
        first[ok], final[ok] = np.minimum(first[ok], t), np.maximum(final[ok], t)
    n, st, sd, stt, std = sums
    with np.errstate(invalid="ignore", divide="ignore"):
        b = (std - st * sd / n) / (stt - st**2 / n)
        a = sd / n - b * st / n
    fitted = (
        (n >= min_surveys) & ((final - first) * 365.25 >= min_span) & np.isfinite(b)
    )
    return np.where(fitted, a, np.nan), np.where(fitted, b, np.nan), last


def pieces(flags, inside, cells):
    """The cells flagged True as a map of the channel, without pieces or holes of fewer than
    MIN_PIECE cells, so the map layers show areas rather than scattered cells."""
    grid = np.zeros(inside.size, bool)
    grid[cells[flags]] = True
    grid = grid.reshape(inside.shape)
    for keep in (True, False):  # the pieces, then the holes in them
        labels, _ = ndimage.label(grid == keep)
        small = np.bincount(labels.ravel()) < MIN_PIECE
        small[0] = False
        grid[small[labels]] = not keep
    return grid


def smoothed(grid, size=5):
    """Each cell with data as the median of the cells with data within size x size cells of it,
    so a few cells' steep trends don't stick up from the map."""
    values = np.where(grid == NODATA, np.nan, grid)
    rows, cols = np.nonzero(grid != NODATA)
    k = size // 2
    around = [
        values[
            np.clip(rows + i, 0, grid.shape[0] - 1),
            np.clip(cols + j, 0, grid.shape[1] - 1),
        ]
        for i in range(-k, k + 1)
        for j in range(-k, k + 1)
    ]
    out = grid.copy()
    out[rows, cols] = np.nanmedian(around, axis=0)
    return out


def rounded(ring, passes=2):
    """A closed ring with its corners cut (Chaikin's method), so cell edges don't show as steps."""
    points = np.array(ring[:-1], float)
    for _ in range(passes):
        following = np.roll(points, -1, axis=0)
        points = np.hstack(
            [0.75 * points + 0.25 * following, 0.25 * points + 0.75 * following]
        ).reshape(-1, 2)
    return [*points.tolist(), points[0].tolist()]


def outline(flags, inside, cells, profile):
    """The cells flagged True as GeoJSON polygons in WGS 84, as pieces() leaves them."""
    grid = pieces(flags, inside, cells).astype("uint8")
    return [
        {
            "type": "Feature",
            "properties": {},
            "geometry": transform_geom(
                profile["crs"],
                "EPSG:4326",
                {
                    "type": "Polygon",
                    "coordinates": [rounded(r) for r in geometry["coordinates"]],
                },
            ),
        }
        for geometry, value in shapes(
            grid, mask=grid == 1, transform=profile["transform"]
        )
        if value == 1
    ]


def chart(
    path,
    shoal,
    series,
    years,
    surveys,
    use,
    steps_years,
    as_of,
    required,
    unit,
    volume_unit,
    options,
    types,
):
    """The shoal's volume, survey by survey, with the forecast, its range and the trigger."""
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.dates as mdates
    import matplotlib.pyplot as plt

    surface, ink, muted, grid_ink = "#fcfcfb", "#0b0b0b", "#52514e", "#e4e3df"
    blue, orange = "#2a78d6", "#eb6834"
    plt.rcParams.update(
        {
            "font.family": "sans-serif",
            "font.size": 11,
            "axes.edgecolor": grid_ink,
            "xtick.color": muted,
            "ytick.color": muted,
        }
    )
    fig, ax = plt.subplots(figsize=(10, 5.4), dpi=200)
    fig.patch.set_facecolor(surface)
    ax.set_facecolor(surface)
    ax.grid(True, axis="y", color=grid_ink, linewidth=1)
    ax.set_axisbelow(True)
    for side in ("top", "right", "left"):
        ax.spines[side].set_visible(False)
    seen = np.isfinite(series)
    days = [surveys[i]["day"] for i in range(len(surveys))]
    thousands = series / 1000
    trend = shoal["trend"]
    if trend:
        when = [as_of + timedelta(days=float(t) * 365.25) for t in steps_years]
        center, low, high = trend.at(steps_years)
        ax.fill_between(
            when,
            np.maximum(low, 0) / 1000,
            high / 1000,
            color=orange,
            alpha=0.16,
            linewidth=0,
        )
        ax.plot(
            when,
            center / 1000,
            color=orange,
            linewidth=2,
            solid_capstyle="round",
            label=f"Forecast made {as_of:%b %Y}, with its {options.confidence:.0%} range",
        )
    before = [i for i in range(len(series)) if seen[i] and use[i]]
    after = [i for i in range(len(series)) if seen[i] and not use[i]]
    ax.plot(
        [days[i] for i in before],
        thousands[before],
        "o",
        color=blue,
        markersize=7,
        markeredgecolor=surface,
        markeredgewidth=2,
        label="Measured: surveys up to the forecast",
    )
    if after:
        ax.plot(
            [days[i] for i in after],
            thousands[after],
            "o",
            color=surface,
            markersize=7,
            markeredgecolor=blue,
            markeredgewidth=2,
            label="Measured: surveys after it",
        )
    top = ax.get_ylim()[1] * 1.12
    ax.set_ylim(0, top)
    if shoal["trigger"]:
        ax.axhline(shoal["trigger"] / 1000, color=ink, linewidth=1.2)
        ax.annotate(
            shoal["trigger_label"],
            (days[0], shoal["trigger"] / 1000),
            xytext=(0, 7),
            textcoords="offset points",
            color=ink,
            fontsize=10,
            bbox={
                "boxstyle": "square,pad=0.15",
                "facecolor": surface,
                "edgecolor": "none",
            },
            zorder=5,
        )
    if shoal.get("reaches") and shoal["trigger"]:
        end = as_of + timedelta(days=float(steps_years[-1]) * 365.25)
        late = (
            shoal["late"].strftime("%b %Y") if shoal["late"] else f"after {end:%b %Y}"
        )
        ax.plot(
            [shoal["reaches"]],
            [shoal["trigger"] / 1000],
            "o",
            color=orange,
            markersize=9,
            markeredgecolor=surface,
            markeredgewidth=2,
            zorder=6,
        )
        ax.annotate(
            f"Reaches it around {shoal['reaches']:%b %Y}\n(range {shoal['early']:%b %Y} to {late})",
            (shoal["reaches"], shoal["trigger"] / 1000),
            xytext=(-12, 28),
            textcoords="offset points",
            ha="right",
            color=ink,
            fontsize=10,
            bbox={
                "boxstyle": "square,pad=0.15",
                "facecolor": surface,
                "edgecolor": "none",
            },
            zorder=6,
        )
    for i in dredgings(series, shoal["cover"], types, options.dredge_volume):
        ax.axvline(days[i], color=muted, linewidth=1, ymax=0.9)
        ax.annotate(
            f"Dredged\n({days[i]:%b %Y} survey)",
            (days[i], top * 0.9),
            xytext=(0, 4),
            textcoords="offset points",
            ha="center",
            va="bottom",
            color=muted,
            fontsize=9.5,
        )
    for note in options.note:
        when, _, text = note.partition("=")
        day = date.fromisoformat(when if len(when) > 7 else when + "-15")
        ax.annotate(
            text.replace("\\n", "\n"),
            (day, top * 0.78),
            ha="center",
            va="bottom",
            color=muted,
            fontsize=9.5,
        )
    if after:
        ax.axvline(as_of, color=orange, linewidth=1, alpha=0.6, ymax=0.9)
    ax.yaxis.set_major_formatter(
        matplotlib.ticker.FuncFormatter(lambda v, _: f"{v:,.0f}K" if v else "0")
    )
    ax.set_ylabel(f"{volume_unit.capitalize()} above {required:g} {unit}", color=muted)
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
    ax.legend(
        loc="upper left",
        frameon=False,
        fontsize=10,
        labelcolor=ink,
        bbox_to_anchor=(0, -0.08),
        ncol=3,
    )
    title = (
        options.chart_title or f"{shoal['name']}: material above {required:g} {unit}"
    )
    subtitle = (
        f"Stations {station_text(shoal['from'])} to {station_text(shoal['to'])}"
        + (f" \u00b7 {options.chart_subtitle}" if options.chart_subtitle else "")
    )
    ax.text(
        0,
        1.13,
        title,
        transform=ax.transAxes,
        color=ink,
        fontsize=13.5,
        fontweight="bold",
    )
    ax.text(0, 1.065, subtitle, transform=ax.transAxes, color=muted, fontsize=10.5)
    fig.tight_layout()
    fig.subplots_adjust(top=0.84, bottom=0.17)
    fig.savefig(path, facecolor=surface)
    plt.close(fig)
    print("chart:", path)


def backtest(
    depths,
    years,
    surveys,
    types,
    starts,
    segment_cells,
    required,
    cell_volume,
    options,
    out,
):
    """A forecast from every survey date, for every segment, checked against the surveys 9 to 18
    months later (before any dredging): how far off, and how often inside the range."""
    unit = options.units
    volume_unit = "cubic yards" if unit == "ft" else "cubic meters"
    rows = []
    for start, members in zip(starts, segment_cells):
        series, cover = volumes(
            depths, members, required, cell_volume, options.min_cover
        )
        following = dredgings(series, cover, types, options.dredge_volume)
        for origin in range(len(surveys)):
            if not np.isfinite(series[origin]):
                continue
            use = np.arange(len(surveys)) <= origin
            trend, events = trend_since_dredging(
                series, cover, types, years, use, options
            )
            if trend is None:
                continue
            stop = min([e for e in following if e > origin], default=len(surveys))
            for later in range(origin + 1, stop):
                ahead = years[later] - years[origin]
                if np.isfinite(series[later]) and 0.75 <= ahead <= 1.5:
                    center, low, high = trend.at(years[later])
                    rows.append(
                        [
                            station_text(start),
                            surveys[origin]["date"],
                            surveys[later]["date"],
                            round(float(series[later])),
                            round(float(center)),
                            round(float(low)),
                            round(float(high)),
                            low <= series[later] <= high,
                            round(float(series[origin])),
                        ]
                    )
    with open(out / "backtest.csv", "w", newline="") as f:
        table = csv.writer(f)
        table.writerow(
            [
                "segment",
                "forecast_made",
                "survey",
                "measured",
                "forecast",
                "range_low",
                "range_high",
                "in_range",
                "no_change",
            ]
        )
        table.writerows(rows)
    if not rows:
        return None
    validation_chart(
        out / "backtest.png",
        rows,
        options,
        f"{len(rows)} forecasts, each for a {options.segment:,.0f} {unit} segment: {volume_unit} above {required:g} {unit}",
    )
    measured = np.array([r[3] for r in rows], float)
    error = np.abs(np.array([r[4] for r in rows], float) - measured)
    inside = np.array([r[7] for r in rows])
    big = measured > options.backtest_volume
    summary = {
        "forecasts": len(rows),
        "median_error": float(np.median(error)),
        "in_range": float(inside.mean()),
        "big": int(big.sum()),
        "big_median_error": float(np.median(error[big])) if big.any() else None,
        "big_relative_error": float(np.median(error[big] / measured[big]))
        if big.any()
        else None,
        "big_in_range": float(inside[big].mean()) if big.any() else None,
        # the same check for the naive forecast that nothing changes after the last survey
        "big_no_change_error": float(
            np.median(
                np.abs(np.array([r[8] for r in rows], float)[big] - measured[big])
                / measured[big]
            )
        )
        if big.any()
        else None,
    }
    print(
        f"backtest: {summary['forecasts']} forecasts 9 to 18 months ahead; median error {summary['median_error']:,.0f} per segment; {summary['in_range']:.0%} inside the {options.confidence:.0%} range"
    )
    if big.any():
        print(
            f"  segments holding more than {options.backtest_volume:,.0f}: {summary['big']} forecasts, median error {summary['big_relative_error']:.0%} ({summary['big_median_error']:,.0f}), {summary['big_in_range']:.0%} inside the range; assuming no change would be off by {summary['big_no_change_error']:.0%}"
        )
    return summary


def validation_chart(path, rows, options, subtitle):
    """Every backtest forecast against what the later survey measured."""
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    surface, ink, muted, grid_ink, blue = (
        "#fcfcfb",
        "#0b0b0b",
        "#52514e",
        "#e4e3df",
        "#2a78d6",
    )
    forecast = np.array([r[4] for r in rows], float) / 1000
    measured = np.array([r[3] for r in rows], float) / 1000
    big = measured > options.backtest_volume / 1000
    plt.rcParams.update(
        {
            "font.family": "sans-serif",
            "font.size": 11,
            "axes.edgecolor": grid_ink,
            "xtick.color": muted,
            "ytick.color": muted,
        }
    )
    fig, ax = plt.subplots(figsize=(7.2, 6.4), dpi=200)
    fig.patch.set_facecolor(surface)
    ax.set_facecolor(surface)
    ax.grid(True, color=grid_ink, linewidth=1)
    ax.set_axisbelow(True)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)
    top = max(forecast.max(), measured.max()) * 1.05
    ax.plot([0, top], [0, top], color=ink, linewidth=1.2)
    ax.annotate(
        "Forecast = measured",
        (top * 0.97, top * 0.97),
        ha="right",
        va="bottom",
        color=ink,
        fontsize=10,
        rotation=45,
        rotation_mode="anchor",
        transform_rotates_text=True,
    )
    ax.plot(
        forecast[~big],
        measured[~big],
        "o",
        color=muted,
        alpha=0.35,
        markersize=4,
        markeredgewidth=0,
        label=f"Segments holding up to {options.backtest_volume:,.0f}",
    )
    ax.plot(
        forecast[big],
        measured[big],
        "o",
        color=blue,
        markersize=6,
        markeredgecolor=surface,
        markeredgewidth=1.5,
        label=f"Segments holding more than {options.backtest_volume:,.0f}",
    )
    ax.set_xlim(0, top)
    ax.set_ylim(0, top)
    fmt = matplotlib.ticker.FuncFormatter(lambda v, _: f"{v:,.0f}K" if v else "0")
    ax.xaxis.set_major_formatter(fmt)
    ax.yaxis.set_major_formatter(fmt)
    ax.set_xlabel("Forecast, 9 to 18 months ahead", color=muted)
    ax.set_ylabel("Measured by the later survey", color=muted)
    ax.legend(loc="upper left", frameon=False, fontsize=9.5, labelcolor=ink)
    ax.text(
        0,
        1.1,
        options.chart_title_backtest
        or "Backtest: forecasts vs what later surveys measured",
        transform=ax.transAxes,
        color=ink,
        fontsize=13,
        fontweight="bold",
    )
    ax.text(0, 1.045, subtitle, transform=ax.transAxes, color=muted, fontsize=10)
    fig.tight_layout()
    fig.subplots_adjust(top=0.86)
    fig.savefig(path, facecolor=surface)
    plt.close(fig)
    print("chart:", path)


def main():
    parser = argparse.ArgumentParser(
        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
    )
    parser.add_argument("project")
    parser.add_argument("folder", type=Path)
    parser.add_argument(
        "--shoal",
        action="append",
        default=[],
        metavar="NAME=FROM:TO",
        help="a stretch to forecast as one shoal, by station, such as 'Beaumont=880+00:930+00' (repeat for more)",
    )
    parser.add_argument(
        "--required",
        type=float,
        help="the depth to keep (default: channel.json's maintained depth)",
    )
    parser.add_argument(
        "--allowance",
        type=float,
        default=0.5,
        help="count every depth as this much shallower, as USACE's planning quantities do (default 0.5)",
    )
    parser.add_argument(
        "--trigger-volume",
        type=float,
        help="act when a shoal holds this much above the required depth (default: what it held before it was last dredged)",
    )
    parser.add_argument(
        "--lead-time",
        type=int,
        default=180,
        help="days needed to confirm, contract, mobilize and dredge (default 180)",
    )
    parser.add_argument(
        "--horizon", type=int, default=24, help="months to look ahead (default 24)"
    )
    parser.add_argument(
        "--as-of",
        type=date.fromisoformat,
        help="forecast from the surveys up to this date (default: the last survey)",
    )
    parser.add_argument(
        "--stations",
        help="the stations of the two ends, such as 505+00:950+16 (default: distance in feet from the --start end)",
    )
    parser.add_argument(
        "--start",
        choices=("east", "west", "north", "south"),
        default="east",
        help="the end where stations start (default east)",
    )
    parser.add_argument(
        "--segment",
        type=float,
        default=1000,
        help="length of the segments, along the channel (default 1000 ft)",
    )
    parser.add_argument(
        "--min-volume",
        type=float,
        default=10000,
        help="without --shoal, segments that will hold more than this by the horizon form shoals (default 10,000)",
    )
    parser.add_argument(
        "--min-cover",
        type=float,
        default=0.97,
        help="use a survey for a stretch only if it covers this share of it (default 0.97)",
    )
    parser.add_argument(
        "--dredge-volume",
        type=float,
        default=3000,
        help="a drop by more than half that's also more than this counts as dredging (default 3,000)",
    )
    parser.add_argument(
        "--dredge-drop",
        type=float,
        default=1.5,
        help="for the map: a cell that gets deeper by more than this was dredged (default 1.5)",
    )
    parser.add_argument(
        "--min-surveys",
        type=int,
        default=3,
        help="surveys a trend needs since the last dredging (default 3)",
    )
    parser.add_argument(
        "--min-span",
        type=int,
        default=180,
        help="days a trend needs to span (default 180)",
    )
    parser.add_argument(
        "--confidence",
        type=float,
        default=0.8,
        help="the range's confidence (default 0.8)",
    )
    parser.add_argument(
        "--backtest",
        action="store_true",
        help="also check forecasts from every survey date against later surveys",
    )
    parser.add_argument(
        "--backtest-volume",
        type=float,
        default=10000,
        help="also report the backtest for segments holding more than this (default 10,000)",
    )
    parser.add_argument(
        "--vertical-epsg",
        type=int,
        default=5866,
        help="the depths' vertical CRS for AO (default 5866, MLLW depth); 0 for a custom offset of 0",
    )
    parser.add_argument(
        "--units",
        choices=("ft", "m"),
        default="ft",
        help="units of the depths (default ft)",
    )
    parser.add_argument("--chart-title", help="the charts' title")
    parser.add_argument("--chart-title-backtest", help="the backtest chart's title")
    parser.add_argument(
        "--chart-subtitle",
        help="added to the charts' subtitle, such as the data's source",
    )
    parser.add_argument(
        "--note",
        action="append",
        default=[],
        metavar="DATE=TEXT",
        help="a note on the charts at a date, such as '2020-10=Hurricanes Laura and Delta' (repeat for more)",
    )
    parser.add_argument(
        "--upload-surveys",
        type=int,
        default=1,
        help="also add this many of the latest surveys to the project, for context (default 1)",
    )
    parser.add_argument("--no-upload", action="store_true", help="only write the files")
    options = parser.parse_args()

    folder = options.folder
    surveys, depths, inside, cells, profile = load(folder)
    depths -= options.allowance
    channel = (
        json.loads((folder / "channel.json").read_text())
        if (folder / "channel.json").exists()
        else {}
    )
    required = options.required or channel.get("maintained_ft")
    if not required:
        sys.exit("Give the depth to keep with --required.")
    unit = options.units
    volume_unit, per_cubic = (
        ("cubic yards", 27) if unit == "ft" else ("cubic meters", 1)
    )
    cell = abs(profile["transform"].a)
    cell_volume = cell * cell / per_cubic
    as_of = options.as_of or surveys[-1]["day"]
    use = np.array([s["day"] <= as_of for s in surveys])
    types = [s.get("type", "") for s in surveys]
    if use.sum() < options.min_surveys:
        sys.exit(
            f"Only {use.sum()} surveys up to {as_of}: a trend needs {options.min_surveys}."
        )
    years = np.array([(s["day"] - as_of).days / 365.25 for s in surveys])
    steps = np.arange(0, options.horizon * 30.44 + 14, 14) / 365.25
    horizon_date = as_of + timedelta(days=round(steps[-1] * 365.25))
    lead = timedelta(days=options.lead_time)
    epsg = profile["crs"].to_epsg()
    if not epsg:
        sys.exit("The grids' CRS has no EPSG code, which AO needs.")

    known = (
        tuple(parse_station(s) for s in options.stations.split(":"))
        if options.stations
        else None
    )
    station, length = stations(inside, cells, cell, options.start, known)
    named = []
    for text in options.shoal:
        name, _, span = text.partition("=")
        low, high = (parse_station(s) for s in span.split(":"))
        if not low < high or high <= station.min() or low >= station.max():
            sys.exit(
                f"--shoal {text} isn't within the channel's stations, {station_text(station.min())} to {station_text(station.max())}. "
                "Give its span low to high, and the channel's own station numbers with --stations FIRST:LAST."
            )
        named.append((name, low, high))
    print(
        f"{len(surveys)} surveys, {use.sum()} up to {as_of}; channel {length:,.0f} {unit} long; required depth {required:g} {unit}"
    )

    # Segments along the channel
    first = np.floor(station.min() / options.segment) * options.segment
    starts = np.arange(first, station.max(), options.segment)
    if len(starts) > 1 and station.max() - starts[-1] < options.segment / 4:
        starts = starts[:-1]  # a short piece at the end joins the segment before it
    ends = np.append(starts[1:], station.max() + 1)
    segment_cells = [
        np.flatnonzero((station >= s) & (station < e)) for s, e in zip(starts, ends)
    ]
    out = folder / f"forecast-{as_of}"
    out.mkdir(exist_ok=True)
    segment_rows, horizon_volume = [], []
    for start, members in zip(starts, segment_cells):
        series, cover = volumes(
            depths, members, required, cell_volume, options.min_cover
        )
        trend, events = trend_since_dredging(series, cover, types, years, use, options)
        latest = [i for i in range(len(surveys)) if use[i] and np.isfinite(series[i])]
        now = series[latest[-1]] if latest else np.nan
        at_horizon = trend.at(steps[-1])[0] if trend else now
        horizon_volume.append(at_horizon)
        segment_rows.append(
            [
                station_text(start),
                station_text(
                    min(start + options.segment, station.max())
                    if start != starts[-1]
                    else station.max()
                ),
                surveys[latest[-1]]["date"] if latest else "",
                round(float(now)) if np.isfinite(now) else "",
                round(float(trend.slope)) if trend else "",
                round(float(at_horizon)) if np.isfinite(at_horizon) else "",
                surveys[events[-1]]["date"] if events else "",
            ]
        )
    with open(out / "segments.csv", "w", newline="") as f:
        table = csv.writer(f)
        table.writerow(
            [
                "from",
                "to",
                "last_survey",
                f"now_{volume_unit.replace(' ', '_')}",
                "per_year",
                f"at_horizon_{volume_unit.replace(' ', '_')}",
                "last_dredged",
            ]
        )
        table.writerows(segment_rows)

    # Shoals: the ones named, or runs of segments that will hold more than --min-volume
    if not named:
        hot = np.nan_to_num(np.array(horizon_volume)) >= options.min_volume
        k = 0
        while k < len(starts):
            if hot[k]:
                j = k
                while j + 1 < len(starts) and hot[j + 1]:
                    j += 1
                named.append(
                    (
                        f"Sta. {station_text(starts[k])} to {station_text(min(starts[j] + options.segment, station.max()))}",
                        starts[k],
                        min(starts[j] + options.segment, station.max()),
                    )
                )
                k = j + 1
            else:
                k += 1
    shoals = []
    for name, low, high in named:
        members = np.flatnonzero((station >= low) & (station < high))
        series, cover = volumes(
            depths, members, required, cell_volume, options.min_cover
        )
        trend, events = trend_since_dredging(series, cover, types, years, use, options)
        latest = [i for i in range(len(surveys)) if use[i] and np.isfinite(series[i])]
        shoal = {
            "name": name,
            "from": low,
            "to": high,
            "members": members,
            "series": series,
            "cover": cover,
            "trend": trend,
            "events": events,
        }
        before = [
            i
            for i in range(len(series))
            if np.isfinite(series[i]) and use[i] and events and i < events[-1]
        ]
        if options.trigger_volume:
            shoal["trigger"], shoal["trigger_label"] = (
                options.trigger_volume,
                f"Act at {options.trigger_volume:,.0f} {volume_unit}",
            )
        elif before:
            held = series[before[-1]]
            shoal["trigger"] = float(held)
            shoal["trigger_label"] = (
                f"{held:,.0f} {volume_unit}: what it held before the {surveys[events[-1]]['day']:%b %Y} dredging"
            )
        else:
            shoal["trigger"], shoal["trigger_label"] = None, ""
        shoal["now"] = float(series[latest[-1]]) if latest else float("nan")
        shoal["last_survey"] = surveys[latest[-1]]["day"] if latest else None
        shoal["reaches"] = shoal["early"] = shoal["late"] = None
        if trend and shoal["trigger"]:
            shoal["reaches"], shoal["early"], shoal["late"] = (
                None if t is None else as_of + timedelta(days=float(t) * 365.25)
                for t in reaches(trend, steps, shoal["trigger"])
            )
        shoal["act_by"] = shoal["reaches"] - lead if shoal["reaches"] else None
        d = depths[latest[-1], members] if latest else np.array([np.nan])
        shoal["shoalest"] = (
            float(np.nanmin(d)) if np.isfinite(d).any() else float("nan")
        )
        where = (
            members[np.nanargmin(np.where(np.isnan(d), np.inf, d))]
            if np.isfinite(d).any()
            else members[len(members) // 2]
        )
        row, col = divmod(int(cells[where]), inside.shape[1])
        x, y = rasterio.transform.xy(profile["transform"], row, col)
        shoal["lon"], shoal["lat"] = (
            v[0] for v in transform(profile["crs"], "EPSG:4326", [x], [y])
        )
        shoals.append(shoal)

    with open(out / "shoals.csv", "w", newline="") as f:
        table = csv.writer(f)
        table.writerow(
            [
                "shoal",
                "from",
                "to",
                "last_survey",
                f"now_{volume_unit.replace(' ', '_')}",
                "trend_now",
                "per_year",
                "trigger",
                "reaches_trigger",
                "range_from",
                "range_to",
                "act_by",
                "shoalest",
            ]
        )
        for s in shoals:
            t = s["trend"]
            table.writerow(
                [
                    s["name"],
                    station_text(s["from"]),
                    station_text(s["to"]),
                    s["last_survey"],
                    round(s["now"]) if np.isfinite(s["now"]) else "",
                    round(float(t.at(0)[0])) if t else "",
                    round(t.slope) if t else "",
                    round(s["trigger"]) if s["trigger"] else "",
                    s["reaches"] or "",
                    s["early"] or "",
                    s["late"] or "",
                    s["act_by"] or "",
                    f"{s['shoalest']:.1f}",
                ]
            )
            rate = (
                f"trend line at {float(t.at(0)[0]):,.0f} on {as_of}, {t.slope:,.0f} {volume_unit} a year since {surveys[s['events'][-1]]['date'] if s['events'] else 'the first survey'}"
                if t
                else "no trend yet"
            )
            when = (
                f"reaches {s['trigger']:,.0f} {month(s['reaches'])} (range {month(s['early'], horizon_date)} to {month(s['late'], horizon_date)}), act by {month(s['act_by'])}"
                if s["trigger"]
                else "no trigger: set --trigger-volume"
            )
            print(
                f"  {s['name']}: {s['now']:,.0f} {volume_unit} above {required:g} {unit} ({s['last_survey']}), {rate}; {when}"
            )
    for s in shoals:
        chart(
            out
            / f"chart-{re.sub(r'[^a-z0-9]+', '-', s['name'].lower()).strip('-')}.png",
            s,
            s["series"],
            years,
            surveys,
            use,
            steps,
            as_of,
            required,
            unit,
            volume_unit,
            options,
            types,
        )

    # Map: each cell's own trend, for the forecast grid and the shapes
    a, b, last = cell_trends(
        depths, years, use, options.dredge_drop, options.min_surveys, options.min_span
    )
    in_12 = as_of + timedelta(days=365)
    grid_profile = {
        **profile,
        "nodata": NODATA,
        "dtype": "float32",
        "count": 1,
        "compress": "deflate",
    }
    full = np.full(inside.size, NODATA, dtype="float32")
    forecast = np.maximum(a + b * steps[-1] + options.allowance, 0)
    ok = np.isfinite(forecast)
    full[cells[ok]] = forecast[ok]
    full = smoothed(full.reshape(inside.shape))
    grid_path = out / f"forecast-{horizon_date}.tif"
    with rasterio.open(grid_path, "w", **grid_profile) as tif:
        tif.write(full, 1)
    now = last <= required
    soon = ~now & (a + b * 1.0 <= required)
    later = ~now & ~soon & (a + b * steps[-1] <= required)
    # where the bottom will be above the required depth, at its forecast depth: the material itself
    above = pieces(np.nan_to_num(now | soon | later).astype(bool) & ok, inside, cells)
    above_path = out / f"above-{required:g}{unit}-by-{horizon_date}.tif"
    with rasterio.open(above_path, "w", **grid_profile) as tif:
        tif.write(np.where(above, full, NODATA), 1)
    shape_files = []
    for (label, color), flags, title in zip(
        SHALLOWER,
        (now, soon, later),
        (
            f"now ({as_of:%b %Y})",
            f"by {in_12:%b %Y} (forecast)",
            f"by {horizon_date:%b %Y} (forecast)",
        ),
    ):
        features = outline(np.nan_to_num(flags).astype(bool), inside, cells, profile)
        if features:
            path = out / f"shallower-than-{required:g}{unit}-{label}.geojson"
            path.write_text(
                json.dumps({"type": "FeatureCollection", "features": features})
            )
            shape_files.append(
                (path, color, f"Shallower than {required:g} {unit} {title}")
            )

    if options.backtest:
        backtest(
            depths,
            years,
            surveys,
            types,
            starts,
            segment_cells,
            required,
            cell_volume,
            options,
            out,
        )

    if options.no_upload:
        print(f"done: see {out}")
        return
    sign_in()
    try:
        project = find_project(options.project)
        grids = [
            (
                grid_path,
                f"Forecast for {horizon_date:%b %Y} (trend since dredging)",
                "jet",
            ),
            (
                above_path,
                f"Above {required:g} {unit} by {horizon_date:%b %Y} (forecast)",
                "gold",
            ),
        ]
        recent = (
            [i for i in range(len(surveys)) if use[i]][-options.upload_surveys :]
            if options.upload_surveys
            else []
        )
        grids += [
            (
                folder / surveys[i]["grid"],
                f"Survey {surveys[i]['day']:%b %d, %Y} ({surveys[i].get('type') or 'survey'})",
                "jet",
            )
            for i in recent
        ]
        for path, name, colors in grids:
            file_id, new = upload(project["uuid"], display_copy(path, unit, out))
            if new:
                visualize(
                    project,
                    file_id,
                    grid_params(name, epsg, options.vertical_epsg, colors),
                )
        for path, color, name in shape_files:
            file_id, new = upload(project["uuid"], path)
            if new:
                visualize(project, file_id, shapes_params(name, color))
        for s in shoals:
            t = s["trend"]
            notes = (
                f"{s['now']:,.0f} {volume_unit} above {required:g} {unit} on {s['last_survey']}"
                + (
                    f". The trend since it was last dredged is at {float(t.at(0)[0]):,.0f} on {as_of} and rising about {t.slope:,.0f} a year"
                    if t
                    else ""
                )
                + (
                    f". Reaches {s['trigger']:,.0f} around {month(s['reaches'])} (range {month(s['early'], horizon_date)} to {month(s['late'], horizon_date)}); act by {month(s['act_by'])} with a {options.lead_time}-day lead time"
                    if s["trigger"]
                    else ""
                )
                + f". Shallowest {s['shoalest'] + options.allowance:.1f} {unit}. Forecast from the surveys up to {as_of} by ao_shoaling.py."
            )
            call(
                "POST",
                "/marker",
                json={
                    "name": f"{s['name']}: act by {month(s['act_by'])}"
                    if s["act_by"]
                    else s["name"],
                    "type": "Pin",
                    "org_id": project["org_id"],
                    "survey_uuid": project["uuid"],
                    "marker_coordinates": f"POINT({s['lon']:.8f} {s['lat']:.8f})",
                    "label": "Shoaling forecast",
                    "comments": notes,
                    "source": "ao_shoaling.py",
                },
            )
        print(f"added {len(shoals)} objects; done: see {out}")
    finally:
        sign_out()


if __name__ == "__main__":
    main()