<jupyter_start><jupyter_text>Compare latitudinal grid spacing in GISS-H, GISS-G, and ECCO<jupyter_code>import os
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import pych.ecco as pe
import ecco_v4_py as ecco
from pych.plot import nice_inward_ticks
plt.rcParams.update({'font.size':18,'axes.labelsize':18,'axes.titlesize':22,'text.usetex':True})
fig_dir = '/workspace/results/eccov4r4/equatorial-under-current/figures/'
if not os.path.isdir(fig_dir):
    os.makedirs(fig_dir)<jupyter_output><empty_output><jupyter_text>ECCOv4<jupyter_code>ds = pe.read_mds('/workspace/grids/llc90',iters=None)
grid = ecco.get_llc_grid(ds)<jupyter_output><empty_output><jupyter_text>Meridional Grid Spacing<jupyter_code>lon = -140
maskW,maskS = ecco.calc_section_trsp._parse_section_trsp_inputs(ds,[lon,2],[lon,-2],None,None,'EUC')
lon=-140
maskW,maskS = ecco.calc_section_trsp._parse_section_trsp_inputs(ds,[lon,2],[lon,-2],None,None,'EUC')
ecco_center_points = grid.interp(ds.YC,'Y',boundary='fill').where(maskS,drop=True).squeeze()
ecco_bds = ecco_center_points - ecco_center_points.diff('i')/2 
ecco_dLat = (ds.dxG.where(maskS,drop=True)/1000).squeeze()
ecco_dLatMid = (ecco_dLat+ecco_dLat.shift(i=-1)).dropna('i')/2<jupyter_output><empty_output><jupyter_text>GISS-G<jupyter_code>gissg = xr.open_dataset('/workspace/results/giss-euc/giss_g.nc')
gissg_center_points = gissg.YC.sel(YC=slice(-2.5,2.5))
gissg_bds = gissg_center_points - gissg_center_points.diff('YC')/2 

gissg_dLat = (gissg.dyG.sel(YC=slice(-2.5,2.5)).sel(XG=lon,method='nearest')/1000)
gissg_dLatMid = (gissg_dLat+gissg_dLat.shift(YC=-1)).dropna('YC')/2<jupyter_output><empty_output><jupyter_text>GISS-H<jupyter_code>gissh = xr.open_dataset('/workspace/results/giss-euc/giss_h.nc')
gissh_center_points = gissh.YC.sel(YC=slice(-2.5,2.5))
gissh_bds = gissh_center_points - gissh_center_points.diff('YC')/2 

gissh_dLat = (gissh.dyG.sel(YC=slice(-2.5,2.5)).sel(XG=lon,method='nearest')/1000)
gissh_dLatMid = (gissh_dLat+gissh_dLat.shift(YC=-1)).dropna('YC')/2
gissh.YG<jupyter_output><empty_output><jupyter_text>Plots  Grid spacing<jupyter_code>x=ecco_bds.where(np.abs(ecco_bds)<.45).values
ecco_dLatMid.where(~np.isnan(x)).values
def get_highlight(x,y,cutoff,ax,xextra=None,yextra=None,**kwargs):
    index = np.abs(x)<cutoff
    xx = x.where(index).values
    yy = y.where(~np.isnan(xx)).values
    if xextra is not None:
        xxx = np.concatenate([[-xextra],xx,[xextra]])
        yyy = np.concatenate([[yextra],yy,[yextra]])
    else:
        xxx=xx
        yyy=yy
    xxx=xxx[~np.isnan(xxx)]
    index = np.argsort(xxx)
    xxx = xxx[index]
    yyy=yyy[~np.isnan(yyy)][index]
    ax.plot(xxx,yyy,**kwargs)
nrows= 3
fig,axs = plt.subplots(nrows,1,figsize=(12,nrows*4))

ax=axs[0]
ax.plot(ecco_center_points,ecco_dLat,marker='.',markersize=10,color='k')
ax.plot(ecco_bds,ecco_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color='k')
get_highlight(ecco_center_points,ecco_dLat,cutoff=.3,ax=ax,
             marker='.',markersize=10,color='r',linestyle='-')
get_highlight(ecco_bds,ecco_dLatMid,cutoff=.45,ax=ax,
             marker='|',markersize=30,markeredgewidth=3,linestyle='',color='r')
get_highlight(ecco_center_points,ecco_dLat,cutoff=.3,ax=ax,
              xextra=ecco_bds[3],yextra=ecco_dLatMid[3],
             marker='',markersize=10,color='r',linestyle='-')
