import glob import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt # ----------------------- # SETTINGS # ----------------------- pattern = "/home/taneja/Documents/kanika/work/NEW_2025/SGP_new_local_filter/sgp_doppler_lidar/sgpdlprofwstats4newsC1.c1/*.nc" files = sorted(glob.glob(pattern)) print("Found files:", len(files)) # Filters used in your plot cbh_min = 300.0 cbh_max = 1500.0 upfrac_thr = 0.6 cloudfreq_thr = 0.7 # Output out_png_sigma = f"sigma_cb_hist_cbh{int(cbh_min)}_{int(cbh_max)}_upf{upfrac_thr:.1f}_cf{cloudfreq_thr:.1f}.png" out_csv = f"sigma_cb_hourly_median_cbh{int(cbh_min)}_{int(cbh_max)}_upf{upfrac_thr:.1f}_cf{cloudfreq_thr:.1f}.csv" def clean_fill(arr, da): """Replace _FillValue/missing_value with NaN.""" arr = np.array(arr, dtype=float) fv = da.attrs.get("_FillValue", None) mv = da.attrs.get("missing_value", None) if fv is not None: arr[arr == float(fv)] = np.nan if mv is not None: arr[arr == float(mv)] = np.nan return arr # ----------------------- # READ FILES AND EXTRACT SIGMA AT CLOUD BASE # ----------------------- dfs = [] skipped = 0 for f in files: try: ds = xr.open_dataset(f, decode_times=True) required = [ "time", "height", "dl_cbh", "w_variance", "cbw_up_fraction", "dl_cloud_frequency", ] if not all(v in ds.variables for v in required): skipped += 1 ds.close() continue time = pd.to_datetime(ds["time"].values) cbh = clean_fill(ds["dl_cbh"].values, ds["dl_cbh"]) upf = clean_fill(ds["cbw_up_fraction"].values, ds["cbw_up_fraction"]) cf = clean_fill(ds["dl_cloud_frequency"].values, ds["dl_cloud_frequency"]) # sigma profile = sqrt(w_variance) sigma_prof = np.sqrt(ds["w_variance"]) # Interpolate sigma to cloud-base height cbh_da = xr.DataArray(cbh, dims=["time"], coords={"time": ds["time"]}) sigma_cb = sigma_prof.interp(height=cbh_da).to_numpy().astype(float) ds.close() # Apply same filters as CBW mask = ( np.isfinite(cbh) & np.isfinite(upf) & np.isfinite(cf) & np.isfinite(sigma_cb) & (cbh >= cbh_min) & (cbh <= cbh_max) & (upf > upfrac_thr) & (cf > cloudfreq_thr) ) if not np.any(mask): continue df = pd.DataFrame({ "time": time[mask], "sigma_cb": sigma_cb[mask], "dl_cbh": cbh[mask], "cbw_up_fraction": upf[mask], "dl_cloud_frequency": cf[mask], }).dropna() if not df.empty: dfs.append(df) except Exception as e: skipped += 1 print("[SKIP]", f.split("/")[-1], "->", e) try: ds.close() except Exception: pass print("Skipped files:", skipped) if not dfs: raise SystemExit("No valid sigma samples found after filtering.") df_all = ( pd.concat(dfs, ignore_index=True) .dropna() .sort_values("time") .drop_duplicates("time", keep="first") .reset_index(drop=True) ) print("Native filtered sigma samples:", len(df_all)) print("Native sigma min/median/max:", df_all["sigma_cb"].min(), df_all["sigma_cb"].median(), df_all["sigma_cb"].max()) # ----------------------- # HOURLY MEDIAN SIGMA # ----------------------- sigma_hourly = ( df_all.set_index("time")["sigma_cb"] .resample("1H") .median() .dropna() ) df_hourly = sigma_hourly.rename("sigma_cb_hourly_median").reset_index() df_hourly.to_csv(out_csv, index=False) print("Saved:", out_csv) sig_for_stats = df_hourly["sigma_cb_hourly_median"].to_numpy() sig_for_stats = sig_for_stats[np.isfinite(sig_for_stats)] print("Hourly median sigma samples:", len(sig_for_stats)) print("Hourly sigma P10/P50/P90:", np.percentile(sig_for_stats, [10, 50, 90])) ## ----------------------- ## PLOT: SIGMA HISTOGRAM ## ----------------------- #plt.figure(figsize=(8, 5)) #plt.hist( # sig_for_stats, # bins=60, # density=True, # alpha=0.7, # label=fr"$\sigma_w$ at CBH (hourly median), n={len(sig_for_stats)}" #) #sig_p = np.percentile(sig_for_stats, [10, 50, 90]) #plt.axvline( # sig_p[1], # lw=2, # linestyle="--", # label=fr"Median $\sigma$={sig_p[1]:.2f}" #) ##plt.plot( ## [], [], " ", ## label=fr"P10={sig_p[0]:.2f}, P90={sig_p[2]:.2f}" ##) #plt.xlabel(r"Cloud-base turbulence sigma $\sigma_w=\sqrt{w\_variance}$ [m s$^{-1}$]") #plt.ylabel("Density") ##plt.title("Context: turbulence intensity at cloud base\n(same filters as CBW)") #plt.grid(True, alpha=0.3) #plt.legend() #plt.tight_layout() #plt.savefig(out_png_sigma, dpi=300) #plt.show() #print("Saved plot:", out_png_sigma) from scipy.stats import lognorm from scipy.stats import lognorm import numpy as np import matplotlib.pyplot as plt # ----------------------- # PLOT: LOGNORMAL FIT + MEDIAN LINE + IQR SHADING # ----------------------- sig = sig_for_stats.copy() sig = sig[np.isfinite(sig)] sig = sig[sig > 0] shape, loc, scale = lognorm.fit(sig, floc=0) mu_ln = np.log(scale) sigma_ln = shape median = np.exp(mu_ln) minus_1s = np.exp(mu_ln - sigma_ln) plus_1s = np.exp(mu_ln + sigma_ln) # IQR from fitted lognormal q25 = lognorm.ppf(0.25, shape, loc, scale) q75 = lognorm.ppf(0.75, shape, loc, scale) sig_p = np.percentile(sig, [10, 50, 90]) # x starts from ~0 x = np.linspace(1e-4, np.percentile(sig, 99.7), 800) pdf = lognorm.pdf(x, shape, loc, scale) plt.figure(figsize=(8, 5)) # lognormal curve plt.plot(x, pdf, lw=2, label=fr"Lognormal fit (median={median:.2f}, $-1\sigma$={minus_1s:.2f}, $+1\sigma$={plus_1s:.2f})") # median line only up to curve median_pdf = lognorm.pdf(median, shape, loc, scale) plt.vlines( median, ymin=0, ymax=median_pdf, lw=2, linestyle="--", # label=fr"Median $\sigma$={median:.2f}" ) # IQR shading x_iqr = np.linspace(q25, q75, 300) pdf_iqr = lognorm.pdf(x_iqr, shape, loc, scale) plt.fill_between( x_iqr, 0, pdf_iqr, alpha=0.25, label=fr"IQR: {q25:.2f}–{q75:.2f}" ) # extra legend text #plt.plot([], [], " ", label=fr"$-1\sigma$={minus_1s:.2f}, $+1\sigma$={plus_1s:.2f}") #plt.plot([], [], " ", label=fr"P10={sig_p[0]:.2f}, P90={sig_p[2]:.2f}") plt.xlim(0, np.percentile(sig, 99.7)) plt.ylim(bottom=0) plt.xlabel(r"Cloud-base turbulence sigma $\sigma_w=\sqrt{w\_variance}$ [m s$^{-1}$]") plt.ylabel("Density") plt.title("SGP") plt.grid(True, alpha=0.3) plt.legend() plt.tight_layout() plt.savefig("sigma_cb_lognormal_fit_IQR_cbh300_1500_upf0.6_cf0.7.png", dpi=300) plt.show()