<jupyter_start><jupyter_text>Totten Ice Shelf geometry<jupyter_code>import sys
import cartopy
import shapefile
import numpy as np
import pyproj
import xarray as xr
import pandas as pd
import matplotlib
import cmocean.cm as cmo
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

sys.path.append("..")
%matplotlib inline
# %matplotlib widget
# %matplotlib notebook
%config InlineBackend.print_figure_kwargs={'bbox_inches':None}
%load_ext autoreload
%autoreload 2

# from geometry import ModelGeometry
ds = xr.open_dataset('../../data/BedMachine/BedMachineAntarctica_2020-07-15_v02.nc')
ds.mask
%%time
ds.mask.where(ds.mask>0).plot(color='darkblue')
ds.mask.where(ds.mask==3).plot(color='lightblue')
from geometry import xlim, ylim
for glacier in list(xlim.keys()):
    plt.figure()
    ds_ = ds['mask'].sel({'x':xlim[glacier], 'y':ylim[glacier]})
    ds_.plot()<jupyter_output><empty_output><jupyter_text>Bedmachine topography data<jupyter_code>ds
%%time
plt.figure(figsize=(10,10))
(ds.surface-ds.thickness).where(ds.mask==3).plot()<jupyter_output><empty_output><jupyter_text>selecting Totten Ice Shelf<jupyter_code>xlim = slice(2.19e6,2.322e6)
ylim = slice(-1e6,-1.3e6)
lim = {'x':xlim, 'y':ylim}
mask = ds['mask'].sel(lim)<jupyter_output><empty_output><jupyter_text>create grounding line and ice front masks<jupyter_code>grl = xr.where(mask-mask.shift(x= 1)==1, mask, 0) + \
      xr.where(mask-mask.shift(x=-1)==1, mask, 0) + \
      xr.where(mask-mask.shift(y= 1)==1, mask, 0) + \
      xr.where(mask-mask.shift(y=-1)==1, mask, 0)
grl = grl/grl
grl.name = 'grounding_line'

isf = xr.where(mask-mask.shift(x= 1)==3, mask, 0) + \
      xr.where(mask-mask.shift(x=-1)==3, mask, 0) + \
      xr.where(mask-mask.shift(y= 1)==3, mask, 0) + \
      xr.where(mask-mask.shift(y=-1)==3, mask, 0)
isf = isf/isf
isf.name = 'iceshelf_front'
# pick points separating grounding line from pinning points
line = [(2.1937e6,-1.2807e6),
        (2.2077e6,-1.2057e6), 
        (2.2234e6,-1.2121e6), 
        (2.2277e6,-1.2121e6), 
        (2.2334e6,-1.1725e6), 
        (2.2429e6,-1.1593e6), 
        (2.2502e6,-1.1075e6), 
        (2.2627e6,-1.1068e6), 
        (2.2728e6,-1.0617e6),
        (2.2974e6,-1.1226e6),
        (2.3199e6,-1.0849e6),
       ]
x_, y_ = [], []
for (x,y) in line:
    x_.append(x)
    y_.append(y)
grl_ = ModelGeometry.remove_points(grl, line)

line = [(2.3030e6,-1.1333e6),
        (2.3179e6,-1.1074e6),
       ]
isf_ = ModelGeometry.remove_points(isf, line, direction='above')
plt.figure(figsize=(10,6), constrained_layout=True)
mask.plot(alpha=.5, add_colorbar=False)
plt.plot(x_,y_)
grl_.where(grl_>0).plot(cmap='Blues', add_colorbar=False)
isf_.where(isf_>0).plot(cmap='Reds' , add_colorbar=False)
X, Y = np.meshgrid(mask.x.values, mask.y.values)
x = xr.DataArray(data=X, dims=['y','x'], coords={'y':mask.y,'x':mask.x})
y = xr.DataArray(data=Y, dims=['y','x'], coords={'y':mask.y,'x':mask.x})

