import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import scipy
import cftime
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy
import xesmf as xe
import matplotlib as mpl
plt.rc('font', size = 18)
with xr.open_dataset('/tigress/GEOCLIM/LRGROUP/RECCAP2_coastal_ocean_repo/final_datasets/area_mask_wide_narrow_plusAP.nc') as ds:
area_wide = ds.area_wide # km2
area_narrow = ds.area_narrow ## pulled from area_wide and masked with new arabian peninsula mask
narrow_mask = ds.narrow_mask ### regridded AP with nearest neighbor to preserve more of the mask and better match alizee's published figure
wide_mask = ds.wide_mask
area_glob = xr.zeros_like(area_wide)
area_glob = area_glob + area_wide.mean("longitude")
Only load what data you need; the sea ice, pco2, fco2 multimodel files are 30G each
with xr.open_dataset('seaicefrac_models_timemean.nc') as ds:
seaice_model = ds.siconc
with xr.open_dataset('seaicefrac_obs_timemean.nc') as ds:
seaice_obs = ds.siconc
with xr.open_dataset('pco2_atm_reference.nc') as ds:
pco2_air = ds.pco2_air
pco2_air_trend = ds.pco2_air_trend
with xr.open_dataset('fCO2_models_timemean.nc') as ds:
fco2_model = ds.fco2
with xr.open_dataset('fCO2_obs_timemean.nc') as ds:
fco2_obs = ds.fco2
fco2_obskw = ds.fco2_alt_kw
with xr.open_dataset('pCO2_models_timemean.nc') as ds:
pco2_model = ds.pco2
with xr.open_dataset('pCO2_obs_timemean.nc') as ds:
pco2_obs = ds.pco2
with xr.open_dataset('fn2o_models_timemean.nc') as ds:
fn2o_model = ds.fn2o
with xr.open_dataset('fn2o_obs.nc') as ds:
fn2o_obs = ds.fn2o
with xr.open_dataset('fch4_weber.nc') as ds:
fch4_weber = ds.fch4
with xr.open_dataset('globalmean_models.nc') as ds:
mean_fco2_thinshelf = ds.mean_fco2_thinshelf ### mean values in mol/m2/s
mean_fco2_wideshelf = ds.mean_fco2_wideshelf
mean_fco2_open = ds.mean_fco2_open
mean_pco2_thinshelf = ds.mean_pco2_thinshelf ### mean values in uatm
mean_pco2_wideshelf = ds.mean_pco2_wideshelf
mean_pco2_open = ds.mean_pco2_open
total_area_wide = ds.area_wide ## area in m2
total_area_narrow = ds.area_narrow
total_area_open = ds.area_open
with xr.open_dataset('globalmean_obs.nc') as ds:
mean_fco2obs_thinshelf = ds.mean_fco2_thinshelf
mean_fco2obs_thinshelf_modelkw = ds.mean_fco2_thinshelf_momkw ### mean values in mol/m2/s
mean_fco2obs_wideshelf = ds.mean_fco2_wideshelf
mean_fco2obs_wideshelf_modelkw = ds.mean_fco2_wideshelf_momkw
mean_fco2obs_open = ds.mean_fco2_open
mean_pco2obs_thinshelf = ds.mean_pco2_thinshelf ### mean values in uatm
mean_pco2obs_wideshelf = ds.mean_pco2_wideshelf
mean_pco2obs_open = ds.mean_pco2_open
plt.rc('font', size = 16)
import matplotlib as mpl
import matplotlib.patches as mpatches
mpl.rc('hatch', color='k', linewidth=5)
fig, axes = plt.subplots(1,1, figsize=(10,10), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
# ----------------------------
axes.coastlines()
gl=axes.gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(wide_mask*0+.5).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
(narrow_mask*0+1).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
axes.set_title('Wide, Narrow coastal mask')
axes.set_global()
axes.add_patch(mpatches.Rectangle(xy=[-100,17], width=23, height=30,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[63,0], width=37, height=25,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[-80,-58], width=30, height=28,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[120,25], width=30, height=25,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
fig.tight_layout()
# fig.savefig('global_wide_narrow_masks.pdf', format='pdf', bbox_inches='tight')
fig, axes = plt.subplots(1,1, figsize=(6,6), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=140)})
import cartopy.feature as cfeature
# ----------------------------
axes.coastlines()
gl=axes.gridlines()
axes.add_feature(cfeature.LAND,facecolor='lightgrey')
# (wide_mask*0+.5).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# (narrow_mask*0+1).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# axes.set_extent([-100,-70, 17,40], crs=ccrs.PlateCarree())##x0, x1, y0, y1
# axes.text(.01,.8,'Tropical\nNorthwest Atlantic',transform=axes.transAxes,fontsize=24)
# (wide_mask*0+.5).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# (narrow_mask*0+1).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# axes.set_extent([63,100, 0,25], crs=ccrs.PlateCarree())##x0, x1, y0, y1
# axes.text(.01,.1,'Northern Indian Ocean',transform=axes.transAxes,fontsize=22)
# (wide_mask*0+.5).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# (narrow_mask*0+1).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
# axes.set_extent([-80,-50, -58,-30], crs=ccrs.PlateCarree())##x0, x1, y0, y1
# axes.text(.4,.93,'Patagonia',transform=axes.transAxes,fontsize=24)
(wide_mask*0+.5).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
(narrow_mask*0+1).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='Blues',vmin=0, vmax=1, add_colorbar=False)
axes.set_extent([117,150, 25,50], crs=ccrs.PlateCarree())##x0, x1, y0, y1
axes.text(.01,.8,'Subtropical \nNorthwest Pacific',transform=axes.transAxes,fontsize=24)
fig.tight_layout()
fig.savefig('nw_pac_wide_narrow_masks.pdf', format='pdf', bbox_inches='tight')
plt.rc('font', size = 16)
fig, axes = plt.subplots(3,1, figsize=(10,15), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
## pCO2
mnsiconc = seaice_model.load()
mn_pco2 = pco2_model #.where(mnsiconc<0.5) #- mean_pco2_atm
med_mn_pco2 = xr.DataArray(np.nanmedian(mn_pco2, axis=2), coords = [pco2.latitude, pco2.longitude], dims=['latitude', 'longitude'])
rms_cmodel=((mn_pco2-mn_pco2.mean('model'))**2).mean('model')**.5
mn_siconc_obs = seaice_obs
mn_pco2_obs = pco2_obs #.where(siconc_obs_mask<.5)
med_mn_pco2_obs = xr.DataArray(np.nanmedian(mn_pco2_obs, axis=0), coords = [pco2_model.latitude, pco2_model.longitude], dims=['latitude', 'longitude'])
rms_cobs=((mn_pco2_obs-mn_pco2_obs.mean('source'))**2).mean('source')**.5
mpl.rc('hatch', color='magenta', linewidth=5)
axes[0].coastlines()
gl=axes[0].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
med_mn_pco2_obs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'$\mu$atm'},levels=np.arange(290,420,10), cmap='YlGnBu')
rms_cobs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), colors='none', hatches=['','//'], levels=[-1,25,1e4], add_colorbar=False)
axes[0].set_title('Mean pCO2, 1998-2018\n4 pCO$_2$ Products Median')
axes[0].text(0,.9,'a', transform=axes[0].transAxes, fontsize=20, fontweight='bold')
axes[0].set_global()
axes[1].coastlines()
gl=axes[1].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(med_mn_pco2).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'$\mu$atm'},levels=np.arange(290,420,10), cmap='YlGnBu')
rms_cmodel.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), colors='none', hatches=['','\\\\'], levels=[-1,33,1e10], add_colorbar=False)
axes[1].set_title('15 Model Median')
axes[1].text(0,.9,'b', transform=axes[1].transAxes, fontsize=20, fontweight='bold')
axes[1].set_global()
mpl.rc('hatch', color='k', linewidth=5)
axes[2].coastlines()
gl=axes[2].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(med_mn_pco2-med_mn_pco2_obs).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'$\mu$atm'},levels=np.arange(-50,51,10), cmap='RdBu_r')
rms_cobs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), colors='none', hatches=['','//'], levels=[-1,25,1e4], add_colorbar=False)
axes[2].set_title('Models - Products')
axes[2].text(0,.9,'c', transform=axes[2].transAxes, fontsize=20, fontweight='bold')
axes[2].set_global()
fig.tight_layout()
plt.savefig('mean_pco2_map',format='pdf')
plt.rc('font', size = 18)
plt.rc('font', size = 16)
import matplotlib as mpl
# mpl.rc('hatch', color='k', linewidth=5)
fig, axes = plt.subplots(3,1, figsize=(10,15), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
# ----------------------------
# # FLUX
# mn_flx = fco2_model.mean('time')
# med_mn_flx = xr.DataArray(np.nanmedian(mn_flx, axis=0), coords = [fco2_model.latitude, fco2_model.longitude], dims=['latitude', 'longitude'])
# rms_fmodel=((mn_flx-mn_flx.mean('model'))**2).mean('model')**.5
# mn_flx_obs = fco2_obs.mean('time')
# med_mn_flx_obs = xr.DataArray(np.nanmedian(mn_flx_obs, axis=0), coords = [fco2_obs.latitude, fco2_obs.longitude], dims=['latitude', 'longitude'])
# rms_fobs=((mn_flx_obs-mn_flx_obs.mean('source'))**2).mean('source')**.5
axes[0].coastlines()
gl=axes[0].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
med_mn_flx_obs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'mol m$^{-2}$ yr$^{-1}$'}, levels=np.arange(-3,3.1,.5), cmap='PuOr_r')
rms_fobs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), colors='none', hatches=['','////'], levels=[-1,.6,100], add_colorbar=False)
axes[0].set_title('Mean Flux Density, 1998-2018\n4 pCO$_2$ Products Median')
axes[0].text(0,.9,'a', transform=axes[0].transAxes, fontsize=20, fontweight='bold')
axes[0].set_global()
axes[1].coastlines()
gl=axes[1].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(med_mn_flx).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'mol m$^{-2}$ yr$^{-1}$'},levels=np.arange(-3,3.1,.5), cmap='PuOr_r')
rms_fmodel.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), colors='none', hatches=['','\\\\\\'], levels=[-1,.95,100], add_colorbar=False)
axes[1].set_title('15 Model Median')
axes[1].text(0,.9,'b', transform=axes[1].transAxes, fontsize=20, fontweight='bold')
axes[1].set_global()
axes[2].coastlines()
gl=axes[2].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(med_mn_flx-med_mn_flx_obs).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'mol m$^{-2}$ yr$^{-1}$'},levels=np.arange(-3,3.1,.5), cmap='RdBu_r')
rms_fobs.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), colors='none', hatches=['','////'], levels=[-1,.6,100], add_colorbar=False)
rms_fmodel.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), colors='none', hatches=['','\\\\\\'], levels=[-1,.95,100], add_colorbar=False)
axes[2].set_title('Models - Products')
axes[2].text(0,.9,'c', transform=axes[2].transAxes, fontsize=20, fontweight='bold');
axes[2].set_global()
plt.tight_layout()
fig.savefig('mean_flux_maps.pdf',format='pdf',dpi=300,bbox_inches='tight',)
import matplotlib.patches as mpatches
fig, axes = plt.subplots(1,1, figsize=(10,10), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
axes.coastlines()
gl=axes.gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
fch4_weber.where(~np.isnan(wide_mask)).plot.contourf(ax=axes,levels=np.arange(0,.41,.02), cmap='YlGnBu',transform=ccrs.PlateCarree(),cbar_kwargs = {'shrink':.4, 'label':'g CH$_4$ m$^{-2}$ yr$^{-1}$',})
# cs.cmap.set_under('magenta')
# axes.text(.4,.15,'%s TgCH$_4$ yr$^{-1}$'%int_ch4.values.round(2),transform=axes.transAxes)
axes.add_patch(mpatches.Rectangle(xy=[-100,17], width=23, height=30,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[63,0], width=37, height=25,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[-80,-58], width=30, height=28,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.add_patch(mpatches.Rectangle(xy=[120,25], width=30, height=25,
facecolor='none', edgecolor='k',linewidth=2,
transform=ccrs.PlateCarree()) )
axes.set_title('CH$_4$ Flux, Weber CH$_4$ product',fontsize=20)
axes.set_global()
# fig.savefig('global_ch4.pdf', format='pdf', dpi=300, bbox_inches='tight')
# with xr.open_dataset('/tigress/GEOCLIM/LRGROUP/shared_data/bathymetry/etop5.nc') as ds:
# bathy0 = ds.bath.rename({'X':'longitude', 'Y':'latitude'})
# regridder = xe.Regridder(ds, wide_mask, "bilinear", periodic = True, ignore_degenerate = True)
# bathy = regridder(bathy0)
fig, axes = plt.subplots(1,1, figsize=(6,6), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=115)})
import cartopy.feature as cfeature
# ----------------------------
axes.coastlines()
gl=axes.gridlines()
axes.add_feature(cfeature.LAND,facecolor='lightgrey')
# fch4_weber.where(~np.isnan(wide_mask)).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='YlGnBu',levels=np.arange(0,.41,.01), add_colorbar=False)
# C = bathy.where(np.abs(bathy.latitude-29)<12).where(np.abs(bathy.longitude+85)<15).plot.contour(ax=axes, transform=ccrs.PlateCarree(),colors=['k'], levels=[-1000,-200, -50], linestyles=[':','--','-'], )
# axes.set_extent((-100,-70, 17,40), crs=ccrs.PlateCarree())
# axes.text(.01,.9,'Tropical Northwest Atlantic',transform=axes.transAxes,fontsize=20)
# fch4_weber.where(~np.isnan(wide_mask)).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='YlGnBu',levels=np.arange(0,.41,.01), add_colorbar=False)
# C = bathy.where(np.abs(bathy.latitude-12.5)<12.5).where(np.abs(bathy.longitude-81.5)<18.5).plot.contour(ax=axes, transform=ccrs.PlateCarree(),colors=['k'], levels=[-1000,-200, -50], linestyles=[':','--','-'], )
# axes.set_extent((63,100, 0,25), crs=ccrs.PlateCarree())
# axes.text(.1,.1,'Northern Indian Ocean',transform=axes.transAxes,fontsize=24, backgroundcolor='w', )
# fch4_weber.where(~np.isnan(wide_mask)).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='YlGnBu',levels=np.arange(0,.41,.01), add_colorbar=False)
# C = bathy.where(np.abs(bathy.latitude+45)<15).where(np.abs(bathy.longitude+65)<15).plot.contour(ax=axes, transform=ccrs.PlateCarree(),colors=['k'], levels=[-1000,-200, -50], linestyles=[':','--','-'], )
# axes.set_extent((-80,-50, -58,-30), crs=ccrs.PlateCarree())
# axes.text(.4,.93,'Patagonia',transform=axes.transAxes,fontsize=24)
fch4_weber.where(~np.isnan(wide_mask)).plot(ax=axes, transform=ccrs.PlateCarree(),cmap='YlGnBu',levels=np.arange(0,.41,.01), add_colorbar=False)
C = bathy.where(np.abs(bathy.latitude-37.5)<12.5).where(np.abs(bathy.longitude-135)<17).plot.contour(ax=axes, transform=ccrs.PlateCarree(),colors=['k'], levels=[-1000,-200, -50], linestyles=[':','--','-'], )
axes.set_extent((117,150, 25,50), crs=ccrs.PlateCarree())
axes.text(.01,.87,'Subtropical Northwest \nPacific',transform=axes.transAxes,fontsize=22)
labels=['1000','200','50']
for i in range(len(labels)):
C.collections[i].set_label(labels[i])
plt.legend(loc=(.01,.65))
fig.tight_layout()
fig.savefig('nw_pac_ch4.pdf', format='pdf', dpi=300, bbox_inches='tight')
model_n2o_mean = fn2o_model
fig, axes = plt.subplots(3,1, figsize=(10,15), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
axes[0].coastlines()
gl=axes[0].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 13, 'color': 'gray'}
gl.xlabel_style = {'size': 13, 'color': 'gray'}
fn2o_obs.plot.contourf(ax=axes[0], levels = np.arange(0,.11,.01),cmap='YlGnBu', transform=ccrs.PlateCarree(),cbar_kwargs = {'shrink':.6, 'label':'gN m$^{-2}$ yr$^{-1}$',})
# axes[0].text(.4,.15,'%s TgN yr$^{-1}$'%int_n2o_yang.values.round(2),transform=axes[1].transAxes)
axes[0].set_title('N$_2$O flux, Yang N$_2$O product',fontsize=20);
axes[0].set_global()
axes[0].text(0,1,'a',transform=axes[0].transAxes, fontweight='bold')
axes[1].coastlines()
gl=axes[1].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 13, 'color': 'gray'}
gl.xlabel_style = {'size': 13, 'color': 'gray'}
model_n2o_mean.median('model').where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1],levels = np.arange(0,.11,.01),cmap='YlGnBu', transform=ccrs.PlateCarree(),cbar_kwargs = {'shrink':.6, 'label':'gN m$^{-2}$ yr$^{-1}$',})
rms = ((model_n2o_mean-model_n2o_mean.mean('model'))**2).where(~np.isnan(wide_mask)).mean('model')**.5
rms.plot.contourf(ax=axes[1],transform=ccrs.PlateCarree(),levels=[0,0.016,100],colors='none',hatches=['','////'], add_colorbar=False)
axes[1].set_title('N$_2$O flux, 5-model median',fontsize=20)
axes[1].set_global()
axes[1].text(0,1,'b',transform=axes[1].transAxes, fontweight='bold')
axes[2].coastlines()
gl=axes[2].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 13, 'color': 'gray'}
gl.xlabel_style = {'size': 13, 'color': 'gray'}
(model_n2o_mean.mean('model') - fn2o_obs).plot.contourf(ax=axes[2], levels = np.arange(-.04,.041,.01),cmap='RdBu_r', transform=ccrs.PlateCarree(),cbar_kwargs = {'shrink':.6, 'label':'gN m$^{-2}$ yr$^{-1}$',})
rms.plot.contourf(ax=axes[2],transform=ccrs.PlateCarree(),levels=[0,0.016,100],colors='none',hatches=['','////'], add_colorbar=False)
axes[2].set_title('N$_2$O flux, Models - product',fontsize=20);
axes[2].set_global()
axes[2].text(0,1,'c',transform=axes[2].transAxes, fontweight='bold')
fig.tight_layout()
plt.savefig('n2o_models_obs_map.pdf',format='pdf',dpi=300,bbox_inches='tight')
# # # # START // LOAD DATA
with xr.open_dataset('mean_pco2_fgco2_latitude_models_simask.nc') as ds:
mean_fgco2_latitude = ds.fgco2_zonalmean[0:11,:] * 86400* 365#
mean_fgco2_latitude = mean_fgco2_latitude.where(mean_fgco2_latitude!=0)
mean_fgco2_latitude_clim = mean_fgco2_latitude.groupby('time.month').mean('time')
djf_fgco2_latitude = mean_fgco2_latitude_clim[:,[0,1,11],:].mean('month')
jja_fgco2_latitude = mean_fgco2_latitude_clim[:,5:8,:].mean('month')
diff_fgco2_latitude = djf_fgco2_latitude-jja_fgco2_latitude
diff_fgco2_latitude_median = np.nanmedian(diff_fgco2_latitude,0)
# diff_fgco2_latitude_25 = np.nanpercentile(djf_fgco2_latitude[0:11,:]-jja_fgco2_latitude[0:11,:],25,axis=0)
# diff_fgco2_latitude_75 = np.nanpercentile(djf_fgco2_latitude[0:11,:]-jja_fgco2_latitude[0:11,:],75,axis=0)
# # # # # #
with xr.open_dataset('mean_pco2_fgco2_latitude_obs_simask_1998-2018.nc') as ds:
mean_fgco2obs_latitude = ds.mean_fgco2_latitude
mean_fgco2obs_latitude_median = np.nanmedian(mean_fgco2obs_latitude.mean('time'),0)
mean_fgco2obs_latitude_clim = mean_fgco2obs_latitude.groupby('time.month').mean('time')
djf_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,[0,1,11],:].mean('month')
jja_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,5:8,:].mean('month')
diff_fgco2obs_latitude = djf_fgco2obs_latitude[0:11,:]-jja_fgco2obs_latitude[0:11,:]
diff_fgco2obs_latitude_median = np.nanmedian(diff_fgco2obs_latitude,0)
# # # # # END // LOAD DATA
# # # # #
fig,axs=plt.subplots(2,1, figsize=(10,10))
for k in range(10): mean_fgco2_latitude_clim.mean('month')[k,:].plot(ax = axs[0],label=diff_fgco2_latitude[k].model.values)
mean_fgco2_latitude_clim.mean('month')[10,:].plot(ax = axs[0],label=diff_fgco2_latitude[10].model.values, linestyle='--')
# axs[0].plot(mean_fgco2_latitude.latitude, mean_fgco2_latitude_median, color='g', linewidth=2, label='Models')
axs[0].plot(mean_fgco2_latitude.latitude, mean_fgco2obs_latitude_median, color='k', linewidth=2, label='pCO$_2$ Product Median\nand range')
axs[0].fill_between(mean_fgco2_latitude.latitude, mean_fgco2obs_latitude.mean('time').min('source'),mean_fgco2obs_latitude.mean('time').max('source'), color='grey',alpha=.2)
axs[0].set_title('a. Annual wide coastal flux density', fontweight='bold')
axs[0].set_xticks([])
diff_fgco2_latitude[0:10,:].plot(ax = axs[1], hue='model')
diff_fgco2_latitude[10,:].plot(ax = axs[1],linestyle='--')
axs[1].plot(mean_fgco2_latitude.latitude, diff_fgco2_latitude_median, linewidth=2, color='g')
axs[1].plot(mean_fgco2obs_latitude.latitude, diff_fgco2obs_latitude_median, linewidth=2, color='k')
axs[1].fill_between(mean_fgco2_latitude.latitude, diff_fgco2obs_latitude.min('source'),diff_fgco2obs_latitude.max('source'), color='grey',alpha=.2)
axs[1].set_title('b. Seasonal flux amplitude (DJF-JJA)', fontweight='bold'); axs[1].legend(loc=(0,0))
# # # #
axs[0].legend(loc=(1.05,-1))
# models_mean.add_legend(bbox_to_anchor=(0.5, -.05),ncol=2)
axs[0].set_ylim(-7,5); axs[1].set_ylim(-7,5);
axs[0].set_xlim(-80,90); axs[1].set_xlim(-80,90);
axs[0].set_ylabel('[mol m$^{-2}$ yr$^{-1}$]'); axs[1].set_ylabel('[mol m$^{-2}$ yr$^{-1}$] ');
axs[0].set_xlabel('')
axs[0].plot(np.arange(-80,90),np.zeros(170),color='grey'); axs[1].plot(np.arange(-80,90),np.zeros(170),color='grey',);
axs[0].axvline(0,color='grey',linewidth=.8,linestyle=':'); axs[1].axvline(0,color='grey',linewidth=.8,linestyle=':');
axs[1].text(5,-6,'Stronger sink \nin fall/winter',fontsize=12)
axs[1].text(5,3,'Stronger sink \nin spring/summer',fontsize=12)
axs[1].text(-40,-6,'Stronger sink \nin spring/summer',fontsize=12)
axs[1].text(-75,3,'Stronger sink \nin fall/winter',fontsize=12)
fig.savefig('mean_flux_latitude_modelcolors.pdf', format='pdf', bbox_inches='tight')
# # # # START // LOAD DATA
with xr.open_dataset('mean_pco2_fgco2_latitude_models_simask.nc') as ds:
mean_fgco2_latitude = ds.fgco2_zonalmean[0:11,:]*86400*365
mean_fgco2_latitude = mean_fgco2_latitude.where(mean_fgco2_latitude!=0)
mean_fgco2_latitude_clim = mean_fgco2_latitude.groupby('time.month').mean('time')
djf_fgco2_latitude = mean_fgco2_latitude_clim[:,[0,1,11],:].mean('month')
jja_fgco2_latitude = mean_fgco2_latitude_clim[:,5:8,:].mean('month')
son_fgco2_latitude = mean_fgco2_latitude_clim[:,8:11,:].mean('month')
mam_fgco2_latitude = mean_fgco2_latitude_clim[:,2:5,:].mean('month')
# # # # # #
with xr.open_dataset('mean_pco2_fgco2_latitude_obs_simask_1998-2018.nc') as ds:
mean_fgco2obs_latitude = ds.mean_fgco2_latitude
mean_fgco2obs_latitude_clim = mean_fgco2obs_latitude.groupby('time.month').mean('time')
djf_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,[0,1,11],:].mean('month')
jja_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,5:8,:].mean('month')
son_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,8:11,:].mean('month')
mam_fgco2obs_latitude = mean_fgco2obs_latitude_clim[:,2:5,:].mean('month')
# # # # END // LOAD DATA
# # # #
fig,axs=plt.subplots(3,2, figsize=(20,16))
mean_fgco2_latitude_clim.mean('month').plot(hue='model', ax = axs[0][0], color='yellowgreen', linewidth=.8, )
axs[0][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(mean_fgco2_latitude_clim.mean('month'),0), color='g', linewidth=2, label='Models')
mean_fgco2obs_latitude.mean('time').plot(hue='source', ax = axs[0][0], color='steelblue', linewidth=.8, )
axs[0][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(mean_fgco2obs_latitude.mean('time'),0), color='steelblue', linewidth=2, label='pCO$_2$ Products')
axs[0][0].set_title('Annual mean flux density')
axs[0][0].legend(loc=(.1,.01),fontsize=20)
axs[0][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[0][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[0][0].set_xlim(-80,90)
(djf_fgco2_latitude-jja_fgco2_latitude).plot(hue='model', ax = axs[0][1], color='yellowgreen', linewidth=.8, )
axs[0][1].plot(mean_fgco2_latitude.latitude, np.nanmedian((djf_fgco2_latitude-jja_fgco2_latitude),0), color='g', linewidth=2,)
(djf_fgco2obs_latitude-jja_fgco2obs_latitude).plot(hue='source', ax = axs[0][1], color='steelblue', linewidth=.8, )
axs[0][1].plot(mean_fgco2_latitude.latitude, np.nanmedian((djf_fgco2obs_latitude-jja_fgco2obs_latitude),0), color='steelblue', linewidth=2,)
axs[0][1].set_title('Seasonal amplitude (DJF - JJA)')
axs[0][1].legend(loc=(1,.9))
axs[0][1].plot(np.arange(-80,80),np.zeros(160),'k')
axs[0][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[0][1].set_xlim(-80,90)
mam_fgco2_latitude.plot(hue='model', ax = axs[1][0], color='yellowgreen', linewidth=.8, )
axs[1][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(mam_fgco2_latitude,0), color='g', linewidth=2, )
mam_fgco2obs_latitude.plot(hue='source', ax = axs[1][0], color='steelblue', linewidth=.8, )
axs[1][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(mam_fgco2obs_latitude,0), color='steelblue', linewidth=2, )
axs[1][0].set_title('MAM flux density')
axs[1][0].legend(loc=(1,.9))
axs[1][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[1][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[1][0].text(65,4,'spring');axs[1][0].text(-75,4,'fall')
axs[1][0].set_xlim(-80,90)
jja_fgco2_latitude.plot(hue='model', ax = axs[1][1], color='yellowgreen', linewidth=.8, )
axs[1][1].plot(mean_fgco2_latitude.latitude, np.nanmedian(jja_fgco2_latitude,0), color='g', linewidth=2, )
jja_fgco2obs_latitude.plot(hue='source', ax = axs[1][1], color='steelblue', linewidth=.8, )
axs[1][1].plot(mean_fgco2_latitude.latitude, np.nanmedian(jja_fgco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[1][1].set_title('JJA flux density')
axs[1][1].legend(loc=(1,.9))
axs[1][1].plot(np.arange(-100,100),np.zeros(200),'k')
axs[1][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[1][1].text(65,4,'summer');axs[1][1].text(-75,4,'winter')
axs[1][1].set_xlim(-80,90)
son_fgco2_latitude.plot(hue='model', ax = axs[2][0], color='yellowgreen', linewidth=.8, )
axs[2][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(son_fgco2_latitude,0), color='g', linewidth=2, )
son_fgco2obs_latitude.plot(hue='source', ax = axs[2][0], color='steelblue', linewidth=.8, )
axs[2][0].plot(mean_fgco2_latitude.latitude, np.nanmedian(son_fgco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[2][0].set_title('SON flux density')
axs[2][0].legend(loc=(1,.9))
axs[2][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[2][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[2][0].text(65,4,'fall');axs[2][0].text(-75,4,'spring')
axs[2][0].set_xlim(-80,90)
djf_fgco2_latitude.plot(hue='model', ax = axs[2][1], color='yellowgreen', linewidth=.8, )
axs[2][1].plot(mean_fgco2_latitude.latitude, np.nanmedian(djf_fgco2_latitude,0), color='g', linewidth=2,)
djf_fgco2obs_latitude.plot(hue='source', ax = axs[2][1], color='steelblue', linewidth=.8, )
axs[2][1].plot(mean_fgco2_latitude.latitude, np.nanmedian(djf_fgco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[2][1].set_title('DJF flux density')
axs[2][1].legend()
axs[2][1].plot(np.arange(-100,100),np.zeros(200),'k')
axs[2][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[2][1].text(65,4,'winter'); axs[2][1].text(-75,4,'summer')
axs[2][1].set_xlim(-80,90)
axs[0][0].set_ylim(-7,5)
axs[0][1].set_ylim(-7,5)
axs[1][0].set_ylim(-7,5)
axs[1][1].set_ylim(-7,5)
axs[2][0].set_ylim(-7,5)
axs[2][1].set_ylim(-7,5)
axs[0][0].axvline(0, color='k',linewidth=.5)
axs[0][1].axvline(0, color='k',linewidth=.5)
axs[1][0].axvline(0, color='k',linewidth=.5)
axs[1][1].axvline(0, color='k',linewidth=.5)
axs[2][0].axvline(0, color='k',linewidth=.5)
axs[2][1].axvline(0, color='k',linewidth=.5)
axs[0][1].text(45,3,'Stronger sink in \nsummer')
axs[0][1].text(45,-6.5,'Stronger sink in \nwinter')
axs[0][1].text(-78,3,'Stronger sink in \nwinter')
axs[0][1].text(-78,-6.5,'Stronger sink in \nsummer')
axs[0][0].text(-80,5.2,'a', fontweight='bold')
axs[0][1].text(-80,5.2,'b', fontweight='bold')
axs[1][0].text(-80,5.2,'c', fontweight='bold')
axs[1][1].text(-80,5.2,'d', fontweight='bold')
axs[2][0].text(-80,5.2,'e', fontweight='bold')
axs[2][1].text(-80,5.2,'f', fontweight='bold')
fig.tight_layout()
fig.savefig('flux_latitude_seasonal.pdf', format='pdf', dpi=300, bbox_inches='tight')
/home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, No handles with labels found to put in legend. /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, No handles with labels found to put in legend. /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, No handles with labels found to put in legend. /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, No handles with labels found to put in legend. /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, /home/hogikyan/anaconda3/envs/myenv_akh/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out, No handles with labels found to put in legend.
# # # # START // LOAD DATA
with xr.open_dataset('mean_pco2_fgco2_latitude_models_simask.nc') as ds:
mean_pco2_latitude = ds.pco2_latitude[0:11,:] ##other models are regional.
mean_pco2_latitude_clim = mean_pco2_latitude#.groupby('time.month').mean('time')
djf_pco2_latitude = mean_pco2_latitude_clim[:,[0,1,11],:].mean('month')
jja_pco2_latitude = mean_pco2_latitude_clim[:,5:8,:].mean('month')
son_pco2_latitude = mean_pco2_latitude_clim[:,8:11,:].mean('month')
mam_pco2_latitude = mean_pco2_latitude_clim[:,2:5,:].mean('month')
mean_pco2_atm = xr.DataArray(384.85)
# # # # # #
with xr.open_dataset('mean_pco2_fgco2_latitude_obs_simask_1998-2018.nc') as ds:
mean_pco2obs_latitude = ds.mean_pco2_latitude
mean_pco2obs_latitude_clim = mean_pco2obs_latitude.groupby('time.month').mean('time')
djf_pco2obs_latitude = mean_pco2obs_latitude_clim[:,[0,1,11],:].mean('month')
jja_pco2obs_latitude = mean_pco2obs_latitude_clim[:,5:8,:].mean('month')
son_pco2obs_latitude = mean_pco2obs_latitude_clim[:,8:11,:].mean('month')
mam_pco2obs_latitude = mean_pco2obs_latitude_clim[:,2:5,:].mean('month')
# # # # END // LOAD DATA
# # # #
fig,axs=plt.subplots(3,2, figsize=(20,15))
mean_pco2_latitude.mean('month').plot(hue='model', ax = axs[0][0], color='yellowgreen', linewidth=.5, )
axs[0][0].plot(mean_pco2_latitude.latitude, np.nanmedian(mean_pco2_latitude.mean('month'),0), color='g', linewidth=2, label='Models')
mean_pco2obs_latitude.mean('time').plot(hue='source', ax = axs[0][0], color='steelblue', linewidth=.5, )
axs[0][0].plot(mean_pco2_latitude.latitude, np.nanmedian(mean_pco2obs_latitude.mean('time'),0), color='steelblue', linewidth=2, label='pCO$_2$ Products')
axs[0][0].set_title('Annual mean pCO$_2$')
axs[0][0].legend(loc=(.1,.01),fontsize=20)
axs[0][0].plot(np.arange(-100,100),np.zeros(200)+mean_pco2_atm.values,'k')
axs[0][0].set_ylabel('$\mu$atm')
axs[0][0].set_xlim(-80,90);axs[0][0].set_ylim(250,500)
(djf_pco2_latitude-jja_pco2_latitude).plot(hue='model', ax = axs[0][1], color='yellowgreen', linewidth=.5, )
axs[0][1].plot(mean_pco2_latitude.latitude, np.nanmedian((djf_pco2_latitude-jja_pco2_latitude),0), color='g', linewidth=2,)
(djf_pco2obs_latitude-jja_pco2obs_latitude).plot(hue='source', ax = axs[0][1], color='steelblue', linewidth=.5, )
axs[0][1].plot(mean_pco2_latitude.latitude, np.nanmedian((djf_pco2obs_latitude-jja_pco2obs_latitude),0), color='steelblue', linewidth=2,)
axs[0][1].set_title('Seasonal amplitude (DJF - JJA)')
axs[0][1].legend(loc=(1,.9))
axs[0][1].plot(np.arange(-100,100),np.zeros(200),'k')
axs[0][1].set_ylabel('$\mu$atm')
axs[0][1].set_xlim(-80,90);axs[0][1].set_ylim(-250,250)
mam_pco2_latitude.plot(hue='model', ax = axs[1][0], color='yellowgreen', linewidth=.5, )
axs[1][0].plot(mean_pco2_latitude.latitude, np.nanmedian(mam_pco2_latitude,0), color='g', linewidth=2, )
mam_pco2obs_latitude.plot(hue='source', ax = axs[1][0], color='steelblue', linewidth=.5, )
axs[1][0].plot(mean_pco2_latitude.latitude, np.nanmedian(mam_pco2obs_latitude,0), color='steelblue', linewidth=2, )
axs[1][0].set_title('MAM pCO$_2$')
axs[1][0].legend(loc=(1,.9))
axs[1][0].plot(np.arange(-100,100),np.zeros(200)+mean_pco2_atm.values,'k')
axs[1][0].set_ylabel('$\mu$atm')
axs[1][0].set_xlim(-80,90);axs[1][0].set_ylim(250,500)
axs[1][0].text(65,480,'spring')
axs[1][0].text(-78,480,'fall')
jja_pco2_latitude.plot(hue='model', ax = axs[1][1], color='yellowgreen', linewidth=.5, )
axs[1][1].plot(mean_pco2_latitude.latitude, np.nanmedian(jja_pco2_latitude,0), color='g', linewidth=2, )
jja_pco2obs_latitude.plot(hue='source', ax = axs[1][1], color='steelblue', linewidth=.5, )
axs[1][1].plot(mean_pco2_latitude.latitude, np.nanmedian(jja_pco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[1][1].set_title('JJA pCO$_2$')
axs[1][1].legend(loc=(1,.9))
axs[1][1].plot(np.arange(-100,100),np.zeros(200)+mean_pco2_atm.values,'k')
axs[1][1].set_ylabel('$\mu$atm')
axs[1][1].set_xlim(-80,90);axs[1][1].set_ylim(250,500)
axs[1][1].text(65,480,'summer')
axs[1][1].text(-78,480,'winter')
son_pco2_latitude.plot(hue='model', ax = axs[2][0], color='yellowgreen', linewidth=.5, )
axs[2][0].plot(mean_pco2_latitude.latitude, np.nanmedian(son_pco2_latitude,0), color='g', linewidth=2, )
son_pco2obs_latitude.plot(hue='source', ax = axs[2][0], color='steelblue', linewidth=.5, )
axs[2][0].plot(mean_pco2_latitude.latitude, np.nanmedian(son_pco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[2][0].set_title('SON pCO$_2$')
axs[2][0].legend(loc=(1,.9))
axs[2][0].plot(np.arange(-100,100),np.zeros(200)+mean_pco2_atm.values,'k')
axs[2][0].set_ylabel('$\mu$atm')
axs[2][0].set_xlim(-80,90);axs[2][0].set_ylim(250,500)
axs[2][0].text(65,480,'fall')
axs[2][0].text(-78,480,'spring')
djf_pco2_latitude.plot(hue='model', ax = axs[2][1], color='yellowgreen', linewidth=.5, )
axs[2][1].plot(mean_pco2_latitude.latitude, np.nanmedian(djf_pco2_latitude,0), color='g', linewidth=2,)
djf_pco2obs_latitude.plot(hue='source', ax = axs[2][1], color='steelblue', linewidth=.5, )
axs[2][1].plot(mean_pco2_latitude.latitude, np.nanmedian(djf_pco2obs_latitude,0), color='steelblue', linewidth=2,)
axs[2][1].set_title('DJF pCO$_2$')
axs[2][1].legend(loc=(1,.9))
axs[2][1].plot(np.arange(-100,100),np.zeros(200)+mean_pco2_atm.values,'k')
axs[2][1].set_ylabel('$\mu$atm')
axs[2][1].set_xlim(-80,90);axs[2][1].set_ylim(250,500)
axs[2][1].text(65,480,'winter')
axs[2][1].text(-78,480,'summer')
# axs[0][0].set_ylim(-7,5)
# axs[0][1].set_ylim(-7,5)
# axs[1][0].set_ylim(-7,5)
# axs[1][1].set_ylim(-7,5)
# axs[2][0].set_ylim(-7,5)
# axs[2][1].set_ylim(-7,5)
axs[0][0].axvline(0, color='k',linewidth=.5)
axs[0][1].axvline(0, color='k',linewidth=.5)
axs[1][0].axvline(0, color='k',linewidth=.5)
axs[1][1].axvline(0, color='k',linewidth=.5)
axs[2][0].axvline(0, color='k',linewidth=.5)
axs[2][1].axvline(0, color='k',linewidth=.5)
axs[0][1].text(45,160,'Stronger sink in \nsummer')
axs[0][1].text(45,-240,'Stronger sink in \nwinter')
axs[0][1].text(-78,160,'Stronger sink in \nwinter')
axs[0][1].text(-78,-240,'Stronger sink in \nsummer')
axs[0][0].text(-80,510,'a', fontweight='bold')
axs[0][1].text(-80,260,'b', fontweight='bold')##amp
axs[1][0].text(-80,510,'c', fontweight='bold')
axs[1][1].text(-80,510,'d', fontweight='bold')
axs[2][0].text(-80,510,'e', fontweight='bold')
axs[2][1].text(-80,510,'f', fontweight='bold')
fig.tight_layout()
fig.savefig('mean_pco2_latitude.pdf', format='pdf', dpi=300, bbox_inches='tight')
### LOADED AT BEGINNING.
blues = ['skyblue','steelblue','darkblue','darkblue','mediumturquoise']
fig, axs = plt.subplots(1,2, figsize=(17,7))
ax1 = axs[1]
ax2 = axs[0]
a=[np.zeros(8)+1e16,mean_fco2_wideshelf.mean('time'), mean_fco2_open.mean('time')]
vp=ax1.violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
b=[mean_fco2_thinshelf.mean('time')[[1,4,6,7]],mean_fco2_wideshelf.mean('time')[[1,4,6,7]], mean_fco2_open.mean('time')[[1,4,6,7]]]
vp=ax1.violinplot(b, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
ax1.scatter(1,mean_fco2_thinshelf.mean('time')[7],color='goldenrod')
ax1.scatter(2,mean_fco2_wideshelf.mean('time')[7],color='goldenrod')
ax1.scatter(3,mean_fco2_open.mean('time')[7],color='goldenrod')
mn_f_glob_thinshelf_obs_timemean = mean_fco2obs_thinshelf.mean('time')
mn_f_glob_thinshelf_obs_timemean_median = np.nanmedian(mn_f_glob_thinshelf_obs_timemean)
ax1.scatter(1,mn_f_glob_thinshelf_obs_timemean[0], color=blues[0], marker='>',s=200)
ax1.scatter(1,mn_f_glob_thinshelf_obs_timemean[1], color=blues[1], marker='^',s=200)
ax1.scatter(1,mn_f_glob_thinshelf_obs_timemean[2], color=blues[2], s=200)
ax1.scatter(1,mn_f_glob_thinshelf_obs_timemean[3], color=blues[3], marker='x',s=300)
ax1.scatter(1,mean_fco2obs_wideshelf_modelkw.mean('time'),edgecolor=blues[3],color='none',linewidth=3,s=200)
ax1P = ax1.twinx()
ax1P.plot(np.arange(.8,1.21,.05), (np.zeros(9)+1)*mn_f_glob_thinshelf_obs_timemean_median, color=blues[3], linewidth=2)
ax1P.errorbar(.7,-.68,yerr=.14,marker='o', color='grey',capsize=10, label='Dai et al. 2022') # # Dai et al
ax1P.errorbar(.65,-1.09,yerr=2.9, marker='d', color='grey', markersize=20, capsize=10,label='Chen et al. 2013')# # Chen et al continental shelves
ax1P.scatter(.6,-.7, marker='>', color='grey', s=200, label='Laruelle et al. 2014')
mn_f_glob_obs_mean = mean_fco2obs_wideshelf.mean('time')
mn_f_glob_obs_mean_median = np.nanmedian(mn_f_glob_obs_mean)
ax1.scatter(2,mn_f_glob_obs_mean[0], color=blues[0], marker='>',s=200)
ax1.scatter(2,mn_f_glob_obs_mean[1], color=blues[1], marker='^',s=200)
ax1.scatter(2,mn_f_glob_obs_mean[2], color=blues[2], s=200)
ax1.scatter(2,mn_f_glob_obs_mean[3], color=blues[3], marker='x',s=300)
ax1.plot(np.arange(1.8,2.21,.05), (np.zeros(9)+1)*mn_f_glob_obs_mean_median, color=blues[3], linewidth=2)
ax1.scatter(2, mean_fco2obs_wideshelf_modelkw.mean('time'), edgecolor = blues[3], color='none',linewidth=3, s=200) #label = 'Coastal-SOM-FFN with MOM6 kw',
# # mn_f_coastalmom6kw
mean_fco2obs_open_mean = mean_fco2obs_open.mean('time')
mean_fco2obs_open_mean_median = np.nanmedian(mean_fco2obs_open_mean)
ax1.scatter(3,mean_fco2obs_open_mean[0], color=blues[0], marker='>',s=200)
ax1.scatter(3,mean_fco2obs_open_mean[1], color=blues[1], marker='^',s=200)
ax1.scatter(3,mean_fco2obs_open_mean[3], color=blues[3], marker='x',s=200)
ax1.plot(np.arange(2.8,3.21,.05), (np.zeros(9)+1)*mean_fco2obs_open_mean_median, color=blues[3], linewidth=2)
ax1.plot(np.arange(-5,-3), np.arange(-5,-3), color='yellowgreen', linewidth=5, label='11 global models')
ax1.plot(np.arange(-5,-3), np.arange(-5,-3), color='goldenrod', linewidth=5, label='4 high-res models')
ax1.set_title('Mean CO$_2$ Flux Density')
ax1.set_xticks([1,2,3]); ax1.set_xticklabels(['Narrow Coastal','Wide Coastal', 'Open Ocean'], rotation=90);
ax1.set_ylabel('mol m$^{-2}$ yr$^{-1}$')
ax1.plot(np.arange(-1,5),np.zeros(6),'k')
ax1.set_ylim(-1.25,0.12); ax1.set_xlim(0.5,3.5)
ax1.legend(loc=(1.05,.8), fontsize=16,)#edgecolor='none')
ax1P.legend(loc=(1.05,.04), fontsize=16,)#edgecolor='none')
ax1P.set_ylim(-1.25,0.12);ax1P.set_yticks([])
# --------------------------------------------
a=[np.zeros(8), mean_pco2_wideshelf.mean('time'), mean_pco2_open.mean('time')]
vp=ax2.violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
b=[mean_pco2_thinshelf[:,[1,4,6,7]].mean('time'), mean_pco2_wideshelf[:,[1,4,6,7]].mean('time'), mean_pco2_open[:,[1,4,6,7]].mean('time')]
vp=ax2.violinplot(b, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
ax2.scatter(1,mean_pco2_thinshelf.mean('time')[7],color='goldenrod')
ax2.scatter(2,mean_pco2_wideshelf.mean('time')[7],color='goldenrod')
ax2.scatter(3,mean_pco2_open.mean('time')[7],color='goldenrod')
mn_pco2_thinshelf_obs_mean = mean_pco2obs_thinshelf.mean('time')
mn_pco2_thinshelf_obs_mean_median = np.nanmedian(mn_pco2_thinshelf_obs_mean)
ax2.scatter(1,mn_pco2_thinshelf_obs_mean[0], color=blues[0], marker='>',s=200)
ax2.scatter(1,mn_pco2_thinshelf_obs_mean[1], color=blues[1], marker='^',s=200)
ax2.scatter(1,mn_pco2_thinshelf_obs_mean[2], color=blues[2], s=200)
ax2.scatter(1,mn_pco2_thinshelf_obs_mean[3], color=blues[3], marker='x',s=300)
ax2.plot(np.arange(.8,1.21,.05), (np.zeros(9)+1)*mn_pco2_thinshelf_obs_mean_median, color=blues[3], linewidth=2)
mn_pco2_glob_obs_mean = mean_pco2obs_wideshelf.mean('time')
mn_pco2_glob_obs_mean_median = np.nanmedian(mn_pco2_glob_obs_mean)
ax2.scatter(2,mn_pco2_glob_obs_mean[0], color=blues[0], marker='>',s=200, label='Carboscope-1')
ax2.scatter(2,mn_pco2_glob_obs_mean[1], color=blues[1], marker='^',s=200, label='CMEMS*')
ax2.scatter(2,mn_pco2_glob_obs_mean[2], color=blues[2], s=200, label='Coastal SOM-FFN')
ax2.scatter(2,mn_pco2_glob_obs_mean[3], color=blues[3], marker='x',s=300, label='Merged SOM-FFN')
ax2.plot(np.arange(1.8,2.21,.05), (np.zeros(9)+1)*mn_pco2_glob_obs_mean_median, color=blues[3], linewidth=2, label='Product Median')
ax2.scatter(2, -500, label = 'Coastal-SOM-FFN-kw', edgecolor = blues[3],color='none',linewidth=3, s=200)
# axs[2].scatter(1.9,pco2_glob_obs[13*12:-12].mean('time'), color='k',s=180)
mean_oo_pco2_obs_mean = mean_pco2obs_open.mean('time')
mean_oo_pco2_obs_mean_median = np.nanmedian(mean_oo_pco2_obs_mean)
ax2.scatter(3,mean_oo_pco2_obs_mean[0], color=blues[0], marker='>',s=200)
ax2.scatter(3,mean_oo_pco2_obs_mean[1], color=blues[1], marker='^',s=200)
ax2.scatter(3,mean_oo_pco2_obs_mean[3], color=blues[3], marker='x',s=300)
ax2.plot(np.arange(2.8,3.21,.05), (np.zeros(9)+1)*mean_oo_pco2_obs_mean_median, color=blues[3], linewidth=2)
ax2.plot(np.arange(0,5),np.zeros(5)+pco2_air.mean('time').values,color='k')
ax2.text(.6,pco2_air.mean('time').values+1,'Atmospheric pCO$_2$')
ax2.set_title('Mean surface pCO$_2$') # #and bias is not improved if we end in 2017 with socat
ax2.set_xticks([1,2,3]); ax2.set_xticklabels(['Narrow Coastal','Wide Coastal', 'Open Ocean'], rotation=90);
ax2.set_xlim(.5,3.5)
ax2.set_ylabel('$\mu$atm')
ax2.legend(loc=(2.23,.32), fontsize=16,)#edgecolor='none')
ax2.text(2.22,.25,'Prior Studies',transform=ax2.transAxes,fontsize=18)
ax2.text(2.22,.73,'pCO$_2$ Products',transform=ax2.transAxes,fontsize=18)
ax2.text(2.22,.94,'Models',transform=ax2.transAxes,fontsize=18)
ax2.set_ylim(340,388);
# fig.tight_layout()
rect = patches.Rectangle((2.22,.01),.65,.98,transform=ax2.transAxes,edgecolor='grey', facecolor='none',clip_on=False)
ax2.add_patch(rect)
ax2.text(-.1,1.01,'a',fontweight='bold',transform=ax2.transAxes)
ax1.text(-.1,1.01,'b',fontweight='bold',transform=ax1.transAxes)
# fig.savefig('./mean_pco2_flux_narrow_wide_open.pdf', bbox_inches='tight', format='pdf')
int_fco2_thinshelf = mean_fco2_thinshelf * total_area_narrow
int_fco2_wideshelf = mean_fco2_wideshelf * total_area_wide
int_fco2obs_thinshelf = mean_fco2obs_thinshelf * total_area_narrow
int_fco2obs_thinshelf_modelkw = mean_fco2obs_thinshelf_modelkw * total_area_narrow
int_fco2obs_wideshelf = mean_fco2obs_wideshelf * total_area_wide
int_fco2obs_wideshelf_modelkw = mean_fco2obs_wideshelf_modelkw * total_area_wide
blues = ['skyblue','steelblue','darkblue','darkblue','mediumturquoise']
fig = plt.figure()
ax0C = fig.add_axes((0,1.6,.8,1.4), );ax1C = ax0C.twinx();ax2C = ax0C.twinx()
a=[np.zeros(8)+1e6,int_fco2_wideshelf.mean('time')*12.011/1e15,np.zeros(8)+1e6]
vp=ax0C.violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
b=[int_fco2_thinshelf.mean('time')[[1,4,6,7]]*12.011/1e15,int_fco2_wideshelf.mean('time')[[1,4,6,7]]*12.011/1e15,np.zeros(8)]
vp=ax0C.violinplot(b, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
ax0C.plot(np.arange(2),np.zeros(2)-50,label='11 global models', color='yellowgreen',linewidth=4)
ax0C.plot(np.arange(2),np.zeros(2)-50,label='4 high res. models', color='goldenrod',linewidth=4)
f_glob_thinshelf_obs_int = int_fco2obs_thinshelf.mean('time')*12.011/1e15
f_glob_thinshelf_obs_int_median = np.nanmedian(f_glob_thinshelf_obs_int)
ax1C.scatter(1,f_glob_thinshelf_obs_int[0], color=blues[0], marker = '>',s=200, label='Carboscope-1')
ax1C.scatter(1,f_glob_thinshelf_obs_int[1], color=blues[1], marker='^',s=200, label='CMEMS*')
ax1C.scatter(1,f_glob_thinshelf_obs_int[2], color=blues[2], s=200, label='Coastal-SOM-FFN')
ax1C.scatter(1,f_glob_thinshelf_obs_int[3], color=blues[3], marker='x',s=300, label='Merged SOM-FFN')
ax1C.scatter(1,int_fco2obs_thinshelf_modelkw.mean('time')*12.011/1e15, edgecolor=blues[3], color='none',linewidth=3,s=200, label='Coastal-SOM-FFN-kw')
ax1C.plot(np.arange(.8,1.21,.05),(np.zeros(9)+1)*f_glob_thinshelf_obs_int_median, color=blues[3], linewidth=2, label='Median')
ax2C.scatter(.75,-.4, color='grey', marker='d', s=100, label='Chen et al. 2013')
ax2C.errorbar(.75,-.25,yerr=.05,marker='o', color='grey',capsize=3, label='Dai et al. 2022')
ax2C.scatter(.7,-.26,color='grey', marker='<',s=100, label='Laruelle et al. 2018')
ax2C.errorbar(.65,-.32,yerr=.08,color='grey',marker='*', capsize=3,markersize=10, label='Regnier et al. 2022')
ax2C.scatter(.65,-.2,color='grey', marker='^',s=100, label='Regnier et al. 2013')
ax2C.errorbar(.6,-.19,yerr=.5,color='grey',marker='>', capsize=3,markersize=10, label='Laruelle et al. 2014')
ax2C.scatter(.57,-.25,color='grey', marker='x',s=100, label='Bauer et al. 2013')
ax2C.scatter(.75,-.1,color='grey', marker='p',s=100, label='Lacroix et al. 2021')
ax2C.scatter(.77,-.15,color='grey', marker='P',s=100, label='Bourgeois et al. 2016')
f_glob_obs_int = int_fco2obs_wideshelf.mean('time')*12.011/1e15
f_glob_obs_int_median = np.nanmedian(f_glob_obs_int)
ax1C.scatter(2,f_glob_obs_int[0], color=blues[0], marker='>',s=200)
ax1C.scatter(2,f_glob_obs_int[1], color=blues[1], marker='^',s=200)
ax1C.scatter(2,f_glob_obs_int[2], color=blues[2], s=200)
ax1C.scatter(2,f_glob_obs_int[3], color=blues[3], marker='x',s=300)
ax1C.plot(np.arange(1.8,2.21,.05), (np.zeros(9)+1)*f_glob_obs_int_median, color=blues[3], linewidth=2)
ax1C.scatter(2, int_fco2obs_wideshelf_modelkw.mean('time')*12.011/1e15, edgecolor = blues[3],linewidth=3, color='none',s=200,)# label = 'Coastal-SOM-FFN-kw')
ax0C.set_title('Total CO$_2$ Flux')
ax0C.set_xticks([1,2]); ax0C.set_xticklabels(['Narrow coast', 'Wide coast', ],);
ax0C.set_ylabel('PgC yr$^{-1}$')
ax0C.legend(loc=(1.03,.87), fontsize=15);ax0C.text(1.03,1,'Model distribution and Median', transform = ax0C.transAxes)
ax1C.legend(loc=(1.03,.41), fontsize=15);ax1C.text(1.03,.78,'pCO$_2$ Products', transform = ax1C.transAxes)
ax2C.legend(loc=(1.68,.41), fontsize=15);ax2C.text(1.68,.94,'Prior Estimates', transform = ax2C.transAxes)
ax0C.plot(np.arange(-1,5),np.zeros(6),'k')
ax0C.set_ylim(-1,.1); ax0C.set_xlim(.5,2.5)
ax1C.set_ylim(-1,.1); #ax1C.set_xlim(.5,2.5)
ax2C.set_ylim(-1,.1); #ax2C.set_xlim(.5,2.5)
rect = patches.Rectangle((1.01,.4),1.3,.65,transform=ax0C.transAxes,edgecolor='grey', facecolor='white',clip_on=False)
ax0C.add_patch(rect)
ax1C.set_yticks([]);ax2C.set_yticks([])
# ---------------------------------
blues = ['skyblue','steelblue','darkblue','darkblue','mediumturquoise']
ax0N = fig.add_axes((0,0,.8,1.4),);ax1N=ax0N.twinx()
ax0N.plot(np.arange(-1,1),np.zeros(2)-1,'yellowgreen',linewidth=4,label='5 global models')
a=[[ 0.35, 0.17, 0.14],[ 0.64,0.60,1.39,],]
vp=ax0N.violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
ax0N.plot(np.arange(-1,1),np.zeros(2)-1,'goldenrod',linewidth=4,label='3 high res. models')
a=[[1e6],[0.64,0.60,1.73,1.29,1.27],]
vp=ax0N.violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
ax1N.scatter(1,0.70, color=blues[0], marker = '>',s=200) ## Yang
ax1N.scatter(1,0.90, color=blues[1], marker='^',s=200, )#label='MARCATS-N2O')
ax1N.scatter(2,1.63, color=blues[0], marker='>',s=200, label='Yang-N2O')
ax1N.scatter(2,3.55, color=blues[1], marker='^',s=200, label='MARCATS-N2O')
ax0N.set_title('Total N$_2$O Flux')
ax0N.set_xticks([1,2]); ax0N.set_xticklabels(['Narrow coast', 'Wide coast', ], );
ax0N.set_ylabel('TgN yr$^{-1}$')
ax0N.legend(loc=(1.02,.82), fontsize=15)
ax0N.text(1.03,.95,'Model distribution \nand Median', transform=ax0N.transAxes)
ax1N.legend(loc=(1.02,.61), fontsize=15)
ax1N.text(1.03,.75,'N$_2$O Products', transform=ax1N.transAxes)
ax0N.set_ylim(0,4); ax0N.set_xlim(.5,2.5);ax1N.set_ylim(0,4);ax1N.set_yticklabels('');
rect = patches.Rectangle((1.01,.6),.6,.45,transform=ax0N.transAxes,edgecolor='grey', facecolor='white',clip_on=False)
ax0N.add_patch(rect)
blues = ['skyblue','steelblue','darkblue','darkblue','mediumturquoise']
# -------------
ax2 = fig.add_axes((2,0,.8,1.4),)
ax2.scatter(1,3.64,color=blues[1], marker='o',s=200, ) ##MARCATS
ax2.errorbar(1, 6.80,yerr = [[6.80-2.3], [8.8-6.8]], color=blues[0], marker='>',markersize=20, ) ## weber
ax2.scatter(1, 2.46,edgecolor=blues[0],color='none', marker='>',s=300, )## weber, diffusive only
ax2.scatter(2,11.02,color=blues[1], marker='o',s=200, label='MARCATS-CH4\n(diffusive)')
ax2.errorbar(2,7.85, yerr=[[7.85-2.50], [9.2-7.85]],color=blues[0], marker='>',markersize=20, label='Weber-CH4\n(diffusive+ebullitive)')
ax2.scatter(2, 3.06,edgecolor=blues[0],color='none', marker='>',s=300, label = 'Weber-CH4\n(diffusive)')## weber, diffusive only
rect = patches.Rectangle((1.01,.6),.65,.41,transform=ax2.transAxes,edgecolor='grey', facecolor='white',clip_on=False)
ax2.add_patch(rect)
ax2.set_title('Total CH$_4$ Flux')
ax2.set_xticks([1,2]); ax2.set_xticklabels(['Narrow coast', 'Wide coast', ],);
ax2.set_ylabel('Tg CH$_4$ yr$^{-1}$')
ax2.legend(loc=(1.02,.61), fontsize=15);ax2.text(1.02,.95,'CH$_4$ Products', transform=ax2.transAxes)
ax2.set_ylim(0,12); ax2.set_xlim(.5,2.5)
ax0C.text(-.15,1.01,'a',fontweight='bold',transform=ax0C.transAxes)
ax1N.text(-.15,1.01,'b',fontweight='bold',transform=ax1N.transAxes)
ax2.text(-.15,1.01,'c',fontweight='bold',transform=ax2.transAxes)
fig.savefig('./flux-co2_n2o_ch4.pdf', bbox_inches='tight', format='pdf')
# # MODELS
from scipy.stats import linregress
fluxtrend_thinshelf = np.zeros(11)
fluxtrend_wideshelf = np.zeros(11)
fluxtrend_open = np.zeros(11)
pco2trend_thinshelf = np.zeros(11)
pco2trend_wideshelf = np.zeros(11)
pco2trend_open = np.zeros(11)
for k in [0]: ##CESM-LR has only through the end of 2017 so we calculate trend 1998-2017, inclusive
fluxtrend_thinshelf[k], a,b,c,d = linregress(np.arange(240),mean_fco2_thinshelf[0:-12,k])
fluxtrend_wideshelf[k], a,b,c,d = linregress(np.arange(240),mean_fco2_wideshelf[0:-12,k])
fluxtrend_open[k], a,b,c,d = linregress(np.arange(240),mean_fco2_open[0:-12,k])
pco2trend_thinshelf[k], a,b,c,d = linregress(np.arange(240),mean_pco2_thinshelf[0:-12,k])
pco2trend_wideshelf[k], a,b,c,d = linregress(np.arange(240),mean_pco2_wideshelf[0:-12,k])
pco2trend_open[k], a,b,c,d = linregress(np.arange(240),mean_pco2_open[0:-12,k])
for k in range (1,11):## all others we calculate 1998-2018, inclusive
fluxtrend_thinshelf[k], a,b,c,d = linregress(np.arange(252),mean_fco2_thinshelf[:,k])
fluxtrend_wideshelf[k], a,b,c,d = linregress(np.arange(252),mean_fco2_wideshelf[:,k])
fluxtrend_open[k], a,b,c,d = linregress(np.arange(252),mean_fco2_open[:,k])
pco2trend_thinshelf[k], a,b,c,d = linregress(np.arange(252),mean_pco2_thinshelf[:,k])
pco2trend_wideshelf[k], a,b,c,d = linregress(np.arange(252),mean_pco2_wideshelf[:,k])
pco2trend_open[k], a,b,c,d = linregress(np.arange(252),mean_pco2_open[:,k])
fluxtrend_thinobs = np.zeros(2)
fluxtrend_wideobs = np.zeros(2)
fluxtrend_openobs = np.zeros(2)
pco2trend_thinobs = np.zeros(2)
pco2trend_wideobs = np.zeros(2)
pco2trend_openobs = np.zeros(2)
for k in range (2):
fluxtrend_thinobs[k], a,b,c,d = linregress(np.arange(252),mean_fco2obs_thinshelf[:,k])
fluxtrend_wideobs[k], a,b,c,d = linregress(np.arange(252),mean_fco2obs_wideshelf[:,k])
fluxtrend_openobs[k], a,b,c,d = linregress(np.arange(252),mean_fco2obs_open[:,k])
pco2trend_thinobs[k], a,b,c,d = linregress(np.arange(252),mean_pco2obs_thinshelf[:,k])
pco2trend_wideobs[k], a,b,c,d = linregress(np.arange(252),mean_pco2obs_wideshelf[:,k])
pco2trend_openobs[k], a,b,c,d = linregress(np.arange(252),mean_pco2obs_open[:,k])
blues = ['skyblue','steelblue','darkblue','darkblue','mediumturquoise']
fig, axs = plt.subplots(1,2, figsize=(20,8))
a=[np.zeros(8)+1e6,fluxtrend_wideshelf*120,fluxtrend_open*120]
vp=axs[1].violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
b=[fluxtrend_thinshelf[[1,4,6,7]]*120,fluxtrend_wideshelf[[1,4,6,7]]*120,fluxtrend_open[[1,4,6,7]]*120]
vp=axs[1].violinplot(b, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
axs[1].scatter(1,fluxtrend_thinobs[0]*120, color=blues[0], marker='>',s=200,label='Carboscope-1')
axs[1].scatter(1,fluxtrend_thinobs[1]*120, color=blues[1], marker='^',s=200,label='CMEMS*')
axs[1].scatter(2,fluxtrend_wideobs[0]*120, color=blues[0],marker='>',s=200)
axs[1].scatter(2,fluxtrend_wideobs[1]*120, color=blues[1], marker='^',s=200)
axs[1].scatter(3,fluxtrend_openobs[0]*120, color=blues[0],marker='>',s=200)
axs[1].scatter(3,fluxtrend_openobs[1]*120, color=blues[1], marker='^',s=200)
axs[1].plot(np.arange(5),np.zeros(5),color='k')
axs[1].set_title('Flux Density Trend 1998-2018')
axs[1].set_xticks([1,2,3]); axs[1].set_xticklabels(['Narrow Coast','Wide Coast', 'Open Ocean'],)# rotation=40);
axs[1].set_ylabel('[mol m$^{-2}$ yr$^{-1}$] decade$^{-1}$')
axs[1].legend(loc=(.5,.1), fontsize=13)
axs[1].set_ylim(-.22,.01); axs[1].set_xlim(.5,3.5)
a=[np.zeros(8),pco2trend_wideshelf*120, pco2trend_open*120]
vp=axs[0].violinplot(a, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('yellowgreen')
for k in [vp['cmedians'],]:
k.set_color('yellowgreen')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
b=[pco2trend_thinshelf[[1,4,6,7]]*120,pco2trend_wideshelf[[1,4,6,7]]*120, pco2trend_open[[1,4,6,7]]*120]
vp=axs[0].violinplot(b, showmedians=True,)
for pc in vp['bodies']:
pc.set_facecolor('goldenrod')
for k in [vp['cmedians'],]:
k.set_color('goldenrod')
k.set_linewidths(3)
for k in [vp['cbars'],vp['cmaxes'],vp['cmins']]: k.set_color('none')
axs[0].scatter(1,pco2trend_thinobs[0]*120, color=blues[0], marker='>',s=200, )
axs[0].scatter(1,pco2trend_thinobs[1]*120, color=blues[1], marker='^',s=200, )
axs[0].scatter(.7,13,marker='<', color='grey',s=200, label='Laruelle et al. 2018')
axs[0].scatter(.7,19.3, marker='o', color='grey', s=200, label='Wang et al. 2017')
axs[0].scatter(2,pco2trend_wideobs[0]*120, color=blues[0], marker='>',s=200)
axs[0].scatter(2,pco2trend_wideobs[1]*120, color=blues[1], marker='^',s=200)
axs[0].scatter(3,pco2trend_openobs[0]*120, color=blues[0], marker='>',s=200)
axs[0].scatter(3,pco2trend_openobs[1]*120, color=blues[1], marker='^',s=200)
axs[0].plot(np.arange(5),np.zeros(5)+trend_pco2_atm*120,color='k')
axs[0].set_title('pCO$_2$ Trend 1998-2018')
axs[0].set_xticks([1,2,3]); axs[0].set_xticklabels(['Narrow Coast','Wide Coast', 'Open Ocean'], )#rotation=40);
axs[0].set_ylabel('[uatm] decade$^{-1}$')
axs[0].set_ylim(12,21); axs[0].set_xlim(.5,3.5)
axs[0].text(1.7,20.3,'Atmospheric pCO$_2$ growth rate',fontsize=18)
axs[1].text(.55,-0.01,'Surface ocean matches atmospheric pCO$_2$ growth rate',fontsize=18)
axs[0].legend(loc=(2.22,.2), fontsize=16,edgecolor='none',facecolor='none',)
axs[1].legend(loc=(1.02,.4), fontsize=16,edgecolor='none',facecolor='none',)
ax1 = axs[1].twinx()
ax1.plot(np.arange(2),np.zeros(2),color='yellowgreen',linewidth=3,label='11 global models')
ax1.plot(np.arange(2),np.zeros(2),color='goldenrod',linewidth=3,label='4 hi-res \nglobal models')
ax1.legend(loc=(1.02,.6),edgecolor='none',facecolor='none',fontsize=16)
ax1.set_ylim(12,21);ax1.set_yticks([])
axs[1].text(1.02,.78,'Models',transform=axs[1].transAxes,fontsize=20)
axs[1].text(1.02,.55,'pCO$_2$ Products',transform=axs[1].transAxes,fontsize=20)
axs[1].text(1.02,.33,'Prior Estimates',transform=axs[1].transAxes,fontsize=20)
rect=patches.Rectangle((1.01,.18),.45,.68,edgecolor='grey', facecolor='none',clip_on=False,alpha=.6, transform=axs[1].transAxes)
axs[1].add_patch(rect)
axs[0].text(-.1,1.03,'c',fontweight='bold',transform=axs[0].transAxes,fontsize=20)
axs[1].text(-.1,1.03,'d',fontweight='bold',transform=axs[1].transAxes,fontsize=20)
fig.savefig('./trends_violins.pdf', bbox_inches='tight', format='pdf')
plt.rc('font', size = 16)
import matplotlib as mpl
mpl.rc('hatch', color='k', linewidth=5)
fig, axes = plt.subplots(3,1, figsize=(10,15), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
# ----------------------------
# # pCO2
with xr.open_dataset('pco2_obsproducts_filledcmems_trend.nc',decode_cf=True,) as ds:
trend_pco2_obs = ds.pco2[0,:,:,0:2]
with xr.open_dataset('trend_pco2_1998-2018.nc', decode_cf=True,) as ds: ##models
trend_pco2 = ds.spco2[0,:,:,[0,1,2,3,4,6,7,8,9,10,11]] #vestige of prior model organization
pc = trend_pco2 * 120
med_pco2 = xr.DataArray(np.nanmedian(pc, axis=2), coords = [pc.latitude, pc.longitude], dims=['latitude', 'longitude'])
trend_atm_dec = trend_pco2_atm*120
axes[0].coastlines()
gl=axes[0].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(trend_pco2_obs[:,:,0]*120).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[uatm] decade$^{-1}$'}, levels=np.arange(trend_atm_dec-10,trend_atm_dec+10.1), cmap='YlGnBu')
axes[0].set_title('Trend in pCO$_2$, 1998-2018\nCarboscope-1')
axes[0].text(0,.9,'a', transform=axes[0].transAxes, fontsize=20, fontweight = 'bold')
axes[0].set_global()
axes[1].coastlines()
gl=axes[1].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(trend_pco2_obs[:,:,1].sel(latitude = slice(-90,75))*120).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[uatm] decade$^{-1}$',}, levels=np.arange(trend_atm_dec-10,trend_atm_dec+10.1), cmap='YlGnBu')
axes[1].set_title('CMEMS*')
axes[1].text(0,.9,'b', transform=axes[1].transAxes, fontsize=20, fontweight = 'bold')
axes[1].set_global()
axes[2].coastlines()
gl=axes[2].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(med_pco2.where(~np.isnan(wide_mask))).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[uatm] decade$^{-1}$',}, levels=np.arange(trend_atm_dec-10,trend_atm_dec+10.1), cmap='YlGnBu')
axes[2].set_title('Model Median')
axes[2].text(0,.9,'c', transform=axes[2].transAxes, fontsize=20, fontweight = 'bold')
axes[2].set_global()
fig.tight_layout()
plt.savefig('pco2_trend_map.pdf', format='pdf', dpi=300, bbox_inches='tight')
plt.rc('font', size = 16)
import matplotlib as mpl
mpl.rc('hatch', color='k', linewidth=5)
fig, axes = plt.subplots(3,1, figsize=(10,15), squeeze=True, subplot_kw={'projection':ccrs.Robinson(central_longitude=370)})
# ----------------------------
# # FLUX
with xr.open_dataset('fgco2_obsproducts_filledcmems_trend.nc',decode_cf=True,) as ds:
trend_fgco2_obs = ds.fgco2_obs[0,:,:,0:2]
with xr.open_dataset('/tigress/GEOCLIM/LRGROUP/RECCAP2_coastal_ocean_repo/regrid/trend_fgco2_1998-2018.nc', decode_cf=True,) as ds:
trend_fgco2 = ds.fgco2[0,:,:,[0,1,2,3,4,6,7,8,9,10,11]]
flx = trend_fgco2 * 86400*365 * 120
med_flx = xr.DataArray(np.nanmedian(flx, axis=2), coords = [trend_fgco2.latitude, trend_fgco2.longitude], dims=['latitude', 'longitude'])
trend_disagree_modelmed = (med_pco2*0+1).where((med_pco2-trend_atm_dec)*med_flx >0,other=-1).where(~np.isnan(wide_mask))
# # trend_fgco2_obs
# trend_disagree_obs = (trend_pco2_obs*0+1).where((trend_pco2_obs*120-trend_atm_dec)*trend_fgco2_obs >0,other=-1).where(~np.isnan(wide_mask))
axes[0].coastlines()
gl=axes[0].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(trend_fgco2_obs[:,:,0]*120).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[mol m$^{-2}$ yr$^{-1}$] decade$^{-1}$'}, levels=np.arange(-.5,.55,.05), cmap='RdBu_r')
(trend_disagree_obs[:,:,0]).plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(),colors='none',hatches=['\\\\',''],levels=[-10,0,10],add_colorbar=False)
axes[0].set_title('Trend in Flux Density, 1998-2018\nCarboscope-1')
axes[0].text(0,.9,'a', transform=axes[0].transAxes, fontsize=20, fontweight = 'bold')
axes[0].set_global()
axes[1].coastlines()
gl=axes[1].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
(trend_fgco2_obs[:,:,1]*120).sel(latitude = slice(-90,75)).where(~np.isnan(wide_mask)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[mol m$^{-2}$ yr$^{-1}$] decade$^{-1}$',}, levels=np.arange(-.5,.55,.05), cmap='RdBu_r')
(trend_disagree_obs[:,:,1]).sel(latitude = slice(-90,75)).plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(),colors='none',hatches=['\\\\',''],levels=[-10,0,10],add_colorbar=False)
axes[1].set_title('Trend in Flux Density, 1998-2018\nCMEMS*')
axes[1].text(0,.9,'b', transform=axes[1].transAxes, fontsize=20, fontweight = 'bold')
axes[1].set_global()
axes[2].coastlines()
gl=axes[2].gridlines(draw_labels=True,)
gl.rotate_labels=False
gl.ylabel_style = {'size': 15, 'color': 'gray'}
gl.xlabel_style = {'size': 15, 'color': 'gray'}
med_flx.where(~np.isnan(wide_mask)).plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6, 'label':'[mol m$^{-2}$ yr$^{-1}$] decade$^{-1}$',}, levels=np.arange(-.5,.55,.05), cmap='RdBu_r')
trend_disagree_modelmed.plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(),colors='none',hatches=['\\\\',''],levels=[-10,0,10],add_colorbar=False)
axes[2].set_title('Trend in Flux Density, 1998-2018\nModel Median')
axes[2].text(0,.9,'c', transform=axes[2].transAxes, fontsize=20, fontweight = 'bold')
axes[2].set_global()
fig.tight_layout()
plt.savefig('flux_trend_map.pdf', format='pdf', dpi=300, bbox_inches='tight')
# # # # # START // LOAD DATA
# a_obs = area_wide * (1-seaice_obs)
# a_obs.load()
# fgco2_obs.load()
# fco2_coastal_mom6kw['time'] = a_obs.time[13*12:-48]
# lat_mean_coastal = (fgco2_obs[2,:]*a_obs/a_obs.sum('longitude')) .sum('longitude')
# lat_mean_coastal_momkw = (fco2_alt_kw*a_obs/a_obs.sum('longitude')) .sum('longitude')
mean_fgco2_latitude_clim = lat_mean_coastal.groupby('time.month').mean('time')
djf_fgco2_latitude = mean_fgco2_latitude_clim[[0,1,11],:].mean('month')
jja_fgco2_latitude = mean_fgco2_latitude_clim[5:8,:].mean('month')
son_fgco2_latitude = mean_fgco2_latitude_clim[8:11,:].mean('month')
mam_fgco2_latitude = mean_fgco2_latitude_clim[2:5,:].mean('month')
# # # # # #
mean_fgco2obs_latitude_clim = lat_mean_coastal_momkw.groupby('time.month').mean('time')
djf_fgco2obs_latitude = mean_fgco2obs_latitude_clim[[0,1,11],:].mean('month')
jja_fgco2obs_latitude = mean_fgco2obs_latitude_clim[5:8,:].mean('month')
son_fgco2obs_latitude = mean_fgco2obs_latitude_clim[8:11,:].mean('month')
mam_fgco2obs_latitude = mean_fgco2obs_latitude_clim[2:5,:].mean('month')
# # # # END // LOAD DATA
# # # #
fig,axs=plt.subplots(3,2, figsize=(20,15))
mean_fgco2_latitude_clim.mean('month').plot( ax = axs[0][0], color='darkblue', label = 'Coastal SOM-FFN')
mean_fgco2obs_latitude_clim.mean('month').plot(ax = axs[0][0], color='mediumturquoise', label = 'Coastal SOM-FFN-kw' )
axs[0][0].set_title('Annual mean flux density')
axs[0][0].legend(loc=(.1,.01),fontsize=20)
axs[0][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[0][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[0][0].set_xlim(-80,90)
axs[0][0].legend()
(djf_fgco2_latitude-jja_fgco2_latitude).plot(ax = axs[0][1], color='darkblue', )
(djf_fgco2obs_latitude-jja_fgco2obs_latitude).plot( ax = axs[0][1], color='mediumturquoise',)
axs[0][1].set_title('Seasonal amplitude (DJF - JJA)')
axs[0][1].legend(loc=(1,.9))
axs[0][1].plot(np.arange(-80,80),np.zeros(160),'k')
axs[0][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[0][1].set_xlim(-80,90)
mam_fgco2_latitude.plot( ax = axs[1][0], color='darkblue', )
mam_fgco2obs_latitude.plot( ax = axs[1][0], color='mediumturquoise', )
axs[1][0].set_title('MAM flux density')
axs[1][0].legend(loc=(1,.9))
axs[1][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[1][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[1][0].text(65,4,'spring');axs[1][0].text(-75,4,'fall')
axs[1][0].set_xlim(-80,90)
jja_fgco2_latitude.plot( ax = axs[1][1], color='darkblue', )
jja_fgco2obs_latitude.plot(ax = axs[1][1], color='mediumturquoise', )
axs[1][1].set_title('JJA flux density')
axs[1][1].legend(loc=(1,.9))
axs[1][1].plot(np.arange(-100,100),np.zeros(200),'k')
axs[1][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[1][1].text(65,4,'summer');axs[1][1].text(-75,4,'winter')
axs[1][1].set_xlim(-80,90)
son_fgco2_latitude.plot( ax = axs[2][0], color='darkblue', )
son_fgco2obs_latitude.plot(ax = axs[2][0], color='mediumturquoise',)
axs[2][0].set_title('SON flux density')
axs[2][0].legend(loc=(1,.9))
axs[2][0].plot(np.arange(-100,100),np.zeros(200),'k')
axs[2][0].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[2][0].text(65,4,'fall');axs[2][0].text(-75,4,'spring')
axs[2][0].set_xlim(-80,90)
djf_fgco2_latitude.plot(ax = axs[2][1], color='darkblue',)
djf_fgco2obs_latitude.plot(ax = axs[2][1], color='mediumturquoise',)
axs[2][1].set_title('DJF flux density')
axs[2][1].legend()
axs[2][1].plot(np.arange(-100,100),np.zeros(200),'k')
axs[2][1].set_ylabel('mol m$^{-2}$ yr$^{-1}$')
axs[2][1].text(65,4,'winter'); axs[2][1].text(-75,4,'summer')
axs[2][1].set_xlim(-80,90)
axs[0][0].set_ylim(-7,5)
axs[0][1].set_ylim(-7,5)
axs[1][0].set_ylim(-7,5)
axs[1][1].set_ylim(-7,5)
axs[2][0].set_ylim(-7,5)
axs[2][1].set_ylim(-7,5)
axs[0][0].axvline(0, color='k',linewidth=.5)
axs[0][1].axvline(0, color='k',linewidth=.5)
axs[1][0].axvline(0, color='k',linewidth=.5)
axs[1][1].axvline(0, color='k',linewidth=.5)
axs[2][0].axvline(0, color='k',linewidth=.5)
axs[2][1].axvline(0, color='k',linewidth=.5)
axs[0][1].text(45,3,'Stronger sink in \nsummer')
axs[0][1].text(45,-6.5,'Stronger sink in \nwinter')
axs[0][1].text(-78,3,'Stronger sink in \nwinter')
axs[0][1].text(-78,-6.5,'Stronger sink in \nsummer')
axs[0][0].text(-80,5.2,'a', fontweight='bold')
axs[0][1].text(-80,5.2,'b', fontweight='bold')
axs[1][0].text(-80,5.2,'c', fontweight='bold')
axs[1][1].text(-80,5.2,'d', fontweight='bold')
axs[2][0].text(-80,5.2,'e', fontweight='bold')
axs[2][1].text(-80,5.2,'f', fontweight='bold')
fig.tight_layout()
fig.savefig('mean_pco2_latitude_coastal_vs_coastalkw.pdf', format='pdf', dpi=300, bbox_inches='tight')