ax.set(ylim=[44,54],xlim=[-2.1,2.1])
ax.text(-0.2,52,'ECCOv4')

ax = axs[1]
g = ax.plot(gissg_center_points,gissg_dLat,marker='.',markersize=10,color='k')
ax.plot(gissg_bds,gissg_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color=g[0].get_color())
get_highlight(gissg_center_points,gissg_dLat,cutoff=.6,ax=ax,
             marker='.',markersize=10,color='r',linestyle='-')
get_highlight(gissg_bds,gissg_dLatMid,cutoff=1.1,ax=ax,
             marker='|',markersize=30,markeredgewidth=3,linestyle='',color='r')
get_highlight(gissg_center_points,gissg_dLat,cutoff=.6,ax=ax,
              xextra=gissg_bds[3],yextra=gissg_dLatMid[3],
             marker='',markersize=10,color='r',linestyle='-')
ax.set(ylim=[108,114],xlim=[-2.1,2.1])
ax.text(-0.15,113,'E2.1-G')#,bbox={'facecolor':'white','boxstyle':'round,pad=.5'})

ax = axs[2]
h = ax.plot(gissh_center_points,gissh_dLat,marker='.',markersize=10,color='k')
ax.plot(gissh_bds,gissh_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color=h[0].get_color())
get_highlight(gissh_center_points,gissh_dLat,cutoff=.6,ax=ax,
             marker='.',markersize=10,color='r',linestyle='-')
get_highlight(gissh_bds,gissh_dLatMid,cutoff=1.1,ax=ax,
             marker='|',markersize=30,markeredgewidth=3,linestyle='',color='r')
get_highlight(gissh_center_points,gissh_dLat,cutoff=.6,ax=ax,
              xextra=gissh_bds[3],yextra=gissh_dLatMid[3],
             marker='',markersize=10,color='r',linestyle='-')
ax.set(ylim=[108,114],xlim=[-2.1,2.1])
ax.text(-0.15,113,'E2.1-H')#,bbox={'facecolor':'white','boxstyle':'round,pad=.5'})


fig.suptitle(f'Grid cell width in meridional direction near the equator\n'+
              'Dots indicate the grid cell center, Vertical lines indicate grid cell edges',y=.95)
[ax.set_ylabel('Grid cell width, km') for ax in axs]
[ax.set_xticklabels([r'%.1f$^\circ$%s' % (x,y) for x,y in \
                    zip([2,2,1.5,1,.5,0,.5,1,1.5,2],['S']*5+['']+['N']*4)]) for ax in axs]
[nice_inward_ticks(ax,xminor_skip=.1,yminor_skip=1) for ax in axs];
axs[2].set_xlabel('Latitude')
fig.savefig(f'{fig_dir}/compare_grid_cell_width_equator.png',
            bbox_inches='tight',dpi=300,pad_inches=1)<jupyter_output><empty_output><jupyter_text>Vertical Layer Thickness<jupyter_code>def autolabel(rects, labels=None, ypos='center'):
    """
    Attach a text label above each bar in *rects*, displaying its height.

    *ypos* indicates which side to place the text w.r.t. the center of
    the bar. It can be one of the following {'center', 'right', 'left'}.
    """

    ha = {'center': 'center', 'right': 'left', 'left': 'right'}
    offset = {'center': 0, 'right': 1, 'left': -1}

    for i,rect in enumerate(rects):
        width = rect.get_width()
        label = labels[i] if labels is not None else '{}'.format(width)
        ax.annotate(label,
                    xy=(width,rect.get_y() - rect.get_height() / 3),
                    xytext=(40,offset[ypos]*4-12),  # use 3 points offset
                    textcoords="offset points",  # in both directions
                    ha=ha[ypos], va='bottom',
                   fontsize=14)
        
    return ax
max_depth=450
ncols = 3
fig,axs = plt.subplots(1,ncols,figsize=(15,10))

labels = ('ECCOv4r4','E2.1-G','E2.1-H')
for d,ax,lbl in zip([ds,gissg,gissh],axs,labels):
    if d.Z[0]<0:
        z=-d.Z
    else:
        z= d.Z
    if 'k' in d.dims:
        ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in zip(d.k.where(z<max_depth,drop=True).values,d.drF.where(z<max_depth,drop=True).values)]
    else:
        ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in enumerate(d.drF.where(z<max_depth,drop=True).values)]
        
    bars = ax.barh(z.where(z<max_depth,drop=True),d.drF.where(z<max_depth,drop=True),height=3);
    autolabel(bars,labels=ll)

    ax.invert_yaxis()
    ax.set(ylim=[450,0],xlim=[0,110],ylabel='Depth, m',
           xlabel=f'Vertical Grid Cell Thickness, m');
    
    ax.text(70,25,lbl)#bbox={'facecolor':'white','boxstyle':'round,pad=.5'})
    nice_inward_ticks(ax,xminor_skip=5,yminor_skip=10)
    ax.axhline(y=400,xmin=0,xmax=110,color='gray',alpha=.8,linestyle='--')

