#!/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()