Simple normalized difference

nd.py

# nd.py a.tiff b.tiff difference.tiff
# quick and dirty (opinionated and not robust)
# assumes same-size single-band inputs
# output always uint16
# error always 2**15

from sys import argv
import rasterio as rio
from pathlib import Path
import numpy as np

A_path, B_path, dst_path = (Path(x) for x in argv[1:])

assert A_path.exists() and B_path.exists()
assert not dst_path.exists()

with rio.open(A_path) as A_src, rio.open(B_path) as B_src:
    profile = A_src.profile
    assert profile["count"] == 1
    A_nd = A_src.profile["nodata"]
    B_nd = B_src.profile["nodata"]

    A = A_src.read(1).astype(np.float32)
    B = B_src.read(1).astype(np.float32)

np.seterr(divide="ignore")
ratio = (A - B) / (A + B)

mask = np.where((A == A_nd) | (B == B_nd))
ratio[mask] = 0

profile.update({"interleave": "pixel", "nodata": 0})

ratio = (((ratio + 1) / 2) * (2**16 - 1)).clip(0, 2**16 - 1).astype(np.uint16)

with rio.open(dst_path, "w", **profile) as dst:
    dst.write(ratio, 1)