import netCDF4,glob,numpy as np,matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt, cartopy.crs as ccrs, cartopy.feature as cf
fs=sorted(glob.glob('data/2026-09-29/2026-09-29/*.nc'))
H=[15,20,25,30]; out={h:[] for h in H}; n=0
for f in fs:
    d=netCDF4.Dataset(f); d.set_auto_mask(False)
    if str(d.getncattr('bad')).strip()!='0': d.close(); continue
    n+=1
    z=d['alt_refrac'][0]/1e3; t=d['dry_temp'][0]
    la=d['lat_tp'][0]; lo=d['lon_tp'][0]; ihb=(d['impact_opt'][0]-d['roc'][0]-d['undulation'][0])/1e3
    lat0,lon0=float(d['lat'][0]),float(d['lon'][0])
    d.close()
    v=(z>-9e3)&(t>100)&(t<400)
    if v.sum()<2: continue
    z,t=z[v],t[v]; o=np.argsort(z); z,t=z[o],t[o]
    # tangent-point location: from lev1b grid (impact height ~ MSL height of tangent point, offset by refraction; small for >15 km)
    w=(la>-900)&(lo>-900)&(ihb>-9e3)
    if w.sum()>=2:
        ih=ihb[w]; oo=np.argsort(ih); ih=ih[oo]; lab=la[w][oo]; lob=np.degrees(np.unwrap(np.radians(lo[w][oo])))
    for h in H:
        if z[0]<=h<=z[-1]:
            T=np.interp(h,z,t)
            if w.sum()>=2: y,x=np.interp(h,ih,lab),(np.interp(h,ih,lob)+180)%360-180
            else: y,x=lat0,lon0
            out[h].append((y,x,T))
print('profiles',n,{h:len(out[h]) for h in H})
for h in H:
    a=np.array(out[h]); print(h,a[:,2].min(),np.percentile(a[:,2],[1,50,99]),a[:,2].max())
fig,axs=plt.subplots(2,2,figsize=(16,9.5),subplot_kw=dict(projection=ccrs.Robinson()))
for ax,h in zip(axs.flat,H):
    ax.set_global(); ax.add_feature(cf.LAND.with_scale('50m'),facecolor='#e8e8e8',zorder=0); ax.coastlines('50m',lw=0.5,color='0.3')
    gl=ax.gridlines(draw_labels=True,lw=0.3,color='0.6',alpha=0.6); gl.top_labels=gl.right_labels=False; gl.xlabel_style=gl.ylabel_style={'size':8}
    a=np.array(out[h]); vmin,vmax=np.percentile(a[:,2],[1,99])
    sc=ax.scatter(a[:,1],a[:,0],c=a[:,2],s=7,cmap='RdYlBu_r',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree(),edgecolor='k',lw=0.1)
    cb=plt.colorbar(sc,ax=ax,orientation='horizontal',pad=0.07,fraction=0.045,aspect=40,extend='both'); cb.set_label('Dry temperature (K)',fontsize=9)
    ax.set_title(f'{h} km MSL  (n = {len(a)})',fontsize=11)
fig.suptitle('YM08 RO dry temperature at fixed MSL heights (bad = 0), 2026-09-29',fontsize=14)
plt.savefig('/mnt/user-data/outputs/dry_temp_msl_heights.png',dpi=150,bbox_inches='tight')