[ax.set_ylabel('') for ax in axs[1:]];

fig.suptitle('Vertical Grid Cell Thickness in Upper Water Column\nBlue bar is located at grid cell center',y=.98)
fig.savefig(f'{fig_dir}/compare_grid_layer_thickness.png',
            bbox_inches='tight',dpi=300,pad_inches=1)<jupyter_output><empty_output><jupyter_text>Now compare to ADCP<jupyter_code>adcp = xr.open_dataset('/mnt/data/datasets/tao_triton/adcp_xyzt_dy.cdf')
adcp = adcp.sel(lon=220,lat=0)
adcp = adcp.rename({'U_1205':'uvel','V_1206':'vvel'})
def replace_nans(ds,fld):
    ds[fld] = xr.where(ds[fld]==ds.missing_value,np.NAN,ds[fld])
    return ds[fld]
for f in adcp.data_vars:
    adcp[f] = replace_nans(adcp,f)
adcp.uvel.sel(time=slice('2007','2007')).plot()
adcp.depth
drf = np.array([5,7.5]+[5]*(len(adcp.depth)-2))
adcp['drF'] = xr.DataArray(drf,coords=adcp.depth.coords,dims=adcp.depth.dims)
ncols = 2
fig,axs = plt.subplots(1,ncols,figsize=(15,10))

labels = ('ECCOv4r4','GTMBA ADCP\n'+r'140$^\circ$W, 0$^\circ$N'+'\nMeasures every 5m')
for d,ax,lbl,xmax in zip([ds,adcp],axs,labels,[110,10]):
    if 'Z' in d.coords:
        if d.Z[0]<0:
            z=-d.Z
        else:
            z= d.Z
        if 'k' in d.dims:
            ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in zip(d.k.where(z<max_depth,drop=True).values,d.drF.where(z<max_depth,drop=True).values)]
        else:
            ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in enumerate(d.drF.where(z<max_depth,drop=True).values)]
    else:
        z=d.depth
        ll = ['']*len(z)
        
    bars = ax.barh(z.where(z<max_depth,drop=True),d.drF.where(z<max_depth,drop=True),height=3);
    autolabel(bars,labels=ll)

    ax.invert_yaxis()
    ax.set(ylim=[450,0],xlim=[0,xmax],ylabel='Depth, m',
           xlabel=f'Vertical Grid Cell Thickness, m');
    
    if xmax == 110:
        ax.text(70,70,lbl)#bbox={'facecolor':'white','boxstyle':'round,pad=.5'})
    elif xmax == 10:
        ax.text(6,70,lbl)#bbox={'facecolor':'white','boxstyle':'round,pad=.5'})
    else:
        raise NotImplementedError('xmax is what?')
    nice_inward_ticks(ax,xminor_skip=5,yminor_skip=10)
    ax.axhline(y=400,xmin=0,xmax=xmax,color='gray',alpha=.8,linestyle='--')

[ax.set_ylabel('') for ax in axs[1:]];

fig.suptitle('Vertical Grid Cell Thickness in Upper Water Column\nBlue bar is located at grid cell center',y=.98)
fig.savefig(f'{fig_dir}/compare_ecco_adcp_vertical.png',
            bbox_inches='tight',dpi=300,pad_inches=1)
max_depth = 90
min_depth = 60
np.array(ll)[]
ind =adcp.depth.where((z<max_depth),drop=True)
ind
np.array(ll)
ncols = 2
fig,axs = plt.subplots(1,ncols,figsize=(15,10))