mask_a = xr.where(mask==3, 1, 0)  # ice shelf
mask_b = xr.where(grl_==1, 1, 0)  # grounding line
mask_c = xr.where(isf ==1, 1, 0)  # ice shelf front
mask_c.plot()
%%time
# distance to grounding line
di = ModelGeometry.distance_to_line(mask_a,mask_c)
dg = ModelGeometry.distance_to_line(mask_a,mask_b)
divnorm = matplotlib.colors.DivergingNorm(vmin=-3000., vcenter=0, vmax=500)
kwargs = {'cbar_kwargs':{'orientation':'horizontal'}}
f, ax = plt.subplots(2, 3, figsize=(15,10), constrained_layout=True, sharey=True, sharex=True)

ds.surface.where(ds.surface>0).sel(lim).plot(ax=ax[0,0], **kwargs)
ds.bed.sel(lim).plot(ax=ax[0,1], cmap=cmo.balance, norm=divnorm, **kwargs)
ds.thickness.where(ds.thickness>0).sel(lim).plot(ax=ax[0,2], cmap='magma_r', **kwargs)
da_ = (ds.surface-ds.thickness).where(ds.mask==3).sel(lim)
da_.name = 'draft [meters]'
da_.plot(ax=ax[1,0], cmap='plasma', **kwargs, vmax=0)
grl_.where(grl_>0).plot(ax=ax[1,0], cmap='Blues', add_colorbar=False)
isf_.where(isf_>0).plot(ax=ax[1,0], cmap='Reds'   , add_colorbar=False)
di.name = 'to ice shelf front [km]'
(di/1e3).plot(ax=ax[1,1], **kwargs)
dg.name = 'to grounding line [km]'
(dg/1e3).plot(ax=ax[1,2], **kwargs)

f.suptitle('Totten Ice Shelf', fontsize=18)
plt.savefig('../../results/Bedmachine/TottenIS_geometry')<jupyter_output><empty_output><jupyter_text>add lat/lon coords and write to netcdf<jupyter_code>da_.name='draft'
grl_.name='grl'
isf_.name='isf'
da = xr.merge([da_, xr.where(mask==3,1,0), xr.where(grl_==1,1,0), xr.where(isf_==1,1,0)])

project = pyproj.Proj("epsg:3031")
xx, yy = np.meshgrid(da.x, da.y)
lons, lats = project(xx, yy, inverse=True)
dims = ['y','x']
da = da.assign_coords({'lat':(dims,lats), 'lon':(dims,lons)})

da.to_netcdf(f'../../results/Totten.nc')<jupyter_output><empty_output><jupyter_text>projection of bed data<jupyter_code>colors_undersea = cmo.deep_r(np.linspace(0, 1, 256))
colors_land = matplotlib.cm.cividis_r(np.linspace(0.1, .5, 256))
all_colors = np.vstack((colors_undersea, colors_land))
terrain_map = matplotlib.colors.LinearSegmentedColormap.from_list('terrain_map', all_colors)
divnorm = matplotlib.colors.DivergingNorm(vmin=-6000., vcenter=0, vmax=2000)

# example plot
X = np.random.rand(5,5)*4000-2000
plt.pcolormesh(X, cmap=terrain_map, norm=divnorm)
plt.colorbar()
%%time
plt.figure(figsize=(10,10))
ax = plt.axes(projection=ccrs.Orthographic(central_latitude=-90, central_longitude=0))
ax.pcolormesh(ds.x,ds.y,ds.bed, cmap=terrain_map, norm=divnorm,
              transform=ccrs.SouthPolarStereo(true_scale_latitude=-71),
             )
ax.coastlines()
ax.gridlines()
plt.colorbar()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19286871961102106, 'fraction_numerical': 0.05780659103187466, 'mean_word_length': 5.282805429864253, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 32, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '405083', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2565, 'n_tokens_neox': 2392, 'n_words': 350, 'path': 'src/notebooks/Bedmachine.ipynb', 'repo_name': 'AJueling/melt', 'size': '5553'}