← LAMPIRAN MENULAMPIRAN 14 / 18
S1D4NG.5KR1P51.404 · ANNEX FILE
Lampiran 14. Program python CAPE dan CIN 23-25 Maret 2025 di Kabupaten Kudus
//⚠ 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 14. Program python CAPE dan CIN 23-25 Maret 2025 di Kabupaten Kudus
| import cfgrib import xarray as xr import numpy as np import matplotlib.pyplot as plt import matplotlib.dates as mdates import matplotlib.patches as mpatches from matplotlib.path import Path as MplPath from pathlib import Path import pandas as pd import geopandas as gpd from scipy.interpolate import make_interp_spline from shapely.geometry import MultiPolygon, Polygon from shapely.ops import unary_union OUTPUT_DIR = Path(r"~\Skripsi\intisari\Cape Cin\CC4") KUDUS_GRIB = r"~\Skripsi\era\single level\Kudus Single level.grib" JAWA_GRIB = r"~\Skripsi\era\single level\Jawa single.grib" KABUPATEN_GEOJSON = r"~\Skripsi\GeoJson-Indonesia-38-Provinsi\Kabupaten\38 Provinsi Indonesia - Kabupaten.json" KEC_KUDUS_GEOJSON = r"~\Skripsi\indonesia-district\id33_jawa_tengah\id3319_kudus\id3319_kudus.geojson" PERIODS = { "Sebelum": { "label": "Sebelum Kejadian Banjir\n(22 Maret 2025, 14.00\u201317.30 WIB)", "start_utc": "2025-03-22T07:00", "end_utc": "2025-03-22T10:00", }, "Saat": { "label": "Saat Kejadian Banjir\n(24\u201325 Maret 2025)", "start_utc": "2025-03-24T16:00", "end_utc": "2025-03-25T03:00", }, "Pasca": { "label": "Pasca Kejadian Banjir\n(28 Maret 2025)", "start_utc": "2025-03-28T10:00", "end_utc": "2025-03-28T17:00", }, } JAWA_LON = (109.20, 112.35) JAWA_LAT = (-7.80, -5.70) KUDUS_LON = (110.72, 111.02) KUDUS_LAT = (-7.00, -6.58) print("Loading shapefiles...") gdf_kab = gpd.read_file(KABUPATEN_GEOJSON) JAWA_PROVS = [ "Jawa Tengah", "Jawa Timur", "Jawa Barat", "DKI Jakarta", "Daerah Istimewa Yogyakarta", "Banten" ] gdf_kab_jawa = gdf_kab[gdf_kab["WADMPR"].isin(JAWA_PROVS)].copy() gdf_villages = gpd.read_file(KEC_KUDUS_GEOJSON) gdf_kec_kudus = gdf_villages.dissolve(by="district", aggfunc="first").reset_index() # Gabungan polygon Kudus untuk clip kudus_polygon = unary_union(gdf_kec_kudus.geometry) print(f" Kabupaten Jawa: {len(gdf_kab_jawa)}") print(f" Kecamatan Kudus: {len(gdf_kec_kudus)}") print("Loading GRIB data...") ds_jawa_cape = xr.open_dataset(JAWA_GRIB, engine="cfgrib", backend_kwargs={"filter_by_keys": {"shortName": "cape"}}) ds_kudus_cape = xr.open_dataset(KUDUS_GRIB, engine="cfgrib", backend_kwargs={"filter_by_keys": {"shortName": "cape"}}) ds_kudus_cin = xr.open_dataset(KUDUS_GRIB, engine="cfgrib", backend_kwargs={"filter_by_keys": {"shortName": "cin"}}) ds_kudus_cin = ds_kudus_cin.sel(step=pd.Timedelta(hours=1)) ds_kudus_cin = ds_kudus_cin.assign_coords(time=ds_kudus_cin["valid_time"]) ds_kudus_cin = ds_kudus_cin.drop_vars(["valid_time", "step"]) print("Cropping data to domain...") cape_jawa_dom = ds_jawa_cape["cape"].sel( latitude=slice(JAWA_LAT[1], JAWA_LAT[0]), longitude=slice(JAWA_LON[0], JAWA_LON[1]) ) cape_kudus_dom = ds_kudus_cape["cape"].sel( latitude=slice(KUDUS_LAT[1], KUDUS_LAT[0]), longitude=slice(KUDUS_LON[0], KUDUS_LON[1]) ) cin_kudus_dom = ds_kudus_cin["cin"].sel( latitude=slice(KUDUS_LAT[1], KUDUS_LAT[0]), longitude=slice(KUDUS_LON[0], KUDUS_LON[1]) ) print(f" Jawa grid: {cape_jawa_dom.shape}") print(f" Kudus CAPE grid: {cape_kudus_dom.shape}") print(f" Kudus CIN grid: {cin_kudus_dom.shape}") def get_period_mean(data_array, start, end): subset = data_array.sel(time=slice(start, end)) if subset.sizes.get("time", 0) == 0: mid = pd.Timestamp(start) + (pd.Timestamp(end) - pd.Timestamp(start)) / 2 return data_array.sel(time=mid, method="nearest") return subset.mean(dim="time") def mask_data_to_polygon(lon, lat, data, polygon): """Mask data di luar polygon, return NaN di luar polygon.""" from shapely.geometry import mapping lon2d, lat2d = np.meshgrid(lon, lat) masked = data.copy() for i in range(len(lat)): for j in range(len(lon)): from shapely.geometry import Point if not polygon.contains(Point(lon2d[i, j], lat2d[i, j])): masked[i, j] = np.nan return masked def smooth_contourf(lon, lat, data, levels, cmap, vmin, vmax): """Buat contourf dengan interpolasi untuk tampilan lebih halus.""" # Interpolasi ke grid lebih halus from scipy.interpolate import griddata lon_fine = np.linspace(lon.min(), lon.max(), 100) lat_fine = np.linspace(lat.min(), lat.max(), 80) lon2d, lat2d = np.meshgrid(lon, lat) lon2d_f, lat2d_f = np.meshgrid(lon_fine, lat_fine) data_fine = griddata((lon2d.ravel(), lat2d.ravel()), data.ravel(), (lon2d_f, lat2d_f), method="cubic") return lon_fine, lat_fine, data_fine def plot_jawa_kudus(ax, data, lon, lat, title, cmap, vmin, vmax, gdf_admin, clip_polygon=None): """Plot peta Jawa dengan data ter-clip ke polygon Kudus.""" # Smooth lon_f, lat_f, data_f = smooth_contourf(lon, lat, data, None, cmap, vmin, vmax) # Mask ke polygon Kudus if clip_polygon is not None: data_f = mask_data_to_polygon(lon_f, lat_f, data_f, clip_polygon) levels = np.linspace(vmin, vmax, 20) cf = ax.contourf(lon_f, lat_f, data_f, levels=levels, cmap=cmap, vmin=vmin, vmax=vmax, extend="both") # Batas admin gdf_admin.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=0.35, zorder=5) # Garis polygon Kudus gpd.GeoDataFrame(geometry=[clip_polygon]).plot(ax=ax, facecolor="none", edgecolor="red", linewidth=1.5, linestyle="--", zorder=6) ax.set_xlim(JAWA_LON[0], JAWA_LON[1]) ax.set_ylim(JAWA_LAT[0], JAWA_LAT[1]) ax.set_aspect("auto") ax.set_title(title, fontsize=10, fontweight="bold", pad=6) ax.set_xlabel("Longitude (\u00b0BT)", fontsize=9) ax.set_ylabel("Latitude (\u00b0LS)", fontsize=9) ax.tick_params(labelsize=8) return cf print("\nGenerating temporal plot...") # Rata-rata spasial: pakai domain yang SAMA (area Kudus) untuk Jawa & Kudus cape_jawa_ts = cape_jawa_dom.mean(dim=["latitude", "longitude"]) cape_kudus_ts = cape_kudus_dom.mean(dim=["latitude", "longitude"]) cin_kudus_ts = cin_kudus_dom.mean(dim=["latitude", "longitude"]) cape_time = cape_jawa_ts.time.values cape_jawa_vals = cape_jawa_ts.values cape_kudus_vals = cape_kudus_ts.values cin_time = cin_kudus_ts.time.values cin_vals = cin_kudus_ts.values valid_mask = ~np.isnan(cin_vals) cin_time_valid = cin_time[valid_mask] cin_vals_valid = cin_vals[valid_mask] cin_time_num = mdates.date2num(pd.to_datetime(cin_time_valid)) cin_time_smooth = np.linspace(cin_time_num[0], cin_time_num[-1], 200) k_order = min(3, len(cin_time_valid) - 1) cin_spline = make_interp_spline(cin_time_num, cin_vals_valid, k=k_order) cin_smooth_vals = cin_spline(cin_time_smooth) cin_smooth_time = mdates.num2date(cin_time_smooth) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 8), sharex=True, gridspec_kw={"height_ratios": [1, 1]}) ax1.plot(cape_time, cape_jawa_vals, color="orangered", linewidth=1.5, label="CAPE Jawa (area Kudus)", alpha=0.9) ax1.plot(cape_time, cape_kudus_vals, color="darkred", linewidth=1.5, linestyle="--", label="CAPE Kudus", alpha=0.9) ax1.fill_between(cape_time, cape_jawa_vals, alpha=0.15, color="orangered") ax1.set_ylabel("CAPE (J/kg)", fontsize=11) ax1.set_title("Variasi Temporal CAPE dan CIN\n22\u201328 Maret 2025", fontsize=13, fontweight="bold") ax1.grid(True, alpha=0.3, linestyle="--") flood_times = { "Sebelum": pd.Timestamp("2025-03-22T07:00"), "Saat": pd.Timestamp("2025-03-24T16:00"), "Pasca": pd.Timestamp("2025-03-28T10:00"), } colors_flood = ["green", "red", "blue"] labels_flood = ["Sebelum", "Saat Banjir", "Pasca"] for (name, t), c, lbl in zip(flood_times.items(), colors_flood, labels_flood): ax1.axvline(x=t, color=c, linestyle=":", linewidth=1.8, alpha=0.7, label=lbl) ax1.legend(loc="upper right", fontsize=8, framealpha=0.9) ax2.plot(cin_smooth_time, cin_smooth_vals, color="steelblue", linewidth=1.8, label="CIN Kudus", alpha=0.9) ax2.scatter(pd.to_datetime(cin_time_valid), cin_vals_valid, color="navy", s=20, zorder=5, label="Data asli (step 1h)", alpha=0.7) ax2.fill_between(cin_smooth_time, cin_smooth_vals, alpha=0.15, color="steelblue") ax2.set_ylabel("CIN (J/kg)", fontsize=11) ax2.set_xlabel("Waktu (UTC)", fontsize=11) ax2.grid(True, alpha=0.3, linestyle="--") for (name, t), c in zip(flood_times.items(), colors_flood): ax2.axvline(x=t, color=c, linestyle=":", linewidth=1.8, alpha=0.7) ax2.legend(loc="upper right", fontsize=8, framealpha=0.9) ax2.xaxis.set_major_formatter(mdates.DateFormatter("%d/%m")) ax2.xaxis.set_major_locator(mdates.DayLocator()) plt.xlim(cape_time[0], cape_time[-1]) plt.tight_layout() fname = "Temporal_CAPE_CIN.png" fig.savefig(OUTPUT_DIR / fname, dpi=200, bbox_inches="tight") plt.close(fig) print(f" Saved: {fname}") import shutil script_src = Path(__file__).resolve() script_dst = OUTPUT_DIR / "cape_cin_analysis.py" shutil.copy2(script_src, script_dst) print(f" Script copied: {script_dst}") |