labels = ('ECCOv4r4','GTMBA ADCP\n'+r'140$^\circ$W, 0$^\circ$N')
for d,ax,lbl,xmax in zip([ds,adcp],axs,labels,[15,10]):
    if 'Z' in d.coords:
        if d.Z[0]<0:
            z=-d.Z
        else:
            z= d.Z
        if 'k' in d.dims:
            ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in zip(d.k.where(z<max_depth,drop=True).values,d.drF.where(z<max_depth,drop=True).values)]
        else:
            ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in enumerate(d.drF.where(z<max_depth,drop=True).values)]
    else:
        z=d.depth
        ll = [f'{int(x+1)}: {Z:.1f} m' for x,Z in enumerate(d.drF.where(z<max_depth,drop=True).values)]
        
    ind = ~np.isnan(d.drF.where(z<max_depth,drop=True).where(z>min_depth).values)
    ll = list(np.array(ll)[ind])
    bars = ax.barh(z.where(z<max_depth,drop=True)[ind],d.drF.where(z<max_depth,drop=True)[ind],height=3);
    ax = autolabel(bars,labels=ll)

    ax.invert_yaxis()
    ax.set(ylim=[90,60],xlim=[0,xmax],ylabel='Depth, m',
           xlabel=f'Vertical Grid Cell Thickness, m');
    
    if xmax == 110:
        ax.text(70,71,lbl)
    elif xmax == 10:
        ax.text(1,63,lbl)
    elif xmax == 15:
        ax.text(1,62.5,lbl)
    else:
        raise NotImplementedError('xmax is what?')
    nice_inward_ticks(ax,xminor_skip=5,yminor_skip=10)

[ax.set_ylabel('') for ax in axs[1:]];

fig.suptitle('Vertical Grid Cell Thickness in Upper Water Column\nBlue bar is located at grid cell center',y=.98)
fig.savefig(f'{fig_dir}/compare_ecco_adcp_vertical_zoom.png',
            bbox_inches='tight',dpi=300,pad_inches=1)<jupyter_output><empty_output><jupyter_text>ECCO v4r4<jupyter_code>plt.rcParams.update({'text.usetex':False})
ds.XC.sel(tile=8)[:,10].values
ds.YC.sel(tile=8)[0,:].values
ds.dxC.sel(tile=7).plot()
ds.dxC.dims
ecco_center_points = ds.XC.sel(tile=8)
ecco_dLon = ds.dxC.sel(tile=8).sel()
nrows= 3
fig,axs = plt.subplots(nrows,1,figsize=(12,nrows*4))

ax=axs[0]
ax.plot(ecco_center_points,ecco_dLat,marker='.',markersize=10,color='k')
ax.plot(ecco_bds,ecco_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color='k')
ax.set(ylim=[44,54],xlim=[-2.1,2.1])
ax.text(-0.2,52,'ECCOv4')

ax = axs[1]
g = ax.plot(gissg_center_points,gissg_dLat,marker='.',markersize=10,color='k')
ax.plot(gissg_bds,gissg_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color=g[0].get_color())
ax.set(ylim=[108,114],xlim=[-2.1,2.1])
ax.text(-0.15,113,'E2.1-G')#,bbox={'facecolor':'white','boxstyle':'round,pad=.5'})

ax = axs[2]
h = ax.plot(gissh_center_points,gissh_dLat,marker='.',markersize=10,color='k')
ax.plot(gissh_bds,gissh_dLatMid,marker='|',markersize=30,markeredgewidth=3,linestyle='',color=h[0].get_color())
ax.set(ylim=[108,114],xlim=[-2.1,2.1])
ax.text(-0.15,113,'E2.1-H')#,bbox={'facecolor':'white','boxstyle':'round,pad=.5'})


fig.suptitle(f'Grid cell width in meridional direction near the equator\n'+
              'Dots indicate the grid cell center, Vertical lines indicate grid cell edges',y=.95)
[ax.set_ylabel('Grid cell width, km') for ax in axs]
[ax.set_xticklabels([r'%.1f$^\circ$%s' % (x,y) for x,y in \
                    zip([2,2,1.5,1,.5,0,.5,1,1.5,2],['S']*5+['']+['N']*4)]) for ax in axs]
[nice_inward_ticks(ax,xminor_skip=.1,yminor_skip=1) for ax in axs];
axs[2].set_xlabel('Latitude')
fig.savefig(f'{fig_dir}/compare_grid_cell_width_equator.png',
            bbox_inches='tight',dpi=300,pad_inches=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1970916069978593, 'fraction_numerical': 0.03845869934302798, 'mean_word_length': 4.885317115551694, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '108334', 'license': 'MIT', 'n_tokens_mistral': 6238, 'n_tokens_neox': 5902, 'n_words': 686, 'path': 'notebooks/eccov4_giss_compare_equatorial_grid.ipynb', 'repo_name': 'timothyas/eccov4-euc', 'size': '13547'}