← LAMPIRAN MENULAMPIRAN 15 / 18
S1D4NG.5KR1P51.404 · ANNEX FILE
Lampiran 15. Program python perkembangan awan konvektif
//⚠ S1D4NG.5KR1P51.404 — BIOHAZARD ZONE·//⚠ NO SHARP · NO NARCOTICS · NO WEIRD CARGO·//⚠ THE HUM IS LOUDER WHEN YOU NOTICE IT·//⚠ ERROR 404 RECOVERED — STILL WRONG·//⚠ LEVEL 0 · DO NOT TRUST THE CORRIDOR·
Lampiran 15. Program python perkembangan awan konvektif
| import bz2 from pathlib import Path import matplotlib.pyplot as plt import matplotlib.colors as mcolors import numpy as np import geopandas as gpd HIMA24_DIR = Path(r"~\Himawari\Data TIR IR13 Himawari\HIMA 24") HIMA25_DIR = Path(r"~\Himawari\Data TIR IR13 Himawari\hima 25") GEOJSON_JATENG = Path(r"~\Himawari\indonesia-district\id33_jawa_tengah\id33_jawa_tengah_district.geojson") GEOJSON_JATIM = Path(r"~\Himawari\indonesia-district\id35_jawa_timur\id35_jawa_timur_district.geojson") OUT_DIR = Path(r"~\Skripsi\output") GRID_SHAPE = (6000, 6000) RAW_DTYPE = ">u2" FULL_LON = (80, 200) FULL_LAT = (60, -60) DOMAIN = dict(lon_min=110.1, lon_max=111.8, lat_min=-7.5, lat_max=-5.9) KUDUS = (110.87, -6.79) VMIN, VMAX = -90, 40 PANELS_6 = [ ("202503241400", "21:00 WIB"), ("202503241700", "00:00 WIB"), ("202503242000", "03:00 WIB"), ("202503242300", "06:00 WIB"), ("202503250200", "09:00 WIB"), ("202503250400", "11:00 WIB"), ] def bmkg_cmap_25white(): colors_list = [ (-90, "#000000"), (-80, "#1a1a4e"), (-70, "#0000cc"), (-60, "#0066ff"), (-50, "#00ccff"), (-40, "#00ffcc"), (-30, "#00ff66"), (-20, "#66ff00"), (-10, "#ccff00"), ( 0, "#ffff00"), ( 10, "#ff8800"), ( 20, "#ff4400"), ( 25, "#ffffff"), ( 40, "#ffffff"), ] vals = [c[0] for c in colors_list] hexs = [c[1] for c in colors_list] rgbs = [mcolors.to_rgb(h) for h in hexs] norm_vals = [(v - VMIN) / (VMAX - VMIN) for v in vals] cmap = mcolors.LinearSegmentedColormap.from_list( "bmkg_25white", list(zip(norm_vals, rgbs)), N=256 ) cmap.set_over("#ffffff") cmap.set_under("#000000") return cmap def find_file(ts): for d in [HIMA24_DIR, HIMA25_DIR]: for f in d.glob("*.bz2"): if ts in f.name: return f return None def read_tir(path): with bz2.open(path, "rb") as f: raw = f.read() arr = np.frombuffer(raw, dtype=RAW_DTYPE).reshape(GRID_SHAPE).astype(np.float32) p50 = np.nanpercentile(arr, 50) scale = 0.1 if 1000 <= p50 <= 4000 else (0.01 if 10000 <= p50 <= 40000 else 1.0) return arr * scale - 273.15 def crop_domain(tbb2d, domain): lon = np.linspace(FULL_LON[0], FULL_LON[1], tbb2d.shape[1]) lat = np.linspace(FULL_LAT[0], FULL_LAT[1], tbb2d.shape[0]) lon_idx = np.where((lon >= domain["lon_min"]) & (lon <= domain["lon_max"]))[0] lat_idx = np.where((lat >= domain["lat_min"]) & (lat <= domain["lat_max"]))[0] cropped = tbb2d[np.ix_(lat_idx, lon_idx)] extent = [lon[lon_idx[0]], lon[lon_idx[-1]], lat[lat_idx[-1]], lat[lat_idx[0]]] return cropped, extent def load_kabupaten(): gdfs = [] for path in [GEOJSON_JATENG, GEOJSON_JATIM]: gdf = gpd.read_file(path) kab = gdf.dissolve(by="regency").reset_index() gdfs.append(kab) combined = gpd.GeoDataFrame(pd.concat(gdfs, ignore_index=True)) if len(gdfs) > 1 else gdfs[0] return combined def plot_kab_boundaries(ax, kab_gdf): for geom in kab_gdf.geometry: if geom.geom_type == "MultiPolygon": for poly in geom.geoms: plot_polygon(ax, poly) elif geom.geom_type == "Polygon": plot_polygon(ax, geom) def plot_polygon(ax, poly): xs, ys = poly.exterior.xy ax.plot(xs, ys, color="black", linewidth=0.5, zorder=4, solid_capstyle="round") for interior in poly.interiors: ixs, iys = interior.xy ax.plot(ixs, iys, color="black", linewidth=0.4, zorder=4, linestyle="--") def render_panel(ax, ts, wib_label, kab_gdf, kab_in_domain, cmap, show_marker=True): fpath = find_file(ts) if fpath is None: print(f"MISSING: {ts}") return None print(f" {ts} ({wib_label})...") tbb = read_tir(fpath) cropped, extent = crop_domain(tbb, DOMAIN) ax.set_xlim(extent[0], extent[1]) ax.set_ylim(extent[2], extent[3]) ax.set_aspect("equal") im = ax.imshow(cropped, extent=extent, origin="upper", cmap=cmap, vmin=VMIN, vmax=VMAX, interpolation="bilinear", zorder=1) plot_kab_boundaries(ax, kab_gdf) # Label kabupaten for _, row_kab in kab_in_domain.iterrows(): centroid = row_kab.geometry.centroid name = row_kab["regency"] if DOMAIN["lon_min"] <= centroid.x <= DOMAIN["lon_max"] and \ DOMAIN["lat_min"] <= centroid.y <= DOMAIN["lat_max"]: ax.text(centroid.x, centroid.y, name, fontsize=4.5, color="black", fontweight="bold", ha="center", va="center", bbox=dict(boxstyle="round,pad=0.1", facecolor="white", edgecolor="none", alpha=0.6), zorder=6) if show_marker: ax.plot(KUDUS[0], KUDUS[1], "r^", markersize=7, markeredgecolor="black", markeredgewidth=0.7, zorder=10) # Time label below ax.text(0.5, -0.06, wib_label, transform=ax.transAxes, fontsize=10, fontweight="bold", ha="center", va="top") ax.tick_params(labelsize=6) return im import pandas as pd def main(): print("Loading kabupaten boundaries (Jateng + Jatim)...") kab_gdf = load_kabupaten() print(f" Total kabupaten: {len(kab_gdf)}") kab_in_domain = kab_gdf.cx[ DOMAIN["lon_min"]:DOMAIN["lon_max"], DOMAIN["lat_min"]:DOMAIN["lat_max"] ] print(f" In domain: {len(kab_in_domain)}") cmap = bmkg_cmap_25white() print("\nGenerating 2×3 horizontal...") fig3, axes3 = plt.subplots(2, 3, figsize=(20, 12)) im = None for idx, (ts, wib_label) in enumerate(PANELS_6): row, col = divmod(idx, 3) im = render_panel(axes3[row][col], ts, wib_label, kab_gdf, kab_in_domain, cmap) if im: cbar_ax = fig3.add_axes([0.92, 0.12, 0.015, 0.76]) cbar = fig3.colorbar(im, cax=cbar_ax, orientation="vertical", extend="both") cbar.set_label("Cloud Top Temperature (°C)", fontsize=10) cbar.ax.tick_params(labelsize=8) plt.subplots_adjust(left=0.04, right=0.90, top=0.97, bottom=0.04, wspace=0.10, hspace=0.18) out3 = OUT_DIR / "cloud_development_6panel_kudus_v3.png" fig3.savefig(out3, dpi=250, bbox_inches="tight", facecolor="white") print(f"Saved 3×2: {out3}") plt.close() if __name__ == "__main__": main() |