<jupyter_start><jupyter_text>Below just left overs from development<jupyter_code>#include masks psi u v in smoothing and deepening
x_earth = np.empty((rows,columns))
y_earth = np.empty((rows,columns))

x_earth[:,x_0ind] = 0
y_earth[y_0ind,:] = 0

for column in np.arange(columns-x_0ind-1):
    dx = haversine(lon[:,x_0ind+column],lat[:,x_0ind+column],lon[:,x_0ind+column+1],lat[:,x_0ind+column+1])
    x_earth[:,x_0ind+column+1] = x_earth[:,x_0ind+column]+dx
    
for column in np.arange(x_0ind):    
    dx = haversine(lon[:,x_0ind-column],lat[:,x_0ind-column],lon[:,x_0ind-column-1],lat[:,x_0ind-column-1])
    x_earth[:,x_0ind-column-1] = x_earth[:,x_0ind-column]-dist
    
for row in np.arange(rows-y_0ind-1):  
    dy = haversine(lon[y_0ind+row,:],lat[y_0ind+row,:],lon[y_0ind+row+1,:],lat[y_0ind+row+1,:])
    y_earth[y_0ind+row+1,:] = y_earth[y_0ind+row,:]+dy

for row in np.arange(y_0ind):
    dy = haversine(lon[y_0ind-row,:],lat[y_0ind-row,:],lon[y_0ind-row-1,:],lat[y_0ind-row-1,:])
    y_earth[y_0ind-row-1,:] = y_earth[y_0ind-row,:]-dy
    
def calc_angle(lon_u,lat_u,lon_v,lat_v):
    
    DTOR = np.pi/180
    
    a1 = lat_u[1:,:] - lat_u[:-1,:]
    a2 = lon_u[1:,:] - lon_u[:-1,:]
    a2[a2<-180]+=360
    a2[a2>180]-=360
    
    a2 = a2*np.cos(0.5*DTOR*(lat_u[1:,:]+lat_u[:-1,:]))
    angle = np.arctan2(a1,a2)
    
    a1 = lat_v[:,:-1] - lat_v[:,1:]
    a2 = lon_v[:,:-1] - lon_v[:,1:]
    a2[a2<-180]+=360
    a2[a2>180]-=360
    
    a2 = a2*np.cos(0.5*DTOR*(lat_v[:,:-1] + lat_v[:,1:]))
    angle = 0.5*(angle + np.arctan2(a2,-a1))    
    
    return angle

angle = calc_angle(lon_u,lat_u,lon_v,lat_v)

print(np.shape(angle))
plt.close()
plt.contour(angle)
plt.colorbar()
plt.show()
x_rho_da = xr.DataArray(x_rho*1000,name='x_rho',dims=['eta_rho','xi_rho'],attrs={'long_name': 'X-location of RHO-points', 'units': 'meter'})
y_rho_da = xr.DataArray(y_rho*1000,name='y_rho',dims=['eta_rho','xi_rho'],attrs={'long_name': 'Y-location of RHO-points', 'units': 'meter'})
x_psi_da = xr.DataArray(x_psi*1000,name='x_psi',dims=['eta_rho','xi_rho'],attrs={'long_name': 'X-location of PSI-points', 'units': 'meter'})
y_psi_da = xr.DataArray(y_psi*1000,name='y_psi',dims=['eta_psi','xi_psi'],attrs={'long_name': 'Y-location of PSI-points', 'units': 'meter'})
x_u_da = xr.DataArray(x_u*1000,name='x_u',dims=['eta_u','xi_u'],attrs={'long_name': 'X-location of U-points', 'units': 'meter'})
y_u_da = xr.DataArray(y_u*1000,name='y_u',dims=['eta_u','xi_u'],attrs={'long_name': 'Y-location of U-points', 'units': 'meter'})
x_v_da = xr.DataArray(x_v*1000,name='x_v',dims=['eta_v','xi_v'],attrs={'long_name': 'X-location of V-points', 'units': 'meter'})
y_v_da = xr.DataArray(y_v*1000,name='y_v',dims=['eta_v','xi_v'],attrs={'long_name': 'Y-location of V-points', 'units': 'meter'})
%matplotlib notebook
#set up mask_rho depending on water column thickness
wct = (bed+ice).copy()
mask_rho = np.ones_like(wct)
mask_rho[wct<20] = 0

def make_neibhour_mask(mask):
    neib_xi = np.empty_like(mask)
    neib_eta = np.empty_like(mask)

    neib_xi[:,1:-1] = mask[:,2:]+mask[:,:-2]
    neib_xi[:,0] = neib_xi[:,1]
    neib_xi[:,-1] = neib_xi[:,-2]

    neib_eta[1:-1,:] = mask[2:,:]+mask[:-2,:]
    neib_eta[0,:] = neib_eta[1,:]
    neib_eta[-1,:] = neib_eta[-2,:]

    return neib_eta + neib_xi
