#!/home/robert/weather-env/bin/python """Build GOES-R imagery PNGs at FULL NATIVE RESOLUTION -> OUT_DIR. No downscaling: each product is assembled and written at its native pixel grid (CONUS Ch01 1 km = 5000x3000, true color 0.5 km = 10000x6000, SST 2 km = 2500x1500, FD = the NOAAPORT reduced 1808x904). Rendering is PIL-direct (the data is already in the GOES fixed-grid projection, so we colormap the native array straight to pixels) and political borders are reprojected into the same fixed grid with pyproj and drawn with PIL -- this avoids matplotlib downsampling and keeps memory low enough for this 3.2 GB box. We COMPOSITE the latest file per tile position so gaps backfill across scans. Run every 5 min from cron. """ import os, glob, tempfile, time, traceback import numpy as np import netCDF4 from PIL import Image, ImageDraw, ImageFont import matplotlib from matplotlib import colormaps, colors import cartopy.feature as cfeature from pyproj import Transformer ROOT = "/var/noaaport/data/gempak/images/sat/goesr" OUT_DIR = "/var/www/weather/html/goes" MAX_AGE_MIN = 60 _FONT = os.path.join(matplotlib.get_data_path(), "fonts/ttf/DejaVuSans.ttf") _COAST = cfeature.NaturalEarthFeature("physical", "coastline", "50m") _BORD = cfeature.NaturalEarthFeature("cultural", "admin_0_boundary_lines_land", "50m") _STATES = cfeature.NaturalEarthFeature("cultural", "admin_1_states_provinces_lines", "50m") _LAND = cfeature.NaturalEarthFeature("physical", "land", "50m") # ── tile assembly (native grid) ─────────────────────────────────────────────── def _open(fp): d = open(fp, "rb").read() t = tempfile.NamedTemporaryFile(suffix=".nc", delete=False) t.write(d[d.find(b"\x89HDF"):]); t.close() return netCDF4.Dataset(t.name), t.name def _latest_tiles(product_dir): picks, newest = {}, 0.0 for tdir in glob.glob(os.path.join(product_dir, "*")): fs = glob.glob(os.path.join(tdir, "*")) if not fs: continue f = max(fs, key=os.path.getmtime); mt = os.path.getmtime(f) picks[tdir] = (f, mt); newest = max(newest, mt) cut = newest - MAX_AGE_MIN * 60 return [f for f, mt in picks.values() if mt >= cut], newest def _group_ts(product_dir): """All tile files grouped by scan timestamp: {YYYYMMDD_HHMM: [files]}.""" import re groups = {} for tdir in glob.glob(os.path.join(product_dir, "*")): for f in glob.glob(os.path.join(tdir, "*")): m = re.search(r"_(\d{8}_\d{4})$", f) if m: groups.setdefault(m.group(1), []).append(f) return groups def assemble(rel): tiles, _ = _latest_tiles(os.path.join(ROOT, rel)) return _assemble(tiles) def _assemble(tiles): if not tiles: return None ds0, n0 = _open(tiles[0]) prows = int(ds0.getncattr("product_rows")); pcols = int(ds0.getncattr("product_columns")) ntiles = int(ds0.getncattr("number_product_tiles")) pv = ds0.variables.get("goes_imager_projection") or ds0.variables.get("fixedgrid_projection") H = float(pv.getncattr("perspective_point_height")) lon0 = float(pv.getncattr("longitude_of_projection_origin")) sweep = pv.getncattr("sweep_angle_axis") a = float(pv.getncattr("semi_major_axis")); b = float(pv.getncattr("semi_minor_axis")) ds0.close(); os.unlink(n0) canvas = np.full((prows, pcols), np.nan, np.float32) x0 = y0 = dx = dy = None for f in tiles: try: ds, name = _open(f) cmi = ds.variables["Sectorized_CMI"][:] ro = int(ds.getncattr("tile_row_offset")); co = int(ds.getncattr("tile_column_offset")) th, tw = cmi.shape canvas[ro:ro+th, co:co+tw] = np.ma.filled(cmi.astype(np.float32), np.nan) if x0 is None: xx = ds.variables["x"][:]; yy = ds.variables["y"][:] dx = float(xx[1]-xx[0]); dy = float(yy[1]-yy[0]) x0 = float(xx[0])-co*dx; y0 = float(yy[0])-ro*dy ds.close(); os.unlink(name) except Exception as e: print(f" [{rel}] tile {os.path.basename(f)} failed: {e}") # grid origin/step in projection metres (fixed-grid x/y are microradian scan angle) newest = max(os.path.getmtime(f) for f in tiles) return dict(canvas=canvas, shape=(prows, pcols), ncols=pcols, nrows=prows, H=H, lon0=lon0, sweep=sweep, a=a, b=b, x0_m=x0*1e-6*H, dx_m=dx*1e-6*H, y0_m=y0*1e-6*H, dy_m=dy*1e-6*H, npres=len(tiles), ntiles=ntiles, newest=newest) def _resize(arr, shape): if arr.shape == shape: return arr im = Image.fromarray(np.nan_to_num(arr).astype(np.float32)) return np.asarray(im.resize((shape[1], shape[0]), Image.BILINEAR)) # ── border reprojection + drawing (PIL, native pixels) ──────────────────────── def _transformer(a): crs = (f"+proj=geos +h={a['H']} +lon_0={a['lon0']} +sweep={a['sweep']} " f"+a={a['a']} +b={a['b']} +units=m +no_defs") return Transformer.from_crs("EPSG:4326", crs, always_xy=True) def _px(a, lons, lats, tr): X, Y = tr.transform(np.asarray(lons), np.asarray(lats)) return (X - a["x0_m"]) / a["dx_m"], (Y - a["y0_m"]) / a["dy_m"] # cols, rows def _draw_lines(draw, feat, a, tr, color, w): nc, nr = a["ncols"], a["nrows"] for geom in feat.geometries(): for ls in getattr(geom, "geoms", [geom]): xs, ys = ls.xy cols, rows = _px(a, xs, ys, tr) run = [] for c, r in zip(cols, rows): if np.isfinite(c) and np.isfinite(r) and -nc < c < 2*nc and -nr < r < 2*nr: run.append((float(c), float(r))) else: if len(run) > 1: draw.line(run, fill=color, width=w) run = [] if len(run) > 1: draw.line(run, fill=color, width=w) def _fill_land(draw, a, tr, color): for geom in _LAND.geometries(): for poly in getattr(geom, "geoms", [geom]): cols, rows = _px(a, *poly.exterior.xy, tr) pts = [(float(c), float(r)) for c, r in zip(cols, rows) if np.isfinite(c) and np.isfinite(r)] if len(pts) >= 3: draw.polygon(pts, fill=color) def _overlay(img, a, color, land_black=False): tr = _transformer(a) d = ImageDraw.Draw(img) w = max(1, a["ncols"] // 3000 + 1) if land_black: _fill_land(d, a, tr, (0, 0, 0)) _draw_lines(d, _COAST, a, tr, color, w) _draw_lines(d, _BORD, a, tr, color, w) _draw_lines(d, _STATES, a, tr, color, max(1, w-1)) return d def _label(draw, a, text, color): fs = max(14, a["ncols"] // 110) try: font = ImageFont.truetype(_FONT, fs) except Exception: font = ImageFont.load_default() draw.text((fs//2, fs//2), text, fill=color, font=font) def _save(img, outfile): os.makedirs(OUT_DIR, exist_ok=True) out = os.path.join(OUT_DIR, outfile); tmp = out + ".tmp.png" img.save(tmp); os.replace(tmp, out) print(f" {outfile}: {img.size[0]}x{img.size[1]} px ({os.path.getsize(out)} b)") # ── products ────────────────────────────────────────────────────────────────── def build_single(rel, outfile, title): a = assemble(rel) if a is None: print(f" [{rel}] no tiles"); return v = np.sqrt(np.clip(np.nan_to_num(a["canvas"]), 0, 1.0)) # reflectance, gamma img = Image.fromarray((v * 255).astype(np.uint8), "L").convert("RGB") d = _overlay(img, a, (0, 0, 0)) # black borders _label(d, a, f"{title} {time.strftime('%Y-%m-%d %H:%MZ', time.gmtime(a['newest']))} ({a['npres']}/{a['ntiles']})", (255, 255, 0)) _save(img, outfile) def build_ir(rel, outfile, title, vmin=180.0, vmax=310.0): """Inverted-grayscale IR from a brightness-temperature band (cold cloud tops bright/white, warm surface dark). Used for the high-res FD band 13.""" a = assemble(rel) if a is None: print(f" [{rel}] no tiles"); return bt = np.nan_to_num(a["canvas"], nan=vmax) # off-disk -> warm -> black g = 1.0 - np.clip((bt - vmin) / (vmax - vmin), 0, 1) # cold = white img = Image.fromarray((g * 255).astype(np.uint8), "L").convert("RGB") d = _overlay(img, a, (255, 255, 0)) # yellow borders (read on grayscale) _label(d, a, f"{title} {time.strftime('%Y-%m-%d %H:%MZ', time.gmtime(a['newest']))} ({a['npres']}/{a['ntiles']})", (255, 255, 0)) _save(img, outfile) def build_truecolor(outfile="goes_east_conus_truecolor.png"): aB = assemble("TIR/TIRE01"); aR = assemble("TIR/TIRE02"); aV = assemble("TIR/TIRE03") if not (aB and aR and aV): print(" [truecolor] missing channel(s)"); return shp = aR["shape"] # finest = C02 0.5 km gm = lambda c: (np.clip(np.nan_to_num(_resize(c, shp)), 0, 1) ** (1/2.2)).astype(np.float32) R = gm(aR["canvas"]); B = gm(aB["canvas"]); V = gm(aV["canvas"]) G = 0.45*R + 0.1*V + 0.45*B; del V c = 70.0; F = (259*(c+255))/(255*(259-c)) to8 = lambda x: (np.clip(F*(np.clip(x, 0, 1)-0.5)+0.5, 0, 1) * 255).astype(np.uint8) rgb = np.dstack([to8(R), to8(G), to8(B)]); del R, G, B img = Image.fromarray(rgb, "RGB"); del rgb a = aR; a["npres"] = min(aB["npres"], aR["npres"], aV["npres"]) d = _overlay(img, a, (255, 255, 255)) # white borders _label(d, a, f"GOES-East True Color (CONUS, 0.5km) {time.strftime('%Y-%m-%d %H:%MZ', time.gmtime(a['newest']))}", (255, 255, 255)) _save(img, outfile) def _colorbar(img, cmap, vmin, vmax, label): nc = img.size[0] bw, bh = max(18, nc//90), max(120, img.size[1]//4) x = img.size[0] - bw - max(60, nc//40); y = img.size[1] - bh - max(40, nc//50) grad = np.linspace(vmax, vmin, bh)[:, None].repeat(bw, 1) bar = (cmap(colors.Normalize(vmin, vmax)(grad))[..., :3] * 255).astype(np.uint8) img.paste(Image.fromarray(bar, "RGB"), (x, y)) d = ImageDraw.Draw(img) fs = max(12, nc//150) try: font = ImageFont.truetype(_FONT, fs) except Exception: font = ImageFont.load_default() d.rectangle([x, y, x+bw, y+bh], outline=(255, 255, 255)) for frac, val in [(0, vmax), (0.5, (vmin+vmax)/2), (1, vmin)]: d.text((x+bw+4, y+int(frac*bh)-fs//2), f"{val:.0f}", fill=(255, 255, 255), font=font) d.text((x, y-fs-4), label, fill=(255, 255, 255), font=font) if __name__ == "__main__": jobs = [ lambda: build_single("TIR/TIRE01", "goes_east_conus_ch1.png", "GOES-East Ch01 CONUS 1km"), lambda: build_single("TIR/TIRS01", "goes_east_fd_ch1.png", "GOES-East Ch01 Full Disk"), lambda: build_ir("TIR/TIRS13", "goes_east_fd_ir.png", "GOES-East Ch13 (10.3um IR) Full Disk 2km"), build_truecolor, ] for j in jobs: try: j() except Exception: traceback.print_exc()