resolved/scripts/fetch-dated-clouds.py

#!/usr/bin/env python3
"""Downloads MODIS clouds for each day of a date range and makes the renderer's cloud
maps for lighting L2 from them.

    scripts/fetch-dated-clouds.py FIRST [LAST]

FIRST and LAST are dates, YYYY-MM-DD; LAST defaults to FIRST. Each day becomes
data/clouds/YYYY-MM-DD.png, 16384 x 8192, one gray channel: the clouds' reflectance
times 255, as data/clouds/blue-marble.png (scripts/fetch-clouds.sh) holds it. A day
already made is skipped.

Source: NASA's Global Imagery Browse Services (GIBS, https://gibs.earthdata.nasa.gov),
which serves the MODIS cloud products (MOD06 on Terra, MYD06 on Aqua) as images, one
global composite per day. Per day it takes six layers, each as the western and
eastern hemispheres at 8192 x 8192 (GIBS times out on the whole map at once):

- cloud optical thickness, standard and partly cloudy retrievals (PCL), from Terra
  and Aqua: images whose colors stand for values in steps of about 4 %, by GIBS's
  color map (data/clouds/source/thickness.xml, fetched with the first day);
- cloud fraction from Terra and Aqua, only to tell where there is no data at all.

A pixel's optical thickness is the first of Terra's standard retrieval, Aqua's,
Terra's PCL and Aqua's PCL that has one; its reflectance is
0.13 tau / (1 + 0.13 tau) (crates/resolved/src/clouds.rs). A pixel with none of them
but a cloud fraction from either satellite is clear: thickness is retrieved for
every cloudy pixel the cloud mask finds in daylight. A pixel with no data from
either, between the satellites' passes or in polar night, takes the static map's
value, so data/clouds/blue-marble.png must be there first.

About 0.3 GB of downloads and a few minutes per day; nothing is kept but the map
(about 20 MB). Needs numpy and Pillow.
"""

import argparse
import concurrent.futures
import datetime
import io
import re
import sys
import time
import urllib.error
import urllib.request
import xml.etree.ElementTree as ET
from pathlib import Path

import numpy as np
from PIL import Image

Image.MAX_IMAGE_PIXELS = None

WMS = "https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi"
COLORMAP = "https://gibs.earthdata.nasa.gov/colormaps/v1.3/MODIS_VIIRS_Cloud_Optical_Thickness.xml"
# Optical thickness in the order taken.
THICKNESS = [
    "MODIS_Terra_Cloud_Optical_Thickness",
    "MODIS_Aqua_Cloud_Optical_Thickness",
    "MODIS_Terra_Cloud_Optical_Thickness_PCL",
    "MODIS_Aqua_Cloud_Optical_Thickness_PCL",
]
FRACTION = ["MODIS_Terra_Cloud_Fraction_Day", "MODIS_Aqua_Cloud_Fraction_Day"]
SIDE = 8192  # pixels on a hemisphere's side
OUT_W, OUT_H = 2 * SIDE, SIDE
# Hemispheres: the column in the map, and the WMS 1.3.0 bounding box in EPSG:4326's
# axis order, latitude first.
HALVES = [(0, "-90,-180,90,0"), (1, "-90,0,90,180")]
REFLECTANCE_SCALE = 0.13


def download(address, tries=4):
    """The body at `address`, retrying with a growing pause when GIBS times out."""
    for attempt in range(tries):
        try:
            with urllib.request.urlopen(address, timeout=300) as response:
                kind = response.headers.get("Content-Type", "")
                body = response.read()
            if address.startswith(WMS) and not kind.startswith("image/png"):
                raise RuntimeError(f"{address}: {kind}: {body[:300]!r}")
            return body
        except (urllib.error.URLError, TimeoutError) as e:
            if attempt == tries - 1:
                raise
            wait = 10 * (attempt + 1)
            print(f"  {e}; retrying in {wait} s", file=sys.stderr)
            time.sleep(wait)


def colormap(path):
    """Sorted packed RGB colors and the optical thickness each stands for: the
    geometric mean of its bin, whose bounds are about 4 % apart."""
    colors, values = [], []
    for entry in ET.parse(path).getroot().iter("ColorMapEntry"):
        if entry.get("nodata") == "true":
            continue
        r, g, b = (int(x) for x in entry.get("rgb").split(","))
        low, high = (float(x) for x in re.findall(r"[-\d.]+", entry.get("value")))
        colors.append((r << 16) | (g << 8) | b)
        values.append((low * high) ** 0.5)
    order = np.argsort(colors)
    colors, values = np.array(colors)[order], np.array(values)[order]
    if len(np.unique(colors)) != len(colors):
        raise RuntimeError(f"{path}: a color stands for more than one value")
    return colors, values