#first get rid of any dubious ocean point nort of 60S
mask_sa = mask_rho.copy()
mask_sa[(lat_rho>-60)&(neib<=1)] = 0

#neib = make_neibhour_mask(mask_sa)
#mask_sa[(lat_rho>-60)&(neib<=3)] = 0

plt.close()
plt.pcolormesh(mask_sa)
plt.colorbar()
plt.show()
intdirintdirplt.close()
plt.pcolormesh(neib)
plt.colorbar()
plt.show()
#then, more carefully, find and mask ocean point around antarctica that are sorounded by land
Nmodif=1
while Nmodif > 0:
    mask_new = mask_tmp.copy()
    neib = make_neibhour_mask(mask_tmp)
    mask_new[(lat_rho<=-60)&(neib<=1)] = 0
    Nmodif = np.sum(mask_tmp)-np.sum(mask_new)
    print(Nmodif)
    mask_tmp = mask_new.copy()

mask_tmp[(lat_rho<=-60)&(neib<=2)] =0

Nmodif=1
while Nmodif > 0:
    mask_new = mask_tmp.copy()
    neib = make_neibhour_mask(mask_tmp)
    mask_new[(lat_rho<=-60)&(neib<=1)] = 0
    Nmodif = np.sum(mask_tmp)-np.sum(mask_new)
    print(Nmodif)
    mask_tmp = mask_new.copy()
plt.close()
plt.pcolormesh(mask_rho)
plt.colorbar()
plt.show()
%matplotlib notebook

#include the amundsen fix from Milan 2017

#first load and resample to whole roms rho grid
import scipy.io as sio
amundsen_fix_path = os.path.join(os.environ.get('extdir'),'rtopo','Bathymetry_ASE_Millan_et_al_2017.nc')
amu_fix = xr.open_dataset(amundsen_fix_path)
amundsen_fix_latlon_path = os.path.join(os.environ.get('extdir'),'rtopo','lon_lat_romain_grid.mat')
amu_fix_latlon = sio.loadmat(amundsen_fix_latlon_path)

amu_lon = amu_fix_latlon['lone']
amu_lat = amu_fix_latlon['late']
amu_lon[amu_lon>180]-=360.0
amu = resample(amu_lon,amu_lat,lon_rho,lat_rho,amu_fix.BED.values)

#then lookup indices of coorner coordinates and subset the working bathymetry just there 
from scipy.spatial import KDTree

points = np.column_stack((lon_rho.flatten(),lat_rho.flatten()))

tree = KDTree(points)

target = np.column_stack((amu_lon[0,0],amu_lat[0,0]))
llcDist, llcInd = tree.query(target)
llcInd = np.squeeze(llcInd)

target = np.column_stack((amu_lon[-1,-1],amu_lat[-1,-1]))
urcDist, urcInd = tree.query(target)
urcInd = np.squeeze(urcInd)

#establish grid coordinates at rho points only (needed for the amundsen fix)
xi_rho = np.arange((np.size(x,1)-1)/2,dtype=int)
eta_rho = np.arange((np.size(x,0)-1)/2,dtype=int)

xi_2d,eta_2d = np.meshgrid(xi_rho,eta_rho)
xiEta = np.column_stack((xi_2d.flatten(),eta_2d.flatten()))

xi_amu_ll,xi_amu_ur = xiEta[llcInd][0],xiEta[urcInd][0]
eta_amu_ll,eta_amu_ur = xiEta[llcInd][1],xiEta[urcInd][1]

bed[eta_amu_ll:eta_amu_ur+1,xi_amu_ll:xi_amu_ur+1] = amu[eta_amu_ll:eta_amu_ur+1,xi_amu_ll:xi_amu_ur+1]
bed=bed*-1
#have a look at the amundsen patch
plt.close()
plt.pcolormesh(bed[eta_amu_ll-20:eta_amu_ur+21,xi_amu_ll-20:xi_amu_ur+21])
plt.colorbar()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19586092715231787, 'fraction_numerical': 0.039735099337748346, 'mean_word_length': 5.818284424379232, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 12, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '8158', 'license': 'MIT', 'n_tokens_mistral': 2928, 'n_tokens_neox': 2817, 'n_words': 374, 'path': 'notebooks/development/preprocessing/make_grid.ipynb', 'repo_name': 'kuechenrole/tidal_melting', 'size': '6040'}