<jupyter_start><jupyter_text>Figure of the base case's bottom concentration<jupyter_code>import cmocean as cmo
import matplotlib.pyplot as plt
import matplotlib.gridspec as gspec
import matplotlib as mpl
%matplotlib inline
from netCDF4 import Dataset
import numpy as np
import seaborn as sns
import xarray as xr
import canyon_tools.readout_tools as rout
import canyon_tools.savitzky_golay as sg
def calc_rho(RhoRef,T,S,alpha=2.0E-4, beta=7.4E-4):
    """-----------------------------------------------------------------------------
    calc_rho calculates the density using a linear equation of state.
    
    INPUT:
    RhoRef : reference density at the same z as T and S slices. Can be a scalar or a 
             vector, depending on the size of T and S.
    T, S   : should be at least 2D arrays in coordinate order (..., Y , X ) 
    alpha = 2.0E-4 # 1/degC, thermal expansion coefficient
    beta = 7.4E-4, haline expansion coefficient
    OUTPUT:
    rho - Density [...,ny,nx]
    -----------------------------------------------------------------------------"""
    
    #Linear eq. of state 
    rho = RhoRef*(np.ones(np.shape(T)) - alpha*(T[...,:,:]) + beta*(S[...,:,:]))
    return rho

def call_unstag(t):
    UU,VV = rout.unstagger(state.U.isel(T=t),state.V.isel(T=t))
    return(UU,VV)


def call_rho(t,yslice,xslice):
    T = state.Temp.isel(T=t,Y=yslice,X=xslice)
    S = state.S.isel(T=t,Y=yslice,X=xslice)
    rho = calc_rho(RhoRef,T,S,alpha=2.0E-4, beta=7.4E-4)
    return(rho) 

def mask2DCanyon(bathy, sbdepth=-152.5):
    '''Mask out the canyon from the shelf.
    bathy : depths 2D array from the grid file
    sbdepth: shelf depth, always negative float 
    Returns mask'''
    
    bathyMasked = np.ma.masked_less(-bathy, -152.5)
    return(bathyMasked.mask)

def ConcAreaShelf(Tr, hfac, ra, bathy, sbdepth=-152.5):
    '''Tr: Tracer field (nt,nz,ny,nx)
       hfac: fraction of open cell at center (nz,ny,nx)
       ra: array of cell horizontal areas (ny,nx)
       bathy : depths 2D array from the grid file (ny,nx)
       sbdepth: shelf break depth (negative value)
       
       RETURNS:
       ConcArea = concentration at cell closest to bottom times its area (nt,ny,nx). 
       Conc = cocetration near bottom (nt,ny,nx)
       Area = 2D array of cell areas on shelf (ny,nx)
       NOTE: Everything below the canyon rim is masked out'''
    
    ConcArea = np.empty((19,360,360))
    Conc = np.empty((19,360,360))
    ConcFiltered = np.empty((19,360,360))
    Area = np.empty((19,360,360))
    
    BottomInd = np.argmax(hfac[::-1,:,:]>0.0,axis=0) # start looking for first no-land cell from the bottom up.
    BottomInd = np.ones(np.shape(BottomInd))*89 - BottomInd # Get index of unreversed z axis
    
    print(np.shape(BottomInd))
    for tt in range(19):
        for j in range(360):
            for i in range(360):
                TrBottom = Tr[tt,BottomInd[i,j],i,j]
                if TrBottom > 0.0:
                    ConcArea[tt,i,j] = TrBottom*ra[i,j]
                    Conc[tt,i,j] = TrBottom
                    Area[tt,i,j] = ra[i,j]
                else:
                    ConcArea[tt,i,j] = np.NaN
                    Conc[tt,i,j] = np.NaN
                    Area[tt,i,j] = np.NaN
                
            # Filter step noise
            ConcFiltered[tt,:,j] = sg.savitzky_golay(Conc[tt,:,j], 7,3) 
                
    print(np.shape(ConcArea))
    
    maskShelf = np.ma.masked_invalid(ConcFiltered) 
    
    return (np.ma.masked_invalid(ConcArea),
            np.ma.masked_invalid(ConcFiltered),
            np.ma.masked_invalid(Area),
            BottomInd)
# Cross-shelf velocity, Plot1
def Plot1(gs_ax,Tr_slice,zslice,yind,xslice,Trmin=0,Trmax=11): 
    '''gs_ax: grid handle
       Tr_slice: tracer concentration sliced in x-z (nt,nz,nx)
       tslice: advective phase slice to average over.
       yind: cross-shore index of the plot
    '''
    ax = plt.subplot(gs_ax)
    
    csU = np.linspace(Trmin,Trmax,num=31)
    csU2 = np.linspace(Trmin,Trmax,num=16)
    
    Trplot = np.ma.array(np.nanmean(Tr_slice.isel(T=tslice).data,0),mask=MaskC[:,yind])
    mesh = ax.contourf(grid.X[xslice]/1000,grid.Z[zslice],Trplot[zslice,xslice],
                       csU,
                       cmap=cmap_name)
    

    CS = ax.contour(grid.X[xslice]/1000,grid.Z[zslice],
                    Trplot[zslice,xslice],
                    csU2,
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    ax.set_ylabel('Depth (m)',labelpad=0.3)
    ax.text(0.6,0.02,'Tracer Concentration ($\mu$mol/L)',transform=ax.transAxes)
    
# Cross-shelf velocity, Plot1
def Plot1anom(gs_ax,Tr_slice,Tr_sliceNoC,tslice,zslice,yind,xslice,Trmin=0,Trmax=11): 
    '''gs_ax: grid handle
       Tr_slice: tracer concentration sliced in x-z (nt,nz,nx)
       tslice: advective phase slice to average over.
       yind: cross-shore index of the plot
    '''
    ax = plt.subplot(gs_ax)
    
    csU = np.linspace(Trmin,Trmax,num=31)
    csU2 = np.linspace(Trmin,Trmax,num=16)
    
    Trplot = np.ma.array(np.nanmean(Tr_slice.isel(T=tslice).data-Tr_sliceNoC.isel(T=tslice).data,0),
                         mask=MaskCNoC[:,yind,:])
    mesh = ax.contourf(grid.X[xslice]/1000,grid.Z[zslice],Trplot[zslice,xslice],
                       np.linspace(-4,4,17),
                       cmap=cmap_name)
    

    CS = ax.contour(grid.X[xslice]/1000,grid.Z[zslice],
                    Trplot[zslice,xslice],
                    np.linspace(-4,4,9),
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    ax.axhline(y=-95,linestyle='--',color='k')
    
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    ax.set_ylabel('Depth (m)',labelpad=0.3)
    ax.text(0.5,0.02,'Tracer concentration anomaly ($\mu$mol/L)',transform=ax.transAxes)
    
# Cross-shelf velocity, Plot1
def Plot2(gs_ax,Tr_slice,yslice,zslice,xind,Trmin=0,Trmax=11): 
    '''gs_ax: grid handle
       Tr_slice: tracer concentration sliced in y-z (nt,nz,ny)
       tslice: advective phase slice to average over.
       xind: cross-shore index of the plot
    '''
    ax = plt.subplot(gs_ax)
    
    csU = np.linspace(Trmin,Trmax,num=31)
    csU2 = np.linspace(Trmin,Trmax,num=16)
    
    Trplot = np.ma.array(np.nanmean(Tr_slice.isel(T=tslice).data,0),mask=MaskC[:,:,xind])
    mesh = ax.contourf(grid.Y[yslice]/1000,grid.Z[zslice],(Trplot[zslice,yslice]),
                       csU,
                       cmap=cmap_name)
    
    
    CS = ax.contour(grid.Y[yslice]/1000,grid.Z[zslice],
                    (Trplot[zslice,yslice]),
                    csU2,
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    ax.set_ylabel('Depth (m)',labelpad=0.3)
    ax.set_xlabel('Cross-shelf distnce (km)',labelpad=0.3)
    
    ax.text(0.4,0.02,'Tracer concentration ($\mu$mol/L)',transform=ax.transAxes)
    
def Plot2anom(gs_ax,Tr_slice,Tr_sliceNoC,tslice,yslice,zslice,xind,Trmin=0,Trmax=11): 
    '''gs_ax: grid handle
       Tr_slice: tracer concentration sliced in y-z (nt,nz,ny)
       tslice: advective phase slice to average over.
       xind: cross-shore index of the plot
    '''
    ax = plt.subplot(gs_ax)
    
    
    Trplot = np.ma.array(np.nanmean(Tr_slice.isel(T=tslice).data-Tr_sliceNoC.isel(T=tslice).data,0),
                         mask=MaskCNoC[:,:,xind])
    mesh = ax.contourf(grid.Y[yslice]/1000,grid.Z[zslice],(Trplot[zslice,yslice]),
                       np.linspace(-4,4,17),
                       cmap=cmap_name)
    
     
    CS = ax.contour(grid.Y[yslice]/1000,grid.Z[zslice],
                    (Trplot[zslice,yslice]),
                    np.linspace(-4,4,9),
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    ax.set_ylabel('Depth (m)',labelpad=0.3)
    ax.set_xlabel('Cross-shelf distnce (km)',labelpad=0.3)
    
    ax.text(0.2,0.02,'Tracer concentration Anomaly',transform=ax.transAxes)
def Plot3(gs_ax,ConcFilt,tslice,xslice,yslice,xind,yind):
    
    ax = plt.subplot(gs_ax)
    BAC = np.nanmean(ConcFilt[tslice,yslice,xslice],axis=0)
    
    maxCM = np.max(BAC)
    minCM = np.min(BAC)
    
    mesh = ax.contourf(grid.X[xslice]/1000,grid.Y[yslice]/1000,BAC,21, 
                      vmax=maxCM, 
                      vmin=minCM,
                      cmap=cmap_name)
    CS = ax.contour(grid.X[xslice]/1000,grid.Y[yslice]/1000,BAC,
                    11,
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    #cax,kw = mpl.colorbar.make_axes([ax],location='right',anchor=(5,0),shrink=1.1)
    #cb = plt.colorbar(mesh, cax=cax,ticks=[np.linspace(minCM, maxCM,8) ],format='%.2f',**kw)
       
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    
    ax.text(0.6,0.8,'Bottom Concentration ($\mu$mol/L)',transform=ax.transAxes)
    SB = ax.contour(grid.X[xslice]/1000,grid.Y[yslice]/1000,
                     grid.Depth[yslice,xslice],
                    [50,75,100,125,150,200,250,300,350,400,450,500],
                    colors=['0.5','0.5','0.5','0.5','0.1','0.5','0.5','0.5','0.5','0.5','0.5','0.5'],
                    linewidths=[0.75] )
    
    ax.axhline(grid.Y[yind]/1000,color='k',linestyle='--')
    ax.axvline(grid.X[xind]/1000,color='k',linestyle='--')
    
    ax.set_xlabel('Alongshore distance (km)',labelpad=0.3)
    ax.set_ylabel('Cross-shore distance (km)',labelpad=0.3)
def Plot3anom(gs_ax,ConcFilt,tslice,xslice,yslice,xind,yind,minCM,maxCM):
    
    ax = plt.subplot(gs_ax)
    BAC = np.nanmean(ConcFilt[tslice,yslice,xslice],axis=0)
    
    
    mesh = ax.contourf(grid.X[xslice]/1000,grid.Y[yslice]/1000,BAC, 
                      np.linspace(-4,4,17),
                      cmap=cmap_name)
    CS = ax.contour(grid.X[xslice]/1000,grid.Y[yslice]/1000,BAC,
                    np.linspace(-4,4,9),
                    colors='k',
                    linewidths=[0.75] )
    plt.clabel(CS,inline=True, fmt='%.1f')
    
    #cax,kw = mpl.colorbar.make_axes([ax],location='right',anchor=(5,0),shrink=1.1)
    #cb = plt.colorbar(mesh, cax=cax,ticks=[np.linspace(minCM, maxCM,8) ],format='%.2f',**kw)
       
    ax.set_axis_bgcolor((205/255.0, 201/255.0, 201/255.0))
    
    ax.text(0.4,0.02,'Bottom Concentration Anomaly ($\mu$mol/L)',transform=ax.transAxes)
    SB = ax.contour(grid.X[xslice]/1000,grid.Y[yslice]/1000,
                     grid.Depth[yslice,xslice],
                    [150.0],
                    colors=['0.1'],
                    linewidths=[0.75] )
    
    ax.axhline(grid.Y[yind]/1000,color='k',linestyle='--')
    ax.axvline(grid.X[xind]/1000,color='k',linestyle='--')
    
    ax.set_xlabel('Alongshore distance (km)',labelpad=0.3)
    ax.set_ylabel('Cross-shore distance (km)',labelpad=0.3)
# Grid, state and tracers datasets of base case
grid_file = '/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run38/gridGlob.nc'
grid = xr.open_dataset(grid_file)

grid_fileNoC = '/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run42/gridGlob.nc'
gridNoC = xr.open_dataset(grid_fileNoC)

state_file = '/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run38/stateGlob.nc' 
state = xr.open_dataset(state_file)

ptracers_file = '/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run38/ptracersGlob.nc'
ptracers = xr.open_dataset(ptracers_file)

ptracersNoC_file = '/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run42/ptracersGlob.nc'
ptracersNoC = xr.open_dataset(ptracersNoC_file)

#RhoRef = np.squeeze(rdmds('/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run38/RhoRef'))
RhoRef = 999.79998779 # It is constant in all my runs, can't run rdmds
# General input
nx = 360
ny = 360
nz = 90
nt = 19 # t dimension size 

xslice=slice(100,360)
yslice=slice(100,300)

tslice = slice(6,10)

xind = 240
yind = 267 # y index for alongshore cross-section
zind = 27

zslice = slice(0,22)
# plot2
zslice_2 = slice(0,32)
# plot 3
yslice_bac = slice(200,300)
xslice_bac = slice(100,360)

hFacmasked = np.ma.masked_values(grid.HFacC.data, 0)
MaskC = np.ma.getmask(hFacmasked)

hFacmaskedNoC = np.ma.masked_values(gridNoC.HFacC.data, 0)
MaskCNoC = np.ma.getmask(hFacmaskedNoC)

MaskExpand = np.expand_dims(MaskC,0) 
MaskExpand = MaskExpand + np.zeros((ptracers.Tr1).shape)    
ConcAreaWeighted, ConcFilt, Area, BottInd = ConcAreaShelf(ptracers.Tr1.data, grid.HFacC.data, grid.rA.data,
                                                 grid.Depth.data)

ConcAreaWeightedNoC, ConcFiltNoC, AreaNoC,BottIndNoC = ConcAreaShelf(ptracersNoC.Tr1.data,gridNoC.HFacC.data, 
                                                                     gridNoC.rA.data, gridNoC.Depth.data)
cmap_name = cmo.cm.tempo
sns.set_context('paper')
sns.set_style('white')

plt.rcParams['font.size'] = 8.0

f = plt.figure(figsize = (7.4,4.6)) # 19.0cm = 7.4in, 115cm = 4.52in

gs0 = gspec.GridSpec(2, 2, height_ratios=[0.6,1], width_ratios=[1,0.6])

Plot1(gs0[0,0],ptracers.Tr1.isel(Y=267),zslice, yind, xslice) # x-z contours conc

Plot2(gs0[1,1],ptracers.Tr1.isel(X=250),yslice_bac,zslice_2,250) # bottom concentration

Plot3(gs0[1,0],ConcFilt,tslice,xslice,yslice_bac,250,267)



plt.tight_layout()

plt.savefig('fig_BottomConcBaseNew.pdf',format='pdf')

def next_frame(t):
    
    timeIndex = starttime + t 
    data = ConcFilt[timeIndex,yslice_bac,xslice]
    gs0 = gspec.GridSpec(2, 2, height_ratios=[0.6,1], width_ratios=[1,0.6])

    Plot1(gs0[0,0],ptracers.Tr1.isel(Y=267),zslice, yind, xslice) # x-z contours conc

    Plot2(gs0[1,1],ptracers.Tr1.isel(X=250),yslice_bac,zslice_2,250) # bottom concentration

    Plot3(gs0[1,0],ConcFilt,tslice,xslice,yslice_bac,250,267)

    plt.tight_layout()

# frames for beamer animation

cmap_name = cmo.cm.tempo
sns.set_context('paper')
sns.set_style('white')

plt.rcParams['font.size'] = 8.0

f = plt.figure(figsize = (7.4,4.6)) # 19.0cm = 7.4in, 115cm = 4.52in

starttime = 0
data = ConcFilt[starttime,yslice_bac,xslice]
data_all_times = ConcFilt[:,yslice_bac,xslice]

maxCM = 10
minCM = 4.4

minXy = np.min(np.nanmean(data_all_times,axis=2)) 
maxXy = np.max(np.nanmean(data_all_times,axis=2))
minYx = np.min(np.nanmean(data_all_times,axis=1))
maxYx = np.max(np.nanmean(data_all_times,axis=1))

fig = plt.figure()
    
for tt in range(19):
    
    next_frame(tt)
    frame_name2 = ('/ocean/kramosmu/Figures/PoolHighConcentration/BottomConc_basecase_tt%02d.png' %tt) 
    plt.savefig(frame_name2,  format='png')


# Find mask shelf bottom

hold_maskSB = np.empty((ny,nx))

for j in range (nx):
    for i in range(ny):
        if  abs(BottInd[i,j]-BottIndNoC[i,j]) <= 1:
            hold_maskSB[i,j] = 1
        else:
            hold_maskSB[i,j]= 0
    
maskSB = np.ma.masked_equal(hold_maskSB,0)     
bathyMasked = np.ma.masked_less(-gridNoC.Depth, -152.5)

MASK_SB = np.ma.mask_or(bathyMasked.mask, maskSB.mask)

maskExpandSB = np.expand_dims(MASK_SB,0) 
maskExpandSB = maskExpandSB + np.zeros(ConcFilt.shape)   


cmap_name = cmo.cm.balance
sns.set_context('paper')
sns.set_style('white')

plt.rcParams['font.size'] = 8.0

f = plt.figure(figsize = (7.4,4.4)) # 19.0cm = 7.4in, 115cm = 4.52in

gs0 = gspec.GridSpec(2, 2, height_ratios=[0.6,1], width_ratios=[1,0.6])


trmax = 5
trmin= -5

Plot1anom(gs0[0,0],ptracers.Tr1.isel(Y=yind),ptracersNoC.Tr1.isel(Y=yind),
          tslice,zslice,267,xslice,Trmin=trmin,Trmax=trmax) # x-z contours conc

Plot2anom(gs0[1,1],ptracers.Tr1.isel(X=xind),ptracersNoC.Tr1.isel(X=xind),
          tslice,yslice_bac,zslice_2,250,Trmin=trmin,Trmax=trmax) # bottom concentration

Plot3anom(gs0[1,0],np.ma.masked_array(ConcFilt-ConcFiltNoC,mask=maskExpandSB),
          tslice,xslice,yslice_bac,250,267,trmin,trmax)



plt.tight_layout()

plt.savefig('fig_BottomConcAnomBaseNew.eps',format='eps', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16158234126984128, 'fraction_numerical': 0.050285218253968256, 'mean_word_length': 2.7301110083256246, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '646954', 'license': 'Apache-2.0', 'n_tokens_mistral': 6990, 'n_tokens_neox': 6426, 'n_words': 924, 'path': 'PythonScripts/Paper1Figures/Fig_BottomConcBaseCase.ipynb', 'repo_name': 'UBC-MOAD/outputanalysisnotebooks', 'size': '16128'}