def layer(name, day, box):
    """One hemisphere of a layer on `day` as RGBA."""
    query = (
        f"SERVICE=WMS&REQUEST=GetMap&VERSION=1.3.0&LAYERS={name}&CRS=EPSG:4326"
        f"&BBOX={box}&WIDTH={SIDE}&HEIGHT={SIDE}&FORMAT=image/png&TIME={day}"
    )
    return np.asarray(Image.open(io.BytesIO(download(f"{WMS}?{query}"))).convert("RGBA"))


def decode(rgba, colors, values):
    """Optical thickness per pixel, NaN where the layer has none."""
    packed = (
        (rgba[..., 0].astype(np.int32) << 16)
        | (rgba[..., 1].astype(np.int32) << 8)
        | rgba[..., 2].astype(np.int32)
    )
    opaque = rgba[..., 3] > 0
    index = np.clip(np.searchsorted(colors, packed), 0, len(colors) - 1)
    known = colors[index] == packed
    if (opaque & ~known).any():
        raise RuntimeError(f"{int((opaque & ~known).sum())} pixels of colors not in the map")
    return np.where(opaque, values[index], np.nan).astype(np.float32)


def make(day, root, colors, values, static):
    """The map for `day`, and the fractions of Earth's area that took thickness, were
    clear, and had no data."""
    reflectance = np.empty((OUT_H, OUT_W), np.uint8)
    counts = {"thickness": 0.0, "clear": 0.0, "none": 0.0, "area": 0.0, "sum": 0.0}
    names = THICKNESS + FRACTION
    with concurrent.futures.ThreadPoolExecutor(4) as pool:
        for column, box in HALVES:
            images = dict(zip(names, pool.map(lambda n: layer(n, day, box), names)))
            tau = np.full((SIDE, SIDE), np.nan, np.float32)
            for name in THICKNESS:
                tau = np.where(np.isnan(tau), decode(images.pop(name), colors, values), tau)
            seen = np.zeros((SIDE, SIDE), bool)
            for name in FRACTION:
                seen |= images.pop(name)[..., 3] > 0
            has = ~np.isnan(tau)
            none = ~has & ~seen
            x = REFLECTANCE_SCALE * np.nan_to_num(tau)
            r = np.where(has, x / (1.0 + x), 0.0)
            xs = slice(column * SIDE, (column + 1) * SIDE)
            r = np.where(none, static[:, xs] / 255.0, r)
            reflectance[:, xs] = np.round(255.0 * r).astype(np.uint8)
            lat = np.pi * (0.5 - (np.arange(OUT_H) + 0.5) / OUT_H)
            w = np.cos(lat)[:, None]
            counts["thickness"] += float((w * has).sum())
            counts["clear"] += float((w * (~has & seen)).sum())
            counts["none"] += float((w * none).sum())
            counts["sum"] += float((w * r).sum())
            counts["area"] += float(w.sum()) * SIDE
            print(f"  {day}: {['western', 'eastern'][column]} hemisphere done", flush=True)
    if counts["thickness"] == 0.0:
        raise RuntimeError(f"GIBS has no MODIS cloud thickness for {day}")
    part = root / f"{day}.part.png"
    Image.fromarray(reflectance).save(part)
    part.rename(root / f"{day}.png")
    area = counts["area"]
    return (
        counts["thickness"] / area,
        counts["clear"] / area,
        counts["none"] / area,
        counts["sum"] / area,
    )


def days(first, last):
    day = first
    while day <= last:
        yield day
        day += datetime.timedelta(days=1)


def main():
    parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    parser.add_argument("first", type=datetime.date.fromisoformat)
    parser.add_argument("last", type=datetime.date.fromisoformat, nargs="?")
    args = parser.parse_args()
    last = args.last or args.first
    if last < args.first:
        sys.exit("LAST is before FIRST")
    root = Path(__file__).resolve().parent.parent / "data" / "clouds"
    static_path = root / "blue-marble.png"
    if not static_path.is_file():
        sys.exit(f"{static_path} is missing: run scripts/fetch-clouds.sh first")
    (root / "source").mkdir(parents=True, exist_ok=True)
    thickness = root / "source" / "thickness.xml"
    if not thickness.is_file():
        thickness.write_bytes(download(COLORMAP))
    colors, values = colormap(thickness)
    static = None
    for day in days(args.first, last):
        if (root / f"{day}.png").is_file():
            continue
        if static is None:
            static = np.asarray(Image.open(static_path).convert("L"))
            if static.shape != (OUT_H, OUT_W):
                sys.exit(f"{static_path} is {static.shape[1]} x {static.shape[0]}, not {OUT_W} x {OUT_H}")
        print(f"fetching {day}", flush=True)
        start = time.time()
        retrieved, clear, none, mean = make(day, root, colors, values, static)
        print(
            f"{day}: thickness over {100 * retrieved:.1f} % of Earth's area, clear "
            f"{100 * clear:.1f} %, no data {100 * none:.1f} % (the static map's), mean "
            f"reflectance {mean:.3f}, in {time.time() - start:.0f} s",
            flush=True,
        )


if __name__ == "__main__":
    main()