tif_to_png.py

tif_to_png.py

# DGM TIFF to various PNGs (heightmap, hillshade, and shaded heightmap)
#
# disclaimer: vibecoded!
#

import rasterio
import numpy as np
from PIL import Image

SRC = "dgm1_33340_5596_2_sn.tif"
DST_HEIGHT = "heightmap.png"
DST_SHADE = "hillshade.png"
DST_COMBO = "heightmap_shaded.png"

# Sun position for fake lighting
SUN_AZIMUTH_DEG = 315.0   # NW light (classic cartography)
SUN_ALTITUDE_DEG = 45.0   # sun elevation above horizon
Z_FACTOR = 1.0            # vertical exaggeration

with rasterio.open(SRC) as ds:
    elev = ds.read(1).astype(np.float32)
    nodata = ds.nodata
    # pixel size in meters (DGM1 = 1 m, but read it just in case)
    px, py = ds.res

if nodata is not None:
    elev[elev == nodata] = np.nan

lo = np.nanmin(elev)
hi = np.nanmax(elev)
print(f"Elevation range: {lo:.2f} – {hi:.2f} m")

# --- Normalized heightmap (unchanged) ---
norm = (elev - lo) / (hi - lo)
norm = np.nan_to_num(norm, nan=0.0)
img16 = (norm * 65535).astype(np.uint16)
Image.fromarray(img16, mode="I;16").save(DST_HEIGHT)
print(f"Saved {DST_HEIGHT}")

# --- Hillshade (Horn's method) ---
z = np.nan_to_num(elev, nan=lo) * Z_FACTOR

# Gradients (dz/dx, dz/dy) using central differences
dzdx = np.gradient(z, px, axis=1)
dzdy = np.gradient(z, py, axis=0)

slope = np.arctan(np.hypot(dzdx, dzdy))
aspect = np.arctan2(dzdy, -dzdx)  # standard GIS convention

az = np.deg2rad(360.0 - SUN_AZIMUTH_DEG + 90.0)
alt = np.deg2rad(SUN_ALTITUDE_DEG)

shaded = (np.sin(alt) * np.cos(slope) +
          np.cos(alt) * np.sin(slope) * np.cos(az - aspect))
shaded = np.clip(shaded, 0.0, 1.0)

shade8 = (shaded * 255).astype(np.uint8)
Image.fromarray(shade8, mode="L").save(DST_SHADE)
print(f"Saved {DST_SHADE}")

# --- Combined: colorize heightmap by elevation, multiply by shading ---
# Simple grayscale-tinted version (elevation * lighting)
combo = (norm * shaded * 255).astype(np.uint8)
Image.fromarray(combo, mode="L").save(DST_COMBO)
print(f"Saved {DST_COMBO}")