In [94]:
%matplotlib inline
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import math
#import exceptions
# Convert to py file using ipython: !jupyter nbconvert --to script C_budget_per_region.ipynb
In [102]:
fpathRegions = '/users/tams00/Documents/nerc_ssb/c_fluxes/analysis/amm7_subdomains_v2.nc'
shppath = '/users/tams00/Documents/nerc_ssb/c_fluxes/shelf_regions/gebco_200_clip.shp'
outdir = '/users/tams00/Documents/nerc_ssb/c_fluxes/analysis/v0.2concat/'
#fpath = '/users/tams00/Documents/nerc_ssb/c_fluxes/analysis/v0.2concat/DIC_advection_hind_1990_2015.nc'
#varname='DIC_advection_hind'
#name= 'Advection of Inorganic C'
#fpath = '/users/tams00/Documents/nerc_ssb/c_fluxes/analysis/v0.2concat/Organic_C_advection_1990_2015.nc'
#varname='Organic_C_advection'
#name = 'Advection of Organic C'
#fpath = '/users/tams00/Documents/nerc_ssb/c_fluxes/Yuri_ftp/SSB_Cpaper/hindcast/amm7_hindcast_1990_2015_IC.nc'
#yuriDS = xr.open_dataset(fpath)
#varname='IC'
#name='Depth integrated pelagic inorganic carbon'
fpath = '/users/tams00/Documents/nerc_ssb/c_fluxes/Yuri_ftp/SSB_Cpaper/hindcast/amm7_hindcast_1990_2015_benIC.nc'
yuriDS = xr.open_dataset(fpath)
varname='benIC'
name='Benthic inorganic carbon'
#fpath = 'SSB_Cpaper/hindcast/amm7_hindcast_1990_2015_benOC.nc'
#yuriDS = xr.open_dataset(fpath)
#quantDA = yuriDS['burial']
#name=
#quantDA = yuriDS['benPOC']
#name=
#fpath = 'SSB_Cpaper/hindcast/amm7_hindcast_1990_2015_pelbenc.nc'
#varname='pelbenc'
#name='Net benthic to pelagic C flux'
#fpath='../c_fluxes/analysis/v0.2concat/DIC_release_from_benthic_1990_2015.nc'
#varname='DIC_release_from_benthic'
#name = varname.replace('_',' ')
#fpath='../c_fluxes/analysis/v0.2concat/POC_pelagic_to_benthic_1990_2015.nc'
#varname='POC_pelagic_to_benthic'
#name = varname.replace('_',' ')
#fpath = '../c_fluxes/analysis/v0.2/air_sea_flux_1981_1999.nc'
#varname='air_sea_flux'
#name = varname.replace('_',' ')
#fpath = 'SSB_Cpaper/hindcast/amm7_hindcast_1990_2015_airsea.nc'
#varname='O3_fair'
#name = 'Air-sea flux'
DS = xr.open_dataset(fpath)
print(DS)
DS.coords['latitude'] = DS['nav_lat']
DS.coords['longitude'] = DS['nav_lon']
print(DS)
quantDA = DS[varname]
figname = varname+'_1990-2015.png'
print(quantDA.units)
print(quantDA.dims)
In [103]:
(quantDA.mean(dim='time_counter',keep_attrs=True)).plot(robust=True)
Out[103]:
In [97]:
regionsDS = xr.open_dataset(fpathRegions)
regionsDS.rename({'lat':'nav_lat','lon':'nav_lon'},inplace=True)
regionsDS = regionsDS.drop(['Bathymetry','Ocean','DeepOcean','SouthernNEA','NorthernNEA'])
In [98]:
lon2d,lat2d = np.meshgrid(regionsDS.nav_lon,regionsDS.nav_lat)
EARTHRADIUS = 6371000 # https://en.wikipedia.org/wiki/Earth_radius
delNS = np.zeros_like(lon2d)
delNS[1:,:] = lat2d[1:,:] - lat2d[0:-1,:]
delNS[0,:] = delNS[1,:]
delNS = delNS * 2 * math.pi * EARTHRADIUS / 360
delEW = np.zeros_like(lon2d)
delEW[:,1:] = lon2d[:,1:] - lon2d[:,0:-1]
delEW[:,0] = delEW[:,1]
delEW = delEW * 2 * math.pi * EARTHRADIUS / 360 * np.cos(lat2d/180*math.pi)
area = xr.DataArray(delEW * delNS, coords=regionsDS['Shelf'].coords)
In [104]:
molarmC = 12.0107
#print(quantDA)
if 'mgC' in quantDA.units and '/day' in quantDA.units:
print('converting unit from mgC/m2/day to molC/m2/yr)')
quantDA.values = quantDA/molarmC/1000*365
quantDA.attrs['units'] = quantDA.units.replace('mgC','molC')
quantDA.attrs['units'] = quantDA.units.replace('/day','/yr')
elif 'mol' in quantDA.units and '/day' in quantDA.units:
print('converting unit from mmol/m2/day to molC/m2/yr)')
quantDA.values = quantDA/1000*365
quantDA.attrs['units'] = quantDA.units.replace('mmol','mol')
if 'm^3' in quantDA.units: # error in Yuri's stock files?
quantDA.attrs['units'] = quantDA.units.replace('m^3','m^2')
print(quantDA.units)
# Yuri stocks (e.g. IC)
elif 'mol' in quantDA.units and not '/day' in quantDA.units:
print('converting unit from mmol/m2 to mol/m2)')
quantDA.values = quantDA/1000
quantDA.attrs['units'] = quantDA.units.replace('mmol','mol')
if 'm^3' in quantDA.units: # error in Yuri's stock files?
quantDA.attrs['units'] = quantDA.units.replace('m^3','m^2')
print(quantDA.units)
elif 'gC' in quantDA.units and '/day' in quantDA.units:
print('converting unit from gC/m2/day to molC/m2/yr)')
quantDA.values = quantDA/molarmC*365
quantDA.attrs['units'] = quantDA.units.replace('gC','molC')
quantDA.attrs['units'] = quantDA.units.replace('/day','/yr')
else:
print('NOTE: Check units:',quantDA.units )
quant = quantDA.mean(dim='time_counter').values
print( '{0:15} {1:^12} {2:^15} {3:^15} {4:^15} {5:^15} {6:^15}'.format(
'Regions','Area','Average','Total','Std. dev.','Total','Std. dev.' ) )
if '/d' in quantDA.units:
print( '{0:15} {1:^12} {2:^15} {3:^15} {4:^15} {5:^15} {6:^15}'.format(
'','(km2)','molC/m2/yr','molC/yr','molC/yr','gC/yr', 'gC/yr' ) )
else:
print( '{0:15} {1:^12} {2:^15} {3:^15} {4:^15} {5:^15} {6:^15}'.format(
'', '(km2)','molC/m2','molC','molC', 'gC', 'gC' ) )
# Data in molC
for region in regionsDS.data_vars:
mask = regionsDS[region]
areakm2 = float((mask * area ).sum()/1e6)
maskarea = np.expand_dims( (mask*area).values,axis=0 )
try:
ts = (maskarea*quantDA).sum(['x_grid_T','y_grid_T'])
except ValueError:
# Yuris files have different coordinates
ts = (maskarea*quantDA).sum(['x','y'])
mean = float(xr.DataArray.mean(ts))
stddev = float(xr.DataArray.std(ts.groupby('time_counter.year').mean('time_counter')))
print( '{0:15} {1:^12.2e} {2:^15.2e} {3:^15.2e} {4:^15.2e} {5:^15.2e} {6:^15.2e}'.format(
region,areakm2,mean/(areakm2*1e6),mean,stddev,mean*molarmC,stddev*molarmC) )
In [101]:
region = 'Shelf'
mask = regionsDS[region]
maskarea = np.expand_dims( (mask*area).values,axis=0 )
fig, axes = plt.subplots(ncols=2,figsize=(14,4))
try:
# monthly means # CHECK FACTORS!!!!!!!!!!!!!!!!!
( (quantDA*maskarea).sum(['x','y'])/1e12 ).plot(ax=axes[0])
# yearly means
( (quantDA*maskarea).groupby('time_counter.year').mean('time_counter').sum(['x','y'])/1e12 ).plot(ax=axes[1])
for i, ax in enumerate(axes.flat):
ax.set_title(name+': '+region)
ax.set_xlabel('year')
ax.set_ylabel(quantDA.units)
except ValueError:
#except Exception as e: # exceptions.ValueError:
#print('Exception is ',e)
# monthly means # CHECK FACTORS!!!!!!!!!!!!!!!!!
( (quantDA*maskarea).sum(['x_grid_T','y_grid_T'])/1e12 ).plot(ax=axes[0])
# yearly means
( (quantDA*maskarea).groupby('time_counter.year').mean('time_counter').sum(['x_grid_T','y_grid_T'])/1e12 ).plot(ax=axes[1])
for i, ax in enumerate(axes.flat):
ax.set_title(name+': '+region)
ax.set_xlabel('year')
ax.set_ylabel('T molC/yr')
figname_ts = figname.split('.')[0]+'_ts.png'
plt.savefig(outdir+figname_ts,dpi=200)
In [83]:
region = 'Shelf'
mask = regionsDS[region]
maskarea = np.expand_dims( (mask*area).values,axis=0 )
try:
( (quantDA*maskarea).groupby('time_counter.month').mean('time_counter').sum(['x','y']) ).plot()
except Exception as e:
print('Exception is',e)
( (quantDA*maskarea).groupby('time_counter.month').mean('time_counter').sum(['x_grid_T','y_grid_T'])/1e12 ).plot()
plt.title(name+': '+region)
#plt.ylabel(quantDA.units)
plt.ylabel('T molC/yr')
figname_ts = figname.split('.')[0]+'_climatology.png'
plt.savefig(outdir+figname_ts,dpi=200)
In [87]:
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
#import cartopy.io.shapereader as shpreader
import shapefile
#import exceptions
DPI = 150
#****************************************************
# create map, safe to file and display figure
#****************************************************
def create_map(varDA,figtitle,box,limits=[],shppath=[],figname=[]):
plt.figure(figsize=(10,10))
ax = plt.axes(projection=ccrs.PlateCarree())
for coordname in varDA.coords:
if ('lat' in coordname):
ycoordname=coordname
elif ('lon' in coordname):
xcoordname=coordname
if not(xcoordname and ycoordname):
print('varDA.coords:\t',varDA.coords)
sys.exit('Could not find latitude and longitude in coordinates:')
if limits:
im = varDA.plot.pcolormesh(ax=ax,x=xcoordname,y=ycoordname, \
cmap=plt.get_cmap('RdBu_r'), \
vmin=limits[0],vmax=limits[1])
#levels=np.linspace(limits[0],limits[1]),
else:
im = varDA.plot.pcolormesh(ax=ax,x=xcoordname,y=ycoordname, \
cmap=plt.get_cmap('RdYlBu_r'))
if box:
plt.axis(box)
if shppath:
reader = shapefile.Reader(shppath)
for shape in reader.shapes():
pts_arr=np.array(shape.points)
ax.plot(pts_arr[:,0],pts_arr[:,1],color='black',linestyle=':')#,linewidth=0.5)
try:
im.colorbar.set_label(varDA.units)
except AttributeError:
#except KeyError:
#except exceptions.KeyError:
print('No attribute unit in variable to plot')
im.colorbar.set_label('gC/m2/Yr')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle=':')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
ax.coastlines('50m')
ax.set_title(figtitle)
if figname:
fname = figname
else:
fname=figtitle.replace(' ','_')+'.png'
print('saving to ',fname)
plt.savefig(outdir+fname,dpi=DPI,bbox_inches='tight')
#plt.show()
In [92]:
#box=[-20,12,40,65]
box=[-18,12,41,64]
#print(quantDA)
#for Fluxes
#create_map(quantDA.mean(dim='time_counter',keep_attrs=True),figtitle = name,figname=figname,limits=[-25,25],box=box,shppath=shppath)
# for Benthic IC
create_map(quantDA.mean(dim='time_counter',keep_attrs=True),figtitle = name,figname=figname,limits=[0.2,0.35],box=box,shppath=shppath)
#create_map(quantDA.isel(time_counter=311)/30./1000*365,figtitle='Carbon burial',limits=[],box=box,shppath=shppath)
In [86]:
DIC = quantDA
In [128]:
OC = quantDA
In [134]:
box=[-18,12,42,63]
create_map((OC.mean(dim='time_counter',keep_attrs=True)+DIC.mean(dim='time_counter',keep_attrs=True))/1000*365,figtitle = "Adv OC_Adv DIC",figname="Adv_OC_Adv_DIC",limits=[-25,25],box=box,shppath=shppath)