| jupyter | ||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
import sys
print(sys.version)
%matplotlib inline
%run -i ~/Python/pjr3compare the CESM SYN (Synthetic) estimate of fields from individual simulations to the CON (Concurrent) syms
Assume they are on lat/lon grid at same resolution**
configured for the GMD paper
def setfig3b1x1 ():
"""
return fig and axes for a single panel figure
"""
plotproj = ccrs.Mollweide()
plotproj._threshold /= 100.
fig, axes = plt.subplots(ncols=1,
gridspec_kw={'width_ratios': [1]},
subplot_kw={'projection': plotproj},
figsize=(6,4.1),
)
fig.set_dpi(300.0)
return fig, axes;
def pltllbox2(xri, yri,ax=None):
if ax is None:
ax = plt.gca()
if xri[1] < xri[0]:
xri[1] += 360.
regcx = [xri[0],xri[1],xri[1],xri[0],xri[0]]
regcy = [yri[0],yri[0],yri[1],yri[1],yri[0]]
ax.plot(regcx,regcy,color='red',transform=ccrs.PlateCarree())# specify filename characteristics by index
## index 1 is concurrent (CONCUR)
## index 2 is control
## index 3 is used for specifying regional file names used to construct the "SYN" variables
# specify regions (assume lon always specified as west, then east limit)
xreg = np.array([[-150.,-110.],[-110,-70],[-25.,15.],[170.,-120.],[-170.,-90.]])%360.
yreg = np.array([[0.,30.], [-30.,0.], [-30.,0.], [30.,50.], [-50.,-30.] ])
namereg = ['NEP','SEP','SEA','NP','SP']
#xreg = [[0.,360.]]
#yreg = [[-90.,91.]]
if True: #CESM CDNC 600 -1.8 W/m2
case_start1 = "/home/jupyter-haruki/work/CESM_MCB/MCB_R1R2R3_CN375cm/MCB_R1R2R3_CN375cm.cam.h0."
case_start1 = "/home/jupyter-haruki/work/CESM_MCB/MCB_R1R2R3_CN600cm/MCB_R1R2R3_CN600cm.cam.h0."
case_start1 = "/e3sm_prod/phil/climo/cesm/MCB_R1R2R3_CN600cm/fv192x288/MCB_R1R2R3_CN600cm.cam.h0.1-10."
case_end1 = ".nc"
pref1='CESM_CN600'
fstring1 ='%s%s%s'
case_start2 = "/home/jupyter-haruki/work/CESM_MCB/Fixed_SST/Fixed_SST.cam.h0."
case_start2 = "/e3sm_prod/phil/climo/cesm/Fixed_SST/fv192x288/Fixed_SST.cam.h0.1-20."
case_end2 = ".nc"
pref2='Control'
fstring2 ='%s%s%s'
if False: #E3SM CDNC600 very weak forcing
case_start1 = "/e3sm_prod/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_C600_remapped/20221018.v2.LR.F2010.E1_R1-3_CDNC600.eam.h0.y1-5.FORCING.nc"
case_end1 = ""
pref1='E3SM_CN600'
fstring1 ='%s%.0s%.0s'
case_start2 = "/e3sm_prod/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_CNTL_remapped/20220930.v2.LR.F2010.E1_CNTL.eam.h0.y1-14.FORCING.nc"
case_end2 = ""
fstring2 ='%s%.0s%.0s'
pref2='E3SMcontrol'
if False: # E3SM CDNC2000 -1.8W/m2 CI
#case_start1 = "/scratch2/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_C600_remapped/20221018.v2.LR.F2010.E1_R1-3_CDNC600.eam.h0.1-11"
case_start1 = "/e3sm_prod/phil/climo/e3sm/20221123.v2.LR.F2010.E1_R1-3_CDNC2000/fv192x288/20221123.v2.LR.F2010.E1_R1-3_CDNC2000.eam.h0.1-11."
case_end1 = ".nc"
pref1='E3SM_CN2000'
fstring1 ='%s%.0s%.0s'
fstring1 ='%s%s%s'
#case_start2 = "/scratch2/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_CNTL_remapped/20220930.v2.LR.F2010.E1_CNTL.eam.h0.1-14"
case_start2 = "/e3sm_prod/phil/climo/e3sm/20220930.v2.LR.F2010.E1_CNTL/fv192x288/20220930.v2.LR.F2010.E1_CNTL.eam.h0.1-14."
case_end2 = ".nc"
fstring2 ='%s%.0s%.0s'
fstring2 ='%s%s%s'
pref2='E3SMcontrol'
if False: # E3SM 50Tg/yr
#case_start1 = "/scratch2/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_C600_remapped/20221018.v2.LR.F2010.E1_R1-3_CDNC600.eam.h0.1-11"
case_start1 = "/e3sm_prod/phil/climo/e3sm/20230426.v2.LR.F2010.MCB-SSLT-EM.R1-3.test01/20230426.v2.LR.F2010.MCB-SSLT-EM.R1-3.test01_ANN_000101_001112_climo_fv192x288.nc"
case_end1 = ""
pref1='E3SM_50Tgpy'
fstring1 ='%s%.0s%.0s'
#fstring1 ='%s%s%s'
#case_start2 = "/scratch2/PJR/haruki_workdir/E3SM_MCB/F2010.E1_R1-3_CNTL_remapped/20220930.v2.LR.F2010.E1_CNTL.eam.h0.1-14"
case_start2 = "/e3sm_prod/phil/climo/e3sm/20220930.v2.LR.F2010.E1_CNTL/fv192x288/20220930.v2.LR.F2010.E1_CNTL.eam.h0.1-14."
case_end2 = ".nc"
fstring2 ='%s%.0s%.0s'
fstring2 ='%s%s%s'
pref2='E3SMcontrol'
if False: # CESM 7.5 Tg/yr
case_start1 = "/e3sm_prod/phil/climo/cesm/F2010climo.ss_NEP_SEP_SEA.1.5Tg/fv192x288/F2010climo.ss_NEP_SEP_SEA.1.5Tg.cam.h0."
case_end1 = ".1-25.nc"
pref1='CESM_7.5Tgpyr'
fstring1 ='%s%.0s%.0s'
fstring1 ='%s%s%s'
case_start2 = "/e3sm_prod/phil/climo/cesm/Fixed_SST/fv192x288/Fixed_SST.cam.h0.1-20."
case_end2 = ".nc"
fstring2 ='%s%.0s%.0s'
fstring2 ='%s%s%s'
pref2='CESMcontrol'
Varname='<Varname>'
ind1 = fstring1 % (case_start1,Varname,case_end1)
print('example string used for file open',ind1)
def xr_llhplot2 (xrVar, cbar='default', plotproj=None, ax=None, cax=None,
ylabels=None, clevs=None, cmap=None, title=None, cbartitle=None):
"""xr_llhplot xarray lat lon horizontal plot
"""
#print(' entering xr_llhplot', xrVar)
lon=xrVar['lon'].values
lat=xrVar['lat'].values
xv,yv=np.meshgrid(lon,lat)
data_regridded = xrVar.values
#print('aaa',data_regridded.shape, xv.shape, yv.shape)
df = data_regridded.flatten()
dsub = df[np.isfinite(df)] # ignore NaN
zmax = dsub.max()
zmin = dsub.min()
#print('masked interpolated range',zmin,zmax)
dataproj=ccrs.PlateCarree() # data is always assumed to be lat/lon
if ylabels is None: ylabels = True
if clevs is None:
clevs = findNiceContours(np.array([zmin,zmax]),nlevs=10)
#print('clevs',clevs)
if cmap is None:
#print('aaa, grabbing cmap default')
#cmap = mpl.cm.get_cmap()
cmap = plt.get_cmap()
#print('bbb',cmap.N)
#print('cmap',cmap)
extend = 'both'
norm = mpl.colors.BoundaryNorm(clevs,cmap.N,extend=extend)
#print('norm',norm(clevs))
clat = (lat.min()+lat.max())/2.
clon = (lon.min()+lon.max())/2.
if plotproj is None:
plotproj = ccrs.PlateCarree()
plotproj = ccrs.Mollweide()
# if no ax argument, could get current axis, or create it
if ax is None:
#print('grab current axis')
#ax = plt.gca()
ax = plt.axes(projection=plotproj)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=2, color='gray', alpha=0.5)
pl = ax.contourf(xv, yv, data_regridded, levels=clevs, # vmin=zmin, vmax=zmax,
norm=norm, cmap=cmap,
extend=extend, transform=ccrs.PlateCarree())
gl.left_labels=ylabels
gl.right_labels=ylabels
ax.coastlines(linewidth=1,color='blue')
## Find the location of the main plot axes
## has to be done after some plotting is done in projection space
posn = ax.get_position()
# print some registration marks to help in lining up figures
ax2 = fig.add_axes([0,0,0.1,0.1])
ax2.set_position([posn.x0-0.005, posn.y0-0.005, posn.width+0.01, posn.height+0.01])
ax2.patch.set_alpha(0.0)
ax2.scatter([0,0,1,1], [0,1,0,1], c="r", s=100)
ax2.set_axis_off()
ax2.set_xlim([0,1])
ax2.set_ylim([0,1])
if not title is None:
ax.set_title(title)
# Add colorbar to plot
if cbartitle is None:
cbartitle = xrVar.long_name
if cbar == 'default':
if cax is not None:
cax = ax
else:
# create an colorbar axis
cax = fig.add_axes([0,0,0.1,0.1])
## Adjust the positioning and orientation of the colorbar
cax.set_position([posn.x0, posn.y0-0.06, posn.width, 0.04])
cb = plt.colorbar(
pl, orientation='horizontal',ticks=clevs,cax=cax,
label='%s (%s)'%(cbartitle, xrVar.units)
)
cb.ax.tick_params(labelsize=8)
return# plt.style.use('fivethirtyeight')
regtag = ""
weights = None
Varlist = np.array(['RESTOM','FLNT','FSNT','TS','TMQ','PRECT','AEROD_v','CLDLOW','CLDTOT','LWCF','SWCF','TGCLDIWP','TGCLDLWP',
'SHFLX','LHFLX','PBLH','PCONVT','PRECC','PRECS'])
Varlist = np.array(['FLNT','FSNT','TS','PRECC','PRECL','CLDLOW','CLDTOT','LWCF','SWCF','TGCLDIWP','TGCLDLWP',
'SHFLX','LHFLX','PBLH','PRECSC','PRECSL'])
#Varlist = np.array(['TS','TMQ','PRECT'])
#Varlist = np.array(['RESTOM','LWCF','SWCF','FLNT','FSNT'])
#Varlist = np.array(['AEROD_v'])
Varlist = np.array(['FSNT','FSNTC','FLNT','FLNTC'])
Varlist = np.array(['SWCF','TGCLDLWP','CLDLOW'])
#Varlist = np.array(['FSUTOA','FSNT','FSNTC'])
Varlist = np.array(['PRECC'])
Varlist = np.array(['CLDTOT','TGCLDLWP','FSNT','FSNTC','AODVIS'])
Varlist = np.array(['CLDTOT'])
Varlist = np.array(['TGCLDLWP'])
Varlist = np.array(['SWCF'])
#Varlist = np.array(['SOLIN'])
for Varname in Varlist:
print()
print('-------------------------------'+Varname)
ind1 = fstring1 % (case_start1,Varname,case_end1)
print('opening',ind1)
DS1 = xr.open_mfdataset(ind1)
#print('xxx',DS1.time)
DS1 = center_time(DS1)
# define a reference DataSet that we reconcile all datasets to
DSR = DS1.copy()
ind2 = fstring2 % (case_start2,Varname,case_end2)
print('opening ind2',ind2)
#DS2 = xr.open_mfdataset(ind2)
DS2 = xr.open_mfdataset(ind2)
DS2 = center_time(DS2)
DSR, DS1 = reconcile_xr_coords(DSR, DS1)
DSR, DS2 = reconcile_xr_coords(DSR, DS2)
Var1 = xr_getvar(Varname,DS1)
V1 = Var1.mean(dim='time',keep_attrs=True)
Var2 = xr_getvar(Varname,DS2)
V2 = Var2.mean(dim='time',keep_attrs=True)
DV = V1-V2
print('DV range', DV.min().values, DV.max().values)
if 'area' in DS1:
area = DS1['area']
elif 'area' in DS2:
area = DS2['area']
else:
print('calculating weights')
lat = Var1['lat'].values
lon = Var1['lon'].values
area = make_fvarea(lon,lat)
weights = V1.copy()
weights.data =area
weights.attrs['units']='steradians'
print(Varname, V1.attrs['long_name'],'Range V1 and V2 ',V1.min().values, V1.max().values, V2.min().values, V2.max().values)
V1A = V1.weighted(weights).mean()
sV1A = ' (%5.2f)' % V1A
V2A = V2.weighted(weights).mean()
sV2A = ' (%5.2f)' % V2A
DVA = V1A-V2A
sDVA = ' (%5.2f)' % DVA
print('area avgs '+pref1+' %5.2f' % (V1A.values),' '+pref2+' %5.2f' % (V2A.values),' Delta %5.2f' % (DVA.values))
if V1.min().values == V1.max().values:
print('constant field skipping plot ')
else:
clev_rng = {'CDNUMC':np.array([0.,3.e11]),'FSNT':np.array([40.,360]),
'TGCLDLWP':np.array([0.,280.]),'PRECL':np.array([0.,10.]),
'PRECC':np.array([0.,16.]),'SWCF':np.array([-140.,0.]),
'CLDLOW':np.array([0.,90.]),'AODVIS':np.array([0,0.2]),
'XXX':np.array([-45.,45.]),
}
dlev_rng = {'CDNUMC':np.array([0.,3.e11])/2.,'FSNT':np.array([-45.,45.]),
'TGCLDLWP':np.array([-50.,50.]),'PRECL':np.array([-1.,1.]),
'PRECC':np.array([-1.,1.]),'SWCF':np.array([-45.,45.]),
'CLDLOW':np.array([-10.,10.]),'AODVIS':np.array([-.9,.9]),
'CLDTOT':np.array([-10.,10.]),
}
if Varname in clev_rng:
clevs = findNiceContours(clev_rng[Varname],nlevs = 12)
else:
clevs = findNiceContours(np.array([V1.values,V2.values]),nlevs = 12)
if Varname in dlev_rng:
dlevs = findNiceContours(dlev_rng[Varname],nlevs = 15,rmClev=0.,sym=True)
else:
dlevs = findNiceContours(np.array([DV.min().values,DV.max().values]),nlevs = 15, rmClev=0.,sym=True)
#dlevs = [-5.,-2.,-1.,-0.5,-0.2,-0.1,0.1,0.2,0.5,1.,2.,5.]
#print('xxx',dlevs)
dmap = diverge_map()
plconf = '3-1x1'
# good setup for 3 rows of 1 columns
if plconf == '3-1x1':
fig, axes = setfig3b1x1()
xr_llhplot2(DV, ax=axes,clevs=dlevs,cmap=dmap,title=pref1+'-'+pref2+sDVA, ylabels=False)#,cbar=None)
pltllbox2([-150.,-110.],[0.,30.],ax=axes)
pltllbox2([-110.,-70.],[-30.,0.],ax=axes)
pltllbox2([-25.,15.],[-30.,0.],ax=axes)
#plt.rcParams['savefig.facecolor'] = "0.8"
#axes.patch.set_facecolor('xkcd:mint green')
#axes.patch.set_facecolor('xkcd:red')
#fig.patch.set_facecolor('xkcd:mint green') # This changes the grey to white
#plt.tight_layout()
plt.savefig(pref1+'_'+Varname+'-D.pdf',format='pdf',dpi=300,transparent=True)#,facecolor='xkcd:mint green')
plt.show()
print('field processing complete')
SWCF_NCCLIMO_RIGHT = V1casename0 = "Fixed_SST" # reference run
#casename0 = "MCB_R1R2R3_CN600cm"
casename1 = "MCB_R1R2R3_CN600cm_NEP"
casename2 = "MCB_R1R2R3_CN600cm_SEP"
casename3 = "MCB_R1R2R3_CN600cm_SEA"
case_end = ".nc"
pref1='SYN CESM_CN600'
fstring ='%s%s%s'
pref2 = 'Control'
Varname='FSNT'
Varname='SWCF'
#Varname ='SOLIN'
def bld_fname(casename, Varname):
fname = "/e3sm_prod/phil/timeseries/cesm2-mcb-reshaped/"+casename+"/"+casename+".cam.h0.1-10."+Varname+".nc"
return fname
def bld_fname2(casename, Varname):
fname = "/e3sm_prod/phil/timeseries/cesm2-mcb-reshaped/"+casename+"/"+casename+".cam.h0."+Varname+".y1-10.nc"
return fname
def bld_fname3(casename, Varname):
fname = "/e3sm_prod/phil/timeseries/cesm2-mcb-reshaped/"+casename+"/"+casename+".cam.h0.1-20."+Varname+".nc"
return fname
ind0 = bld_fname3(casename0, Varname)
print('ind0',ind0)
ind1 = bld_fname(casename1, Varname)
ind2 = bld_fname(casename2, Varname)
ind3 = bld_fname(casename3, Varname)
DS0 = center_time(xr.open_mfdataset(ind0))
DS1 = center_time(xr.open_mfdataset(ind1))
DS2 = center_time(xr.open_mfdataset(ind2))
DS3 = center_time(xr.open_mfdataset(ind3))
DSR, DS0 = reconcile_xr_coords(DSR, DS0)
DSR, DS1 = reconcile_xr_coords(DSR, DS1)
DS0, DS2 = reconcile_xr_coords(DS0, DS2)
DS0, DS3 = reconcile_xr_coords(DS0, DS3)
# grab only the first 10 years of the control
#DS0 = DS0.isel(time=slice(0,120))
C0 = DS0[Varname].mean(dim="time")
C1 = DS1[Varname].mean(dim="time")
C2 = DS2[Varname].mean(dim="time")
C3 = DS3[Varname].mean(dim="time")
#dlevs=None
#clevs=None
C0A = C0.weighted(weights).mean()
print("C0A",C0A.values)
# produce the synthetic estimate of the forcing, which is the sum of the forcings minus 2*control
CumS = (C1+C2+C3) - 2*C0
CumSA = CumS.weighted(weights).mean()
CSAD = CumSA - C0A
print("CSAD,C0A",CSAD.values,C0A.values)
C1A = C1.weighted(weights).mean()
C1AD = C1A - C0A
#print('C1A,V2A',C1A.values,V2A.values)
#print(C1AD)
sC1AD = ' (%5.2f)' % C1AD
if False:
fig, axes = setfig3b1x1()
xr_llhplot2(C1-C0, ax=axes,clevs=dlevs,cmap=dmap,title=pref1+'-'+pref2+sC1AD)
pltllbox2([-150.,-110.],[0.,30.],ax=axes)
pltllbox2([-110.,-70.],[-30.,0.],ax=axes)
pltllbox2([-25.,15.],[-30.,0.],ax=axes)
C2A = C2.weighted(weights).mean()
C2AD = C2A - C0A
#print(C2AD)
sC2AD = ' (%5.2f)' % C2AD
if False:
fig, axes = setfig3b1x1()
xr_llhplot2(C2-C0, ax=axes,clevs=dlevs,cmap=dmap,title=pref1+'-'+pref2+sC2AD)
pltllbox2([-150.,-110.],[0.,30.],ax=axes)
pltllbox2([-110.,-70.],[-30.,0.],ax=axes)
pltllbox2([-25.,15.],[-30.,0.],ax=axes)
C3A = C3.weighted(weights).mean()
C3AD = C3A - C0A
#print(C3AD)
sC3AD = ' (%5.2f)' % C3AD
if False:
fig, axes = setfig3b1x1()
xr_llhplot2(C3-C0, ax=axes,clevs=dlevs,cmap=dmap,title=pref1+'-'+pref2+sC3AD)
pltllbox2([-150.,-110.],[0.,30.],ax=axes)
pltllbox2([-110.,-70.],[-30.,0.],ax=axes)
pltllbox2([-25.,15.],[-30.,0.],ax=axes)
sCSAD = ' (%5.2f)' % CSAD
fig, axes = setfig3b1x1()
xr_llhplot2(CumS-C0, ax=axes,clevs=dlevs,cmap=dmap,title=pref1+'-'+pref2+sCSAD)
pltllbox2([-150.,-110.],[0.,30.],ax=axes)
pltllbox2([-110.,-70.],[-30.,0.],ax=axes)
pltllbox2([-25.,15.],[-30.,0.],ax=axes)
plt.savefig(pref1+'_'+Varname+'-D.pdf',format='pdf',dpi=300,transparent=True)#,facecolor='xkcd:mint green')