<jupyter_start><jupyter_text>Mapping tracer concentration on shelf bottom Map hcw on shelf by looking at bottom cells and multiplying their concentration by their area.The algorithm: Using hfacC finfd the first occurrence of an open cell, grab the x,y,z index and then use  those indices to fing the corresponding tracer concentration.<jupyter_code>#KRM

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.colors as mcolors

import matplotlib.cm as cm

import matplotlib as mpl

from math import *

from mpl_toolkits.mplot3d.axes3d import Axes3D

import scipy.io

import scipy as spy

%matplotlib inline

from netCDF4 import Dataset

import pylab as pl

import os 

import sys

import seaborn as sns

lib_path = os.path.abspath('../../Building_canyon/BuildCanyon/PythonModulesMITgcm') # Add absolute path to my python scripts
sys.path.append(lib_path)

import ReadOutTools_MITgcm as rout 

sns.set()
sns.set_style('white')
sns.set_context('talk')
#Base case, iso =1 , No 3d diff.
#CanyonGrid='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/gridGlob.nc'
CanyonGrid='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/gridGlob.nc'
CanyonGridOut = Dataset(CanyonGrid)

NoCCanyonGrid='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run07/gridGlob.nc'
NoCCanyonGridOut = Dataset(NoCCanyonGrid)

CanyonState='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/stateGlob.nc'
CanyonStateOut = Dataset(CanyonState)

Ptracer='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/ptracersGlob.nc'
PtracerOut = Dataset(Ptracer)

Ptracer12='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run12/ptracersGlob.nc'
PtracerOut12 = Dataset(Ptracer12)

Ptracer07='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run07/ptracersGlob.nc'
PtracerOut07 = Dataset(Ptracer07)

Ptracer25='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run25/ptracersGlob.nc'
PtracerOut25 = Dataset(Ptracer25)


#PtracerNoC='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run07/ptracersGlob.nc'
#PtracerNoCOut = Dataset(PtracerNoC)

nx = 360
ny = 360
nz = 90
nt = 19 # t dimension size 

xc = rout.getField(CanyonGrid, 'XC') # x coords tracer cells
yc = rout.getField(CanyonGrid, 'YC') # y coords tracer cells
rc = CanyonGridOut.variables['RC']
bathy = rout.getField(CanyonGrid, 'Depth')
bathyNoC = rout.getField(NoCCanyonGrid, 'Depth')

hFacC = rout.getField(CanyonGrid, 'HFacC')
MaskC = rout.getMask(CanyonGrid, 'HFacC')

hFacCNoC = rout.getField(NoCCanyonGrid, 'HFacC')
MaskCNoC = rout.getMask(NoCCanyonGrid, 'HFacC')
rANoC = rout.getField(NoCCanyonGrid, 'rA')
rA = rout.getField(CanyonGrid, 'rA')

z = CanyonStateOut.variables['Z']

time = CanyonStateOut.variables['T']
Tr03 = rout.getField(Ptracer,'Tr1') # 

Tr12 = rout.getField(Ptracer12,'Tr3') # 

Tr07 = rout.getField(Ptracer07,'Tr1') # 
Tr25 = rout.getField(Ptracer25,'Tr1') # 
def ConcArea(Tr, hfac, ra):
    '''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)
       RETURNS:
       ConcArea = concentration at cell closest to bottom times its area (nt,ny,nx)
       Conc = cocnetration near bottom (nt,ny,nx)'''
    
    ConcArea = np.empty((19,360,360))
    Conc = np.empty((19,360,360))
    Area = np.empty((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):
        #print(tt)
        for i in range(360):
            for j in range(360):
                TrBottom = Tr[tt,BottomInd[i,j],i,j]
                ConcArea[tt,i,j] = TrBottom*ra[i,j]
                Conc[tt,i,j] = TrBottom
                Area[i,j] = ra[i,j]
    
    print(np.shape(ConcArea))

    #ConcAreaMasked = np.ma.masked_values(ConcDepths,-2.5)
    #ConcDepths[np.where(np.ma.getmask(ConcDepthsMasked)==True)] = np.nan
    return (ConcArea,Conc, Area)
sns.set_context('talk',font_scale=1.0)
#TrBottom = Tr1[0,BottomInd[i,j],50,50]
#print(Tr1[0,BottomInd[i,j],i,j])
plt.pcolor(xc[1,:],yc[:,1],rA[:,:],cmap='spectral')
cbar = plt.colorbar()
plt.title('rA')
plt.ylabel('Cross-shore distance (m)')
plt.xlabel('Alongshore distance (m)')
ConcArea1,Conc1,Area1=ConcArea(Tr03, hFacCNoC, rANoC)
sns.set_context('talk',font_scale=1.0)
plt.pcolor(xc[1,:],yc[:,1],Conc1[10,:,:],cmap='spectral')
cbar = plt.colorbar()
plt.title('Concentration at bottom, gray means land')
plt.ylabel('Cross-shore distance (m)')
plt.xlabel('Alongshore distance (m)')
ConcArea3,Conc3,Area3=ConcArea(Tr12, hFacCNoC, rANoC)
sns.set_context('talk', font_scale=1.0)

for tt in np.arange(0,19,2):
    
    fig,((ax1),(ax2)) = plt.subplots(1, 2, figsize=(20, 4),  sharey=True, sharex=True)
    
    iso = ax1.contourf(xc[227:,:],yc[227:,:],Conc1[tt,227:,:]-Conc1[0,227:,:],20,
                       vmin=-2.8,vmax=2.8,cmap='RdYlBu_r')
    cb = plt.colorbar(iso, shrink=1,ax=ax1)
    #cb.set_label('')
    cset = ax1.contour(xc[227:,:],yc[227:,:],-bathy[227:,:],[-300,-200, rc[30],-125,-100,-75,-50],
                       colors=['0.50','0.50','k','0.50','0.50','0.50','0.50'],linewidth=1, rstride=2, cstride=2,zdir='z')

    ax1.set_aspect('auto')
    ax1.set_title('Tr1 Shelf bottom C-C(tt=0) day %0.1f ' %(tt/2.0))
    ax1.set_ylabel('Across-shore distance (m)')
    ax1.set_xlabel('Alongshore distance (m)')
    
    iso2 = ax2.contourf(xc[227:,:],yc[227:,:],Conc3[tt,227:,:]-Conc3[0,227:,:],20,
                        vmin=-2.8,vmax=2.8,cmap='RdYlBu_r')
    cb2 = plt.colorbar(iso2, shrink=1,ax=ax2)
    #cb2.set_label('m')
    cset = ax2.contour(xc[227:,:],yc[227:,:],-bathy[227:,:],[-300,-200, rc[30],-125,-100,-75,-50],
                       colors=['0.50','0.50','k','0.50','0.50','0.50','0.50'],linewidth=1, rstride=2, cstride=2,zdir='z')

    ax2.set_aspect('auto')
    ax2.set_title('Tr3 Shelf bottom C-C(tt=0) day %0.1f ' %(tt/2.0))
    #ax2.set_ylabel('Across-shore distance (m)')
    ax2.set_xlabel('Alongshore distance (m)')
    
    plt.tight_layout
    plt.show

#fig.savefig('/ocean/kramosmu/Figures/BathyShelfBoxComplete.eps', format='eps', dpi=1000, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Note that I used the flat shelf mask to find the "shelf bottom". This means that the region over the canyon is the concentration*area over the lid of the canyon.<jupyter_code>ConcArea1noc,Conc1noc,Area1noc=ConcArea(Tr07, hFacCNoC, rANoC)
ConcArea3noc,Conc3noc,Area3noc=ConcArea(Tr25, hFacCNoC, rANoC)
import matplotlib.animation as animation
def Plot(tt,ax1,ax2,ax3):
    ax1.clear()
    iso = ax1.contourf(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,
                       Conc1[tt,227:300,100:],15,
                       vmax=2.2, vmin=-2.2,
                       cmap='Spectral')
    
    cset = ax1.contour(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,-bathy[227:300,100:],
                       [-300,-200, rc[30],-125,-100,-75,-50],
                       colors=['0.50','0.50','k','0.50','0.50','0.50','0.50'],
                       linewidth=1, 
                       rstride=2, 
                       cstride=2,
                       zdir='z')
    
    ax1.clabel(cset, inline=1, fontsize=10)
    ax1.set_aspect('auto')
    ax1.set_title('Base case')
    ax1.set_ylabel('Across-shore distance (km)')
    ax1.set_xlabel('Alongshore distance (km)')
    
    #-------------------------------------------------------------------------------------------------------------------
    ax2.clear()
    iso2 = ax2.contourf(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,
                        Conc1noc[tt,227:300,100:],15,
                        vmax=2.2, vmin=-2.2,
                        cmap='Spectral')
    
    cax,kw = mpl.colorbar.make_axes([ax1,ax2],location='bottom',shrink=1, pad = -0.5)
    cb = plt.colorbar(iso2, cax=cax,**kw)
    cb.set_label('Mol/l')
    
    cset = ax2.contour(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,-bathyNoC[227:300,100:],
                       [-300,-200, rc[30],-125,-100,-75,-50], 
                       colors=['0.50','0.50','k','0.50','0.50','0.50','0.50'],
                       linewidth=1, 
                       rstride=2, 
                       cstride=2,
                       zdir='z')
    ax2.clabel(cset, inline=1, fontsize=10)
    ax2.set_aspect('auto')
    ax2.set_title('Locally-enhanced $k_v$')
    ax2.set_xlabel('Alongshore distance (km)')
    
    #-------------------------------------------------------------------------------------------------------------------
    ax3.clear()
    iso3 = ax3.contourf(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,
                        Conc1[tt,227:300,100:]-Conc1noc[tt,227:300,100:],15,
                        vmax=2.2, vmin=-2.2,
                        cmap='PiYG')
    cax3,kw = mpl.colorbar.make_axes([ax3],location='bottom',shrink=1, pad = -0.5)
    cb3 = plt.colorbar(iso3, cax=cax3,**kw)
    cb3.set_label('Mol/l')
    
    cset = ax3.contour(xc[227:300,100:]/1000.0,yc[227:300,100:]/1000.0,-bathy[227:300,100:],
                       [-300,-200, rc[30],-125,-100,-75,-50],
                       colors=['0.50','0.50','k','0.50','0.50','0.50','0.50'],
                       linewidth=1, 
                       rstride=2, 
                       cstride=2,
                       zdir='z')
    
    ax3.clabel(cset, inline=1, fontsize=10)
    ax3.set_aspect('auto')
    ax3.set_title('High $k_v$ - day %0.1f ' %(tt/2.0))
    ax3.set_xlabel('Alongshore distance (km)')
    
    #plt.tight_layout()
    
## Animation
#N=5


#Empty figures
sns.set_context('talk', font_scale=1.0)
fig,((ax1),(ax2),(ax3)) = plt.subplots(1, 3, figsize=((12, 7)),  sharey=True, sharex=True)

#Initial image
def init():
    Plot(0,ax1,ax2,ax3)
    
    #return[ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8,ax9]
    plt.tight_layout()

def animate(tt):
    Plot(tt,ax1,ax2,ax3)
    
#The animation function (max frames=47)   


anim = animation.FuncAnimation(fig, animate, init_func=init,frames=18, interval = 200,blit=False, repeat=False)

##A line that makes it all work
mywriter = animation.FFMpegWriter()
##Save in current folder
anim.save('TracerShelfBottom_BaseCase.mp4',writer=mywriter,fps=0.2, dpi = 300, bitrate = 1000000, codec = "libx264")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18934854049874347, 'fraction_numerical': 0.07819447129325344, 'mean_word_length': 3.0544670846394983, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '235877', 'license': 'Apache-2.0', 'n_tokens_mistral': 4606, 'n_tokens_neox': 4085, 'n_words': 611, 'path': 'ConcentrationOnShelfBottom/ConcentrationOnShelfBottom.ipynb', 'repo_name': 'UBC-MOAD/outputanalysisnotebooks', 'size': '10346'}