<jupyter_start><jupyter_text>get cube from point pX, pYshorter/quicker version than manually selecting using matplotlib widget<jupyter_code># to skip the above - having to manuall select point...
from datacube_tools import DATACUBETOOLS as dctools
import xarray as xr
import pandas as pd

instance = dctools()

pX = -2.5011e5
pY = -2.87188e6

cube_feature,_ = dctools.find_datacube_catalog_entry_for_point(instance, point_xy = [pX,pY], point_epsg_str = '3413');

# for zarr store modify URL for use in boto open - change http: to s3: and lose s3.amazonaws.com
incubeurl = (
    cube_feature["properties"]["zarr_url"]
    .replace("http:", "s3:")
    .replace(".s3.amazonaws.com", "")
)
print("found the cube")
cube = xr.open_dataset(incubeurl, engine = "zarr", storage_options = {"anon" : True})
print("finished loading cube @ ", pd.Timestamp.now())<jupyter_output><empty_output><jupyter_text>the above is NOT my (ed) handiwork: the below isimport libraries<jupyter_code>import matplotlib.pyplot as plt
import matplotlib as mpl
font = {'size'   : 18}
mpl.rc('font', **font)
import seaborn as sns
import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon
import cartopy
import cartopy.crs as ccrs
import sys
from glob import glob
import os
import numpy as np
import rioxarray
import xarray as xr
import statsmodels.stats.api as sms
from scipy import stats
from cycler import cycler
import datetime as dt
from tqdm import tqdm
from itertools import product
import warnings
warnings.filterwarnings('ignore')

import helper_functions as hf
import importlib
importlib.reload(hf)<jupyter_output><empty_output><jupyter_text>clip single time slice from cube to smaller regioneither a draw ROI, or buffered area around centreline<jupyter_code>%%time
# takes ~10 minutes
# read in geojson of polygon of interest from file
# folder = '../../../qgis_work/polygons_for_itslive/*Region.geojson' # this version is saved my working QGIS folder
folder = '../fromQGIS/*Region.geojson' # this version is saved git directory
paths = glob(folder)

# # read in first file ([0]) in directory only - option to iterate through multiple...later
# # AOI = Region of Interest
# AOI = gpd.read_file(paths[0])
# # this is a simple polygon

# # select second polygon in AOI
# AOI = gpd.GeoDataFrame(geometry = AOI.loc[1], crs = 'epsg:3413').reset_index(drop = True)

# buffer area around centreline and use that to clip cube
centreLine = gpd.read_file('../fromQGIS/isortuarsuupTasia_velocityCentreLine.geojson')
AOI = gpd.GeoDataFrame(geometry = centreLine.loc[[1]].buffer(1000), crs = 'epsg:3413').reset_index(drop = True)

clipped = hf.clipcube(cube,AOI) # this function works - but could certainly be improved
print("finished clipping @ ", pd.Timestamp.now())

# identify time slices where v values are all nan
allnans = np.count_nonzero(np.isnan(clipped.v), axis = (1,2)) == (clipped.v.shape[1] * clipped.v.shape[2])
clipped = clipped.isel(mid_date = np.nonzero(~allnans)[0])
print('removed timeslices with all nans')<jupyter_output><empty_output><jupyter_text>along a line- spatially average velocities within a 3x3 (5x5, or 7x7) window centred on each vertex- construct temporal averages of velocity profiles: error weighted by number and accuracy of measurements in time period<jupyter_code>line = centreline.loc[[1]]]
filterdays = 45
vertexspacing = 300
windowsize = 3

# create centreline xarrays
cl, k_vs = hf.centreline_sample(line = geometry,
                                ds = clipped,
                                sample_freq = vertexspacing,
                                kernal = windowsize,
                                glacier='Name',
                                glaciertype='Lake')
# resampled 
resampled_CL = hf.resample_cl(cl,filterdays = filterdays)
fig, axs = plt.subplots(nrows = 2, ncols = 3, figsize = [20,10], sharex = True, sharey = True)

for j, ind in enumerate([1271, 9875]): # 15059, 286,, 286 , 1271
    
    for i in [0,1,2]:
        # plot the resultant velocity at each grid square in kernal
        underlying = axs[j,i].scatter(np.tile(cumlDist.values, (window**2,1)), k_vs[ind,:,:], c = 'lightgrey', s = 5, zorder = 1, label = 'v estimates from window')# label = 'individual points in kernal');
    
        # plot the smoothed velocity, determined by averaging the vx and vy components in each window, and then computing resultant velocity
        smoothed = cl.isel(time = ind).avg_v.plot(ax = axs[j,i], label = 'mean v: 5x5 window', zorder=2, c = 'k')
        
        crudeavg = axs[j,i].plot(cumlDist, np.nanmean(k_vs[ind,:,:], axis = 0), c = 'k', linestyle = ':', label = 'average v (not using components)', zorder = 1)
            
    # show uncertainty about window averaged velocity, using formal propagation of errors
    verrfill = axs[j,0].fill_between(x = cumlDist,
                                  y1 = cl.isel(time = ind).avg_v.values + cl.isel(time = ind).avg_v_err.values,
                                  y2 = cl.isel(time = ind).avg_v.values - cl.isel(time = ind).avg_v_err.values,
                                  color = 'tab:orange', alpha = 0.5, zorder = 1.5, label = 'mean v ± ITS_LIVE average v_x/y_error')

    # verrfill = axs[j,0].fill_between(x = cumlDist,
    #                               y1 = cl.isel(time = ind).avg_v.values + cl.isel(time = ind).alt_v_err.values,
    #                               y2 = cl.isel(time = ind).avg_v.values - cl.isel(time = ind).alt_v_err.values,
    #                               color = 'tab:purple', alpha = 0.5, zorder = 1.5, label = 'mean v ± ITS_LIVE SE v_x/y_error')

    # stderrfill = axs[j,1].fill_between(x = cumlDist,
    #                               y1 = cl.isel(time = ind).avg_v.values + cl.isel(time = ind).avg_v_err_fromstd.values,
    #                               y2 = cl.isel(time = ind).avg_v.values - cl.isel(time = ind).avg_v_err_fromstd.values,
    #                               color = 'tab:blue', alpha = 0.5, zorder = 1.5, label = 'mean v ± error using using std(v_x/y)')
    
    # CIfill = axs[j,2].fill_between(x = cumlDist,
    #                               y1 = cl.isel(time = ind).avg_v.values + cl.isel(time = ind).avg_v_err_fromCI.values,
    #                               y2 = cl.isel(time = ind).avg_v.values - cl.isel(time = ind).avg_v_err_fromCI.values,
    #                               color = 'tab:pink', alpha = 0.5, zorder = 1.5, label = 'mean v ± uncert from 95% CI of v_x/y')

    CIfill = axs[j,0].fill_between(x = cumlDist,
                                  y1 = cl.isel(time = ind).avg_v.values + cl.isel(time = ind).avg_v_err_fromSE.values,
                                  y2 = cl.isel(time = ind).avg_v.values - cl.isel(time = ind).avg_v_err_fromSE.values,
                                  color = 'tab:cyan', alpha = 0.5, zorder = 3.5, label = 'mean v ± uncert from SE of v_x/y')
    
    for i in [0,1,2]:
        axs[0,i].legend(fontsize = 14, frameon = False, loc = 'upper center')
        td = (pd.to_timedelta(cl.isel(time = ind).date_dt.values)).components
        axs[j,i].set_title('Mid-date: ' + axs[j,i].get_title()[7:17] + ' / date_dt: ' + str(td.days) + ' days')

for ax in axs.flatten():
    ax.set_xlabel("")
    ax.set_ylabel("")

# for ax in axs.flatten()[3:]:
#     ax.set_title("")
    
axs[0,0].set_ylabel('Velocity (m/yr)')
axs[1,0].set_ylabel('Velocity (m/yr)')

for ax in axs.flatten()[-4:]:
    ax.set_xlabel('Distance (m)')

hf.labelFig(fig, pos = 'topleft')
plt.tight_layout()
# fig.savefig('../figures/isortuarsuupTasia/methods_of_spatial_averaging.png', bbox_inches = 'tight', dpi = 300)<jupyter_output><empty_output><jupyter_text>resampling to seasonal/monthly<jupyter_code># filtering velocities along centreline, by image separation: date_dt < 90/60/30 days.
# and resampling to seasonal
# - set **FILTERDAYS**

FILTERDAYS = 45 # np.max(cl.date_dt.values.astype('timedelta64[D]').astype(int))

## resampling to seasonal
ssnl_45 = hf.resample_cl(cl,filterdays = FILTERDAYS)
# for exporting to csv: ssnl_45.to_dataframe().to_csv(*****.csv')

##  resampling to seasonal: unfiltered (i.e. all readings)
sep = cl.date_dt.values.astype('timedelta64[D]').max().astype('int') + 1
ssnl_all = hf.resample_cl(cl, filterdays = sep)<jupyter_output><empty_output><jupyter_text>plotting velocity along centreline- seasonal averages, through time- filtering on image pair separation (date_dt)<jupyter_code>## plot velocity along centreline per season per year

# # for getting color bar nonsense because 
cNorm = mpl.colors.Normalize(vmin = np.min(ssnl_45.time.dt.year.values),# can use `astype('float64')` instead of `.dt.year` to vary colors between seasons
                             vmax = np.max(ssnl_45.time.dt.year.values))

fig, axs = plt.subplots(ncols = 4, nrows = 1, sharex = True, sharey = True, figsize = (18,6))

# need to loop through each season/year to have full control over colormap
for i, s in enumerate(['DJF', 'MAM', 'JJA', 'SON']):
    for t in ssnl_45.isel(time = ssnl_45.time.dt.season.values == s).time.dt.year.values:
        
        # get index
        idx = ((ssnl_45.time.dt.season.values == s) & (ssnl_45.time.dt.year.values == t))
        if np.count_nonzero(idx) == 0:
            continue
        # plot average line
        ssnl_45.isel(time = idx).v.plot(ax = axs.flat[i],
                                                           c = mpl.cm.viridis(cNorm(ssnl_45.isel(time = idx).time.dt.year))) # can use `astype('float64')` instead of `.dt.year` to vary colors between seasons
        
        # to plot the fill between uncertainty
        axs.flat[i].fill_between(x = ssnl_45.isel(time = idx).distance,
                                 y1 = (ssnl_45.isel(time = idx).v + ssnl_45.isel(time = idx).verr).values.flatten(),
                                 y2 = (ssnl_45.isel(time = idx).v - ssnl_45.isel(time = idx).verr).values.flatten(),
                                 alpha = 0.3,
                                #  color = mpl.cm.viridis(cNorm(ssnl_45.isel(time = idx).time.dt.year))) # can use `astype('float64')` instead of `.dt.year` to vary colors between seasons
                                 color = mpl.cm.viridis(cNorm(t)))
        axs.flat[i].set_title(s)
        
for ax in axs.flat:
    ax.set_ylim([0,400])
    ax.set_xlim([-200,15000])
    ax.set_xticks([0,4000,8000,12000])
    ax.set_xticklabels((ax.get_xticks() / 1000).astype(int))
    ax.set_xlabel('Distance (km)')
    ax.set_ylabel('')
    
axs[0].set_ylabel('Velocity (m/yr)')

# add axes labels - do this before adding colorbar axis
hf.labelFig(fig);

################# for color bar
# [left, bottom, width, height]
cbax = fig.add_axes([0.92,0.2,0.15,0.5])
cbax.axis('off')

# cNorm = mpl.colors.Normalize(vmin = minyear,
#                              vmax = maxyear)
cmap = mpl.cm.get_cmap('viridis')
sc = mpl.cm.ScalarMappable(norm = cNorm, cmap = cmap)

# for correctly labelling colorbar
tickys = pd.date_range(np.min(ssnl_45.time.values) - pd.offsets.DateOffset(years=1) ,
              np.max(ssnl_45.time.values) + pd.offsets.DateOffset(years=1) ,
              freq = 'AS')

cb = plt.colorbar(sc,
                  ax=cbax,
                  orientation = 'vertical',
                  fraction = 1,
                  ticks = tickys.year, # values.astype('float64'),
                  label = 'Year'
                  )

cb.set_ticklabels(tickys.year)

plt.suptitle('Seasonal velocity profiles Image pair separation <=' + str(FILTERDAYS)+'days');
# plt.suptitle('Seasonal velocity profiles at Isortuarsuup Tasia. Image pair separation <=' + str(FILTERDAYS)+'days');

fig.subplots_adjust(wspace=0.06, hspace=0)
# fig.savefig('../figures/isortuarsuupTasia/ssnl_cntrln_v_errFromWindowSE_filtered_'+str(FILTERDAYS)+'days.png', bbox_inches='tight') # for ssnl_45
# fig.savefig('../figures/landterminating/seasonal_velocity_along_centreline_filtered_'+str(FILTERDAYS)+'days.png', bbox_inches='tight')
fig, axs = plt.subplots(ncols = 4, nrows = 1, sharex = True, sharey = True, figsize = (18,6))

for i, s in enumerate(['DJF', 'MAM', 'JJA', 'SON']):
    for t in ssnl_45.isel(time = ssnl_45.time.dt.season.values == s).time.dt.year.values:

        # get index
        idx = ((ssnl_45.time.dt.season.values == s) & (ssnl_45.time.dt.year.values == t))
        if np.count_nonzero(idx) == 0:
            continue

        ssnl_45.isel(time = idx)['n'].plot(c = mpl.cm.viridis(cNorm(ssnl_45.isel(time = idx).time.dt.year)), ax = axs[i])
        if i > 0:
            axs[i].set_ylabel('')
        axs[i].set_title(s)
        axs[i].set_yscale('log')
        axs[i].set_xlim([-200,15000])
        axs[i].set_xticks([0,4000,8000,12000])
        axs[i].set_xticklabels((ax.get_xticks() / 1000).astype(int))
        axs[i].set_xlabel('Distance (km)')


axs[0].set_ylabel('Count')

hf.labelFig(fig)

################# for color bar
# [left, bottom, width, height]
cbax = fig.add_axes([0.92,0.2,0.15,0.5])
cbax.axis('off')

# cNorm = mpl.colors.Normalize(vmin = minyear,
#                              vmax = maxyear)
cmap = mpl.cm.get_cmap('viridis')
sc = mpl.cm.ScalarMappable(norm = cNorm, cmap = cmap)

# for correctly labelling colorbar
tickys = pd.date_range(np.min(ssnl_45.time.values) - pd.offsets.DateOffset(years=1) ,
              np.max(ssnl_45.time.values) + pd.offsets.DateOffset(years=1) ,
              freq = 'AS')

cb = plt.colorbar(sc,
                  ax=cbax,
                  orientation = 'vertical',
                  fraction = 1,
                  ticks = tickys.year, # values.astype('float64'),
                  label = 'Year'
                  )

cb.set_ticklabels(tickys.year)

plt.suptitle('Count of velocity measurements used to construct seasonal averages (<=' + str(FILTERDAYS)+'days)');
# plt.suptitle('Seasonal velocity profiles at Isortuarsuup Tasia. Image pair separation <=' + str(FILTERDAYS)+'days');

fig.subplots_adjust(wspace=0.06, hspace=0)

## 3 subplots: 
## 1. all velocity profiles where mid_date is in given month/season
## 2. velocity profiles in given month season where date_dt < FILTERDAYS (~90/60/30 days), with uncertainty
## 3. error weighted average of 1. (red) and 2. (blue). each with propagated uncertainty - which is a bit too narrow. Have included ± 1 standard deviation about blue line

testS = 'JJA'
testY = 2021

fig, axs = plt.subplots(ncols = 3, figsize = [20,6], sharey = True)

################# AXES 0 ##################
# plot velocities along centreline for all timesteps (UNFILTERED by date_dt)
idx = (cl.time.dt.season.values == testS) & (cl.time.dt.year.values == testY)
cl.isel(time = idx).avg_v.plot(x = 'distance', hue = 'time', add_legend = False,
                                linewidth = 0.5, alpha = 0.8, ax = axs[0]);

axs[0].text(x = axs[0].get_xticks()[-2],
            y = axs[0].get_yticks()[-2],
            s = 'n = '+str(cl.isel(time = idx).avg_v.time.shape[0]),
            horizontalalignment = 'right')         

axs[0].set_ylabel("Velocity (m/yr)")
axs[0].set_xlabel("Distance (m)")
axs[0].set_title('Unfiltered - all profiles')
################# AXES 1 ##################

# plot velocities along centreline for all timesteps (filtered by date_dt)
idx = (cl_filtered.time.dt.season.values == testS) & (cl_filtered.time.dt.year.values == testY)
cl_filtered.isel(time = idx).avg_v.plot(x = 'distance', hue = 'time', add_legend = False,
                                         linewidth = 1, alpha = 0.8, ax = axs[1]);

axs[1].text(x = axs[1].get_xticks()[-2],
            y = axs[1].get_yticks()[-2],
            s = 'n = '+str(cl_filtered.isel(time = idx).avg_v.time.shape[0]),
            horizontalalignment = 'right')

## plot velocity uncertainty for each centreline
for i, timey in enumerate(cl_filtered.isel(time = idx).time.values):
    
    idx = cl_filtered.time == timey
    
    if cl_filtered.isel(time = idx).time.shape[0] == 1:
            
        axs[1].fill_between(x = cl_filtered.isel(time = idx).distance,                            
                            y1 = (cl_filtered.isel(time = idx).avg_v +
                                  cl_filtered.isel(time = idx).avg_v_err).values.flatten(),                            
                            y2 = (cl_filtered.isel(time = idx).avg_v -
                                  cl_filtered.isel(time = idx).avg_v_err).values.flatten(),
                            alpha = 0.1)
    else: continue
axs[1].set_ylabel("")
axs[1].set_xlabel("Distance (m)")
axs[1].set_ylim(axs[0].get_ylim());
axs[1].set_title('Filtered profiles: date_dt < '+str(FILTERDAYS)+'days')

################# AXES 2 ##################
############################################# unfiltered #####################################################
# plot error weighted average
ssnl_all.isel(time = ((ssnl_all.time.dt.season.values == testS) & (ssnl_all.time.dt.year.values == testY))).v.plot(x = 'distance',
                                                                                                                         linewidth = 0.5,linestyle = '--',c = 'k',
                                                                                                                         ax = axs[2]); #label = 'unfiltered avg',
# plot uncertainty
Sidx = (ssnl_all.time.dt.season.values == testS) & (ssnl_all.time.dt.year.values == testY)
axs[2].fill_between(x = ssnl_all.isel(time = Sidx).distance,                    
                    y1 = (ssnl_all.isel(time = Sidx).v + 
                          ssnl_all.isel(time = Sidx).verr).values.flatten(),                    
                    y2 = (ssnl_all.isel(time = Sidx).v -
                          ssnl_all.isel(time = Sidx).verr).values.flatten(),
                   label = 'unfiltered',color = 'r')
####################################### filtered ################################################
# plot error weighted average
ssnl_45.isel(time = ((ssnl_45.time.dt.season.values == testS) &
                               (ssnl_45.time.dt.year.values == testY))).v.plot(x = 'distance',
                                                                                                  linewidth = 1,linestyle = ':',c = 'k',
                                                                                                  ax = axs[2]); # label = 'filtered average',
# plot uncertainty
SFidx = (ssnl_45.time.dt.season.values == testS) & (ssnl_45.time.dt.year.values == testY)
axs[2].fill_between(x = ssnl_45.isel(time = SFidx).distance,
                    y1 = (ssnl_45.isel(time = SFidx).v + 
                          ssnl_45.isel(time = SFidx).verr).values.flatten(),
                    y2 = (ssnl_45.isel(time = SFidx).v -
                          ssnl_45.isel(time = SFidx).verr).values.flatten(),
                    alpha = 0.5,
                    label = 'filtered - uncert')

# ± 1 standard deviation of filtered centrelines and plotting either side of seasonal average
SFidx = (ssnl_45.time.dt.season.values == testS) & (ssnl_45.time.dt.year.values == testY)
axs[2].fill_between(x = ssnl_45.isel(time = SFidx).distance,
                    y1 = (ssnl_45.isel(time = SFidx).v + 
                          ssnl_45.isel(time = SFidx).v_std).values.flatten(),
                    y2 = (ssnl_45.isel(time = SFidx).v -
                          ssnl_45.isel(time = SFidx).v_std).values.flatten(),
                    alpha = 0.2, color = 'cornflowerblue',
                    label = 'filtered: ± 1std')

plt.legend()
axs[2].set_ylabel("")
axs[2].set_xlabel("Distance (m)")
axs[2].set_ylim(axs[0].get_ylim());
axs[2].set_title('Seasonal average');#, loc = 'left', y = 1.0,  pad = -30);

fig.suptitle(str(testS) + ' ' +  str(testY))

fig.subplots_adjust(wspace=0.05, hspace=0)
hf.labelFig(fig);
# fig.savefig('../figures/isortuarsuupTasia/filtering_by_datedt_'+str(FILTERDAYS)+'days_and_averaging_'+testS+str(testY)+'.png',bbox_inches='tight')
# fig.savefig('../figures/landterminating/filtering_by_datedt_'+str(FILTERDAYS)+'days_and_averaging_'+testS+str(testY)+'.png',bbox_inches='tight')
## same as #3 above, but for all month/years in record
fig, axs = plt.subplots(nrows = 10, ncols = 4, figsize = [20,20], sharex = True)

# loop through each season/year
theSeasons = ['DJF', 'MAM', 'JJA', 'SON']
yrs = pd.unique(ssnl_all.time.dt.year)

for i, S in enumerate(theSeasons):
    for j, Y in enumerate(yrs):
      if ssnl_all.isel(time = ((ssnl_all.time.dt.season.values == S) & (ssnl_all.time.dt.year.values == Y))).time.shape[0] == 0: # if there are no measurements for given season/year, skip it
            continue
      ############################################# unfiltered #####################################################
      # plot error weighted average
      # ssnl_all.isel(time = ((ssnl_all.time.dt.season.values == S) &
      #                       (ssnl_all.time.dt.year.values == Y))).v.plot(x = 'distance',
      #                                                                             linewidth = 0.5, linestyle = '--', c = 'k',
      #                                                                             label = 'unfiltered avg',
      #                                                                             ax = axs[j,i]);
      
      # get index for season and year
      idx = (ssnl_all.time.dt.season.values == S) & (ssnl_all.time.dt.year.values == Y)

      # plot uncertainty about error weighted average
      axs[j,i].fill_between(x = ssnl_all.isel(time = idx).distance,
                        y1 = (ssnl_all.isel(time = idx).v + 
                              ssnl_all.isel(time = idx).verr).values.flatten(),
                        y2 = (ssnl_all.isel(time = idx).v -
                              ssnl_all.isel(time = idx).verr).values.flatten(),
                        label = 'unfiltered uncert',color = 'r',
                        alpha = 0.5)
      
      ####################################### filtered ################################################
      # plot error weighted average
      # ssnl_45.isel(time = ((ssnl_45.time.dt.season.values == S) &
      #                                (ssnl_45.time.dt.year.values == Y))).v.plot(x = 'distance',
      #                                                                                               linewidth = 2,linestyle = ':',c = 'k',
      #                                                                                               label = 'filtered average',
      #                                                                                               ax = axs[j,i]);
      # get index 
      idx = (ssnl_45.time.dt.season.values == S) & (ssnl_45.time.dt.year.values == Y)
      # plot uncertainty
      axs[j,i].fill_between(x = ssnl_45.isel(time = idx).distance,
                        y1 = (ssnl_45.isel(time = idx).v + 
                              ssnl_45.isel(time = idx).verr).values.flatten(),
                        y2 = (ssnl_45.isel(time = idx).v -
                              ssnl_45.isel(time = idx).verr).values.flatten(),
                        alpha = 0.8,
                        label = 'filtered uncert')
      
      # ± 1 standard deviation of filtered centrelines and plotting either side of seasonal average
      axs[j,i].fill_between(x = ssnl_45.isel(time = idx).distance,
                        y1 = (ssnl_45.isel(time = idx).v + 
                              ssnl_45.isel(time = idx).v_std).values.flatten(),
                        y2 = (ssnl_45.isel(time = idx).v -
                              ssnl_45.isel(time = idx).v_std).values.flatten(),
                        alpha = 0.5, color = 'cornflowerblue',
                        label = 'filtered: ± 1std')

      # removing auto placed by xarray
      axs[j,i].set_title("")
      
      # using annotate to set title - because otherwise bottom row titles do not go into same position as others
      axs[j,i].annotate(str(S)+' '+str(Y), xy=(0.5,0.9),  xycoords='axes fraction',
                        xytext=(0.5, 0.9), textcoords='axes fraction',
                        horizontalalignment='center', verticalalignment='top')
      
      axs[j,i].set_xlabel("") # remove xlabels
      axs[j,i].set_ylabel("")
      
      # add legend to final plot
      if (j == 9) & (i == 1):
            axs[j,i].legend(bbox_to_anchor = (2,1.1))
            
# remove extra axes
axs[0,0].get_xaxis().set_visible(False)
axs[0,0].get_yaxis().set_visible(False)

axs[0,0].annotate('filter image separtion\n <='+str(FILTERDAYS)+' days',
                  xy = (0.5,0.5), xycoords='axes fraction',
                  xytext=(0.5, 0.5), textcoords='axes fraction',
                  horizontalalignment='center', verticalalignment='center')
fig.delaxes(axs[9,2])
fig.delaxes(axs[9,3])

# fig.savefig('../figures/seasonalAvg_unfiltered_vs_filtered_'+str(FILTERDAYS)+'days_wSTD.png',bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>velocities in a box/circle- spatially average velocities within box(es) at specified distance along centreline- take error weighted monthly averages and plot annual cycles.<jupyter_code># specify distances up glaicer to sample around
distanceUpG = [1500,3000,10000]
# and size of radius to buffer around
Radius = 500

P = []
clidx = []
distance = []  

# get x,y coordinates from centreline at distances closest to distanceUpG
for d in distanceUpG:
    
    idx = np.abs(cl.distance - d).argmin(dim = 'distance').values
    clidx.append(idx) # keep track of numerical index from cl.distance
    
    # create shapely Point each vertex that is closest to specified distanceUpG
    P.append(Point(cl['x'][idx].values,
                   cl['y'][idx].values))
                   
    distance.append(cl['distance'][idx].values.tolist()) # keep track of actual distance along centreline (which will approximate distanceUpG)
    
# create geodataframe of areas to sample - with the same number of rows as number of distances specified ↑.
# buffer each point by specified radius, Radius → gives polygons
circles = gpd.GeoDataFrame({'geometry' : gpd.GeoSeries(P).buffer(Radius),
                            'distance' : distance,
                            'clidx' : clidx},
                          crs = 'epsg:3413')
# fig, ax = plt.subplots()
# circles.plot(ax = ax, fc = 'None', column = 'distance', cmap = 'Dark2')
# centreLine.loc[[0],:].plot(ax = ax, color = 'k')
# fig.savefig('../figures/velocity_sample_regions.png', bbox_inches='tight')

from matplotlib.path import Path

# empty dataframes for concat'ing results of each polygon into
rois = pd.DataFrame() # contains spatially averaged velocities within each buffer region, for all timesteps (exc. nans, and where image separation > 30 days)
monthly = pd.DataFrame() # resampled rois to monthly: contains monthly averaged velocities for each each region
 
# get list of all coordinate pairs in clipped grid
xs,ys = np.meshgrid(clipped.x, clipped.y)
xs, ys = xs.flatten(), ys.flatten()
points = np.vstack((xs,ys)).T

cV = [] # empty list for putting ALL (in space and time) velocity measurements into
ddt = []
ts = []

# clip clipped cube to within single circle
for i in tqdm(range(len(circles))):    
    
    ## at present this method results in an inconsistent number of gridsquares being sampled - 
    ## depending on where the sample point is relative to the centre of a grid square
    
    ## (relatively) quick way of clipping, as opposed to: `clipped.rio.clip(circle)`
    x,y = ((circles.loc[[i],'geometry'][i].exterior.coords.xy)) # vertices of polygon...
    verts = list(zip(x,y)) # ..as a list of tuples
    p = Path(verts) # create polygon (matplotlib path array) from vertices
    grid = p.contains_points(points) # returns logical for whether points are in
    xc, yc = map(list, zip(*points[grid])) # index into original list of points, and get coordinate pairs for points that in polygon
    
    # query xarray at point
    circle = clipped.sel(x = pd.unique(xc), y = pd.unique(yc))
    
    #################################################################
    # for getting distribution of velocities WITHIN kernal at all timestemps
    circleV = circle['v'].stack(dimensions = {'xy' :('x','y')}) # this reduces the 'x','y' dims to one.
    circleV = circleV.reset_index(dims_or_levels = 'xy')
    
    # list of image_separation values
    date_dt = circle['date_dt'].values
    ddt.append(date_dt)
    
    # list of image mid_dates
    timestamps = pd.to_datetime(circleV['mid_date'].values)
    ts.append(timestamps)

    cV.append(circleV.values)
    #################################################################    
    
    # average velocities across space
    # averge velocity compoentns in each kernal, ignore nans
    avg_vx = circle.vx.mean(dim = ['x','y'], skipna = True).values
    avg_vy = circle.vy.mean(dim = ['x','y'], skipna = True).values

    # number of valid measurements in each kernal
    n = circle.v.count(dim = ['x','y']).values

    # uncertainty in average velocity compnent value
    avg_vx_err = circle.vx_error / np.sqrt(n)
    avg_vy_err = circle.vx_error / np.sqrt(n)

    # calculate average resultant velocity and average uncertainty
    avg_v, avg_v_err = hf.v_uncert(avg_vx, avg_vy, avg_vx_err, avg_vy_err)

    # create dataframe, index = time, that contains velocities from all image pairs
    roi = pd.DataFrame({'date_dt' : circle.date_dt.values,
                        'avg_v' : avg_v, 'avg_v_err' : avg_v_err,
                        'avg_vx' : avg_vx, 'avg_vy' : avg_vy,
                        'avg_vx_err' : avg_vx_err, 'avg_vy_err' : avg_vy_err},
                       index = circle.mid_date)

    # remove nan-rows
    roi = roi.loc[~np.isnan(roi.avg_v),:]
    


    # add a month number and year column - for seaborn
    roi['Month'] = roi.index.month
    roi['Year'] = roi.index.year
    roi['Region_num'] = i
    roi['distance_upG'] = distanceUpG[i]
    roi['daysep'] = roi.date_dt.round('1D').astype('timedelta64[D]').astype(int)
    
    # add velocity observations in roi to dataframe containing all velocity measurements from all ROIS
    rois = pd.concat([rois, roi])

    
    # resample to MONTHLY w/ error weighted uncertainties
    # error weighted averages of vx and vyROI
    top = roi.avg_vx / (roi.avg_vx_err**2)
    bottom = 1 / (roi.avg_vx_err**2)
    Mavg_vx = (top.resample('M').sum() / bottom.resample('M').sum()).rename('Mavg_vx') # dimensions time

    # vy
    top = roi.avg_vy / (roi.avg_vy_err**2)
    bottom = 1 / (roi.avg_vy_err**2)
    Mavg_vy = (top.resample('M').sum() / bottom.resample('M').sum()).rename('Mavg_vy')

    # monthly average uncertainty
    Mavg_vx_err = (np.sqrt(1/ ((1 / (roi.avg_vx_err**2)).resample('M').sum()))).rename('Mavg_vx_err')
    Mavg_vy_err = (np.sqrt(1/ ((1 / (roi.avg_vy_err**2)).resample('M').sum()))).rename('Mavg_vy_err')

    # number of measurements within each month
    Mn = roi.avg_vx.resample('M').count().rename('MCount') # counting either vx or vy works
    
    # resultant velocity from vx and vy, and error weighted resultant uncertainty
    v_monthly, verr_monthly = hf.v_uncert(Mavg_vx, Mavg_vy, Mavg_vx_err, Mavg_vy_err)
    v_monthly = v_monthly.rename('v_monthly')
    verr_monthly = verr_monthly.rename('verr_monthly')

    # standard deviation of average resultant velocity (avg_v)
    v_std = (roi.avg_v.resample('M').std()).rename('v_std')

    # stack into monthly dataframe
    m = pd.DataFrame([v_monthly, verr_monthly, v_std, Mavg_vx, Mavg_vy, Mavg_vx_err, Mavg_vy_err, Mn]).T

    # and add year and month column
    m['Month'] = m.index.month
    m['Year'] = m.index.year
    m['Region_num'] = i
    m['distance_upG'] = distanceUpG[i]
    
    monthly = pd.concat([monthly, m])
## looks at distribution of velocity estimates for given region, for every month and year.
## prints the error weighted estimate of average velocity in that month - which is the value used in subsequent 'seasonal' plots
## along with the median (image average) and number of valid image-pairs used (i.e. sum of count on left hand axes)
## in red, on right hand axes, distribution of all resultant velocity measurements within the buffered region, with mid_date in that month/year, and date_dt < filter days == 30 days.

# select buffer region
region = 0

######### FOR SELECTION OF MONTHS/YEARS #################
data = rois[rois['Region_num'] == region].reset_index()
# facet grid plot af error weighted average resultant velocity, rows = year, columns = months
g = sns.displot(data = data.loc[data['Year'].isin([2018,2019]) & data['Month'].isin([3,4,5,6,7,8])],
                x = 'avg_v', col = 'Month', row = 'Year', hue = 'Region_num',
                kind = 'kde', common_norm = False, #element = 'step', #kde = True
                facet_kws={'sharey': False, 'sharex': True},
                legend = False)

######### FOR ALL MONTHS AND YEARS #########
# huge facet grid plot af error weighted average resultant velocity, rows = year, columns = months
# g = sns.displot(data = rois[rois['Region_num'] == region].reset_index(),
#                 x = 'avg_v', col = 'Month', row = 'Year', hue = 'Region_num',
#                 kind = 'kde', common_norm = False, #element = 'step', #kde = True
#                 facet_kws={'sharey': False, 'sharex': True},
#                 legend = False)

emptyaxes = []
# loop through each sub figure get year and month from title
for ax in g.axes.flatten():
    y = int(ax.properties()['title'][7:11])
    m = int(ax.properties()['title'][-2:])
    
    # index for rois df
    ind = (rois['Year'] == y) & (rois['Month'] == m) & (rois['Region_num'] == region)
    
    # number of valid image-pairs used
    N = len(rois.loc[ind,:])
    
    # skip if no image-pairs
    if N == 0:
        emptyaxes.append(ax)
        continue    

    # get median value (along time axis)
    med = rois.loc[ind,'avg_v'].median()
    
    # stats from monthly df
    mind = (monthly['Year'] == y) & (monthly['Month'] == m) & (monthly['Region_num'] == region)
    
    avg_v = monthly.loc[mind,'v_monthly'].values[0]
    err_v = monthly.loc[mind,'verr_monthly'].values[0]
    std_v = monthly.loc[mind,'v_std'].values[0]
    
    # smash the stats together and annotate
    string = ('N = ' + str(N) + '\nmean v =\n' + str(np.round(avg_v)) + '±' + str(np.round(err_v)) + '\nmedian =\n' + str(np.round(med)))

    ax.annotate(text = string,
                xy = (0.5,0.7),
                xytext = (0.65,0.5),
                xycoords = 'axes fraction',
                textcoords = 'axes fraction')
        
    # create index into cV on basis of month, year and date_dt
    tindx = ((ts[region].month == m) & (ts[region].year == y) & (ddt[region] < pd.Timedelta(f)))
    
    # create right hand axis
    # ax2 = ax.twinx();
        
    # plot histogram of underlying all underlying values in buffer regio non  secondary axis
    # ax2.hist(cV[region][tindx,:].flatten(),
    #      histtype = 'step', color = 'r', density = True);
    
    ## kde of all underlying values in buffer region - plot on same axes
    vals = (cV[region][tindx,:].flatten())
    vals_clean = vals[~np.isnan(vals)]
    kde = stats.gaussian_kde(vals_clean)
    xx = np.arange(vals_clean.min(), vals_clean.max())
    ax.plot(xx,kde(xx), c = 'r')    
    
    ############ plotting kde of monthly averages - this is the same as seaborn's built in kde plot
    ###### shade standard deviation region
    
    if N == 1:
        continue
    
    vals = rois.loc[ind,'avg_v']
    vals_clean = vals[~np.isnan(vals)]    
    kde = stats.gaussian_kde(vals_clean) # this is the same as the seaborn's built in kde plot
    
    # create shaded region with error weigthed uncertainty about average
    xx = np.arange(avg_v - err_v,
                   avg_v + err_v,
                   step = 0.5)
    yy = kde(xx)    
    ax.fill_between(x = xx,
                    y1 = yy,
                    y2 = np.zeros(len(yy)),
                    color = 'grey',
                    alpha = 0.5)
    
    # use this to create shaded region ± 1 std
    xx = np.arange(avg_v - std_v,
                   avg_v + std_v,
                   step = 0.5)
    yy = kde(xx)    
    ax.fill_between(x = xx,
                    y1 = yy,
                    y2 = np.zeros(len(yy)),
                    color = 'cornflowerblue',
                    alpha = 0.5)
    
    # ax.plot(xx,kde(xx), c = 'k', linewidth = 10, linestyle = ':')  

    # vertical lines at mean and median
    ax.vlines(avg_v,
              ymin = 0,
              ymax = kde(avg_v),
              color = 'k')
    ax.vlines(med,
              ymin = 0,
              ymax = kde(med),
              color = 'k',
              linestyle = ':')
    
    ax.set_xlim([-10,510])
    
for ax in emptyaxes:
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.set_title("")

for ax in g.axes.flatten()[-6:]:
    ax.set_xlabel('Velocity (m/yr)')
    
# g.savefig('../figures/isortuarsuupTasia/kde_of_velocites_used_to_calculate_monthly_averages_plusminus_1std_mean_and_median_with_underlying_pixel_distribution_in_red_at_1500_m_upglacier.png', bbox_inches = 'tight', dpi = 300)
# g.savefig('../figures/isortuarsuupTasia/velocity_kde_for_calc_mnthly_avgs_plusminus_1std_mean_and_median_with_underlying_pixel_dstrbtn_in_red_at_1500_m_upglacier_2018-19_MAMJJA.png', bbox_inches = 'tight', dpi = 300)<jupyter_output><empty_output><jupyter_text>plotting annual cycle<jupyter_code>g = sns.relplot(x = rois.index, y = 'avg_v',
                col = 'Year', row = 'distance_upG',
                #hue = 'daysep',
                data = rois,  facet_kws = {'sharex':'col',
                                           'sharey' : 'row'})

from itertools import product

pairs = list(product(pd.unique(rois['distance_upG']),
                     pd.unique(rois['Year'])))

for i, ax in enumerate(g.axes.flatten()):
    idx = (monthly['Year'] == pairs[i][1]) & (monthly['distance_upG'] == pairs[i][0])
    
    sns.lineplot(x = monthly[idx].index, y = 'v_monthly', color = 'k',
                 data = monthly[idx], ax = ax)
    
    for label in ax.get_xticklabels():
        label.set_rotation(45)
# g.savefig('../figures/isortuarsuupTasia_monthlyvelocities_wUnderlyingPoints_'+str(filterdays)+'dayFilter.png', bbox_inches='tight', dpi = 600)
NC = len(distanceUpG)
fig, axs = plt.subplots(ncols = NC, figsize = [20,6], sharex = True, sharey = True)

for i, ax in enumerate(axs):   
    
    if i == NC-1:
        lgnd = 'full'
    else:
        lgnd = False
    
    sns.lineplot(x = 'Month', y = 'v_monthly',
             hue = 'Year', palette = 'viridis',
             data = monthly[monthly['Region_num'] == i], ax = axs[i], legend = lgnd)
    
    sns.scatterplot(x = 'Month', y = 'v_monthly',
                 hue = 'Year', size = 'MCount', palette = 'viridis',
                 data = monthly[monthly['Region_num'] == i], ax = axs[i], legend = False)
    
    ax.set_xlim([1,12])
    ax.set_xticks(np.arange(1,12+1))
    ax.set_ylim(bottom = 0)
    ax.set_title(str(distanceUpG[i])+' m')
    
    
axs[0].set_ylabel('Velocity (m/yr)')
sns.move_legend(axs[NC-1], loc = 'upper right', bbox_to_anchor=(1.3, 1), facecolor = 'lightgrey')
plt.tight_layout()

hf.labelFig(fig);
# fig.savefig('../figures/landterminating/isortuarsuupTasia_annual_cycle_1500m_3000m_10000m_2013-21_'+str(filterdays)+'_dayfilter.png', bbox_inches='tight', dpi = 600)
## NORMALIZING THE ABOVE VELOCITIES

def nan_normalize(x):
    return (x - np.nanmean(x)) / np.nanstd(x)

# def normalize(x):
#     ind = np.isnan(x)
#     x_noNaN = x[~ind]
#     output = x
#     intermediate = (x_noNaN - (x_noNaN.sum(skipna = True) / x_noNaN.count())) / x_noNaN.std()
#     output[~ind] = intermediate
#     return output

monthly['annual_z'] = None
for yr in pd.unique(monthly['Year']):
    for R in pd.unique(monthly['Region_num']):
        idx = (monthly['Year'] == yr) & (monthly['Region_num'] == R)
        monthly.loc[idx, 'annual_z'] = nan_normalize(monthly.loc[idx, 'v_monthly'])

monthly.head()

cNorm = mpl.colors.Normalize(vmin = monthly['Year'].min(),
                             vmax = monthly['Year'].max())

fig, axs = plt.subplots(ncols = 2, nrows = 3, figsize = [24,16], sharex = True)

for i in range(3):
    
    if i == axs.shape[0] - 1:
        lgnd = 'full'
    else:
        lgnd = False
    
    # plot monthly (error-weighted) average velocities
    sns.lineplot(x = 'Month', y = 'v_monthly',
                 hue = 'Year', palette = 'viridis',
                 data = monthly[monthly['Region_num'] == i], ax = axs[i,0], legend = False).set(title='average monthly velocity\n'+str(distanceUpG[i])+' m')
    
    # error estimate
    # axs[i,0].errorbar(x = monthly[monthly['Region_num'] == i]['Month'],
    #              y = monthly[monthly['Region_num'] == i]['v_monthly'],
    #              yerr = monthly[monthly['Region_num'] == i]['verr_monthly'],
    #              fmt = 'None',
    #              color = 'lightgrey')
    
    # filling in ± 1 std
    for yr in pd.unique(monthly['Year']):
        if yr == 2017:
            lbl = '± 1 std'
        else: lbl = None
        idx = (monthly['Region_num'] == i) & (monthly['Year'] == yr)
        axs[i,0].fill_between(x = monthly.loc[idx]['Month'],
                              y1 = monthly.loc[idx,'v_monthly'] + monthly.loc[idx,'v_std'],
                              y2 = monthly.loc[idx,'v_monthly'] - monthly.loc[idx,'v_std'],
                              color = mpl.cm.viridis(cNorm(yr)),
                              alpha = 0.2,
                              label = lbl)
        axs[i,0].legend()
    
    # plot monthly z-score. each year has mean == 0 and std == 1. 
    # how this handles NaN - is a bit odd. missing months WILL mess with the z-score
    # unsure how best to handle this
    sns.lineplot(x = 'Month', y = 'annual_z',
                 hue = 'Year', palette = 'viridis',
                 data = monthly[monthly['Region_num'] == i], ax = axs[i,1], legend = lgnd).set(title = 'each year has mean = 0 ,std = 1')
    axs[i,1].axhline(0, c = 'k')
    
    # plot average z-score of all years for each month, in red
    monthly.groupby(['Region_num','Month'])['annual_z'].mean()[i].plot(ax = axs[i,1], c = 'r', linewidth = 5)


sns.move_legend(axs[axs.shape[0] - 1,1], loc = 'upper right', bbox_to_anchor=(1.19, 1), facecolor = 'grey')

for ax in axs[:,0]:
    ax.set_ylim([0,310])
    ax.set_xlim([1,12])
    
for ax in axs[:,1]:
    ax.set_ylim([-2.5,2.5])
    ax.set_xlim([1,12])
    
plt.tight_layout()
# fig.savefig('../figures/normalized_average_monthly_velocities_isortuarsuupTasia.png', bbox_inces = 'tight', dpi = 600)
from itertools import product

pairs = list(product(pd.unique(monthly['distance_upG']),
                     pd.unique(monthly['Year'])))

g = sns.catplot(x = 'Month', y = 'avg_v',
                col = 'Year', row = 'distance_upG',
                color = 'lightgrey', # alpha = 0.8,
                data = rois)

# the above creates categorical axes, whereby the xtick lables, do not necessarily equal the xtick values
# the three lines below, creates a new df that has the 'Month' value swapped for the appropriate label value
# this is used as input in the `lineplot` and ensures they line up.
labels = np.array([e.get_text() for e in g.axes.flat[-1].get_xticklabels()]).astype(int)
ticks = g.axes.flat[-1].get_xticks()
data = monthly.replace({'Month': dict(list(zip(labels,ticks)))})


for i, ax in enumerate(g.axes.flatten()):
    idx = (data['Year'] == pairs[i][1]) & (data['distance_upG'] == pairs[i][0])
    sns.lineplot(x = 'Month', y = 'v_monthly', color = 'k',
                 data = data.loc[idx], ax = ax, zorder = 2)
    
    ax.fill_between(x = data.loc[idx,'Month'],
                    y1 = data.loc[idx,'v_monthly'] + data.loc[idx, 'v_std'], # 'verr_monthly' or 'v_std'
                    y2 = data.loc[idx,'v_monthly'] - data.loc[idx, 'v_std'],
                    alpha = 0.8, color = 'cornflowerblue', label = 'std')
    
    ax.set_ylim([0,400])
    
# g.savefig('../figures/isortuarsuupTasia_monthlyvelocities_std_'+str(filterdays)+'dayFilter.png', bbox_inches='tight', dpi = 600)<jupyter_output><empty_output><jupyter_text>plotting (spatial) seasonal means<jupyter_code># Remove the annual, temporal mean to highlight the seasonal variability

# map_proj = ccrs.NorthPolarStereo(central_longitude=-45, true_scale_latitude=None, globe=None)
map_proj = ccrs.Stereographic(central_latitude=70.0, central_longitude=-45, false_easting=0.0, false_northing=0.0, true_scale_latitude=None, scale_factor=None, globe=None)

# plot difference between seasonal average and annual average
fig = (clipped.v.groupby('mid_date.season').mean() - clipped.v.mean(dim = 'mid_date')).plot(
    col = 'season',
    col_wrap = 2,
    robust = True,   # helpful for automatic colorscale that is robust to outliers - i.e. colormap can saturate
    subplot_kws={"projection": map_proj},
    cbar_kwargs={"label": "v anomaly (m/s)"},
    figsize = [10,6]
    # aspect = 1, size = 4
)

for ax in fig.axes.flatten():
    AOI.plot(ax = ax,
         facecolor = 'None',
         edgecolor = 'k')
    
plt.suptitle('seasonal cycle')

# hf.scale_bar(ax,length = 10)
# fig.savefig('../figures/isortuarsuupTasia_seasonalV_2013-2021')
# stream plot
# average vx, vy, and v by season
fig = clipped.groupby('mid_date.season').mean(skipna = True).sortby("y").plot.streamplot(x = "x", y = "y", u = "vx", v = "vy",
                                                                                         hue = "v",
                                                                                         robust = True,
                                                                                         col = "season",
                                                                                         col_wrap = 2,
                                                                                         subplot_kws={"projection": map_proj},
                                                                                         cbar_kwargs={"label": "v (m/s)"},
                                                                                         figsize = [10,6]
                                                                                        )
for ax in fig.axes.flatten():
    AOI.plot(ax = ax,
         facecolor = 'None',
         edgecolor = 'k')

# fig.savefig('../figures/isortuarsuupTasia_seasonalV_2013-21_avgStreamPlot')<jupyter_output><empty_output><jupyter_text>skewer cube at points(s) and export<jupyter_code>lakeX = -248478
lakeY = -2871545

# # create new data-array with at single point
lakecube = cube.sel(x=lakeX, y=lakeY,
              method='nearest') #,
#               # tolerance = 250).to_dataframe() 

lake = lakecube.to_dataframe()['v']
lake = pd.DataFrame(lake)
# lake['ns'] = lake.index.values.astype('float')
# sns.lmplot(data = lake, x = 'ns', y = 'v', scatter_kws = {'s' : 0.2})

landX = -240084
landY = -2838116

landcube = cube.sel(x=landX, y=landY,
                    method='nearest') #.to_dataframe()['v'] #,
land = landcube.to_dataframe()['v']
land = pd.DataFrame(land)
# land['ns'] = land.index.values.astype('float')
# sns.lmplot(data = land, x = 'ns', y = 'v', scatter_kws = {'s' : 0.2})

df = pd.DataFrame({'lake_v' : lake.v.values,
                   'land_v' : land.v.values},
                  index = lake.index)

df['date_dt'] = cube['date_dt'].values

# idx = df['date_dt'] <= pd.Timedelta('45D')
# df = df.loc[idx,:]

idx = df.index >= pd.Timestamp('2017-01-01')
df = df.loc[idx,:]
# sat = sat[idx]

df['yr'] = df.index.year + (df.index.dayofyear-1)/365
df.index = df['yr']

df = df.drop(columns = ['yr','date_dt'])

dflong = df.melt(ignore_index = False)
dflong.head()

idx = np.isnan(dflong['value'])
forLM = dflong.loc[~idx]
forLM.head()
forLM['value'] = forLM.value.astype('float64')
forLM['year'] = forLM.index.values

## get slopes
import scipy

# get lake slopes
lakeV = df.loc[~np.isnan(df['lake_v']),'lake_v']
LakeSlope, LakeIntercept, r, p, se = scipy.stats.linregress(lakeV.index, lakeV.values)

print(LakeSlope)
print(r**2)
print(p)

landV = df.loc[~np.isnan(df['land_v']), 'land_v']

landSlope, landintercept, r, p, se = scipy.stats.linregress(landV.index, landV.values)

print(landSlope)
print(r**2)
print(p)
print(se)

text = 'Rate of change 2017-2021\nLake: ' + str(np.round(LakeSlope,2)) + ' m/yr$^{2}$\nLand: ' + str(np.round(landSlope,2)) + ' m/yr$^{2}$'
# print(text)<jupyter_output><empty_output><jupyter_text>uncertainity workings out<jupyter_code>idx = 11840
v = clipped.isel(mid_date = idx)['v']
v_error = clipped.isel(mid_date = idx)['v_error']
vx = clipped.isel(mid_date = idx)['vx']
vy = clipped.isel(mid_date = idx)['vy']
vy_error = clipped.isel(mid_date = idx)['vy_error']
vx_error = clipped.isel(mid_date = idx)['vx_error']

my_v, my_v_error = hf.v_uncert(vx, vy, vx_error.values, vy_error.values)

differ = (my_v.values - v.values).flatten()
differ = differ[~np.isnan(differ)]

fig, axs = plt.subplots(ncols = 3, figsize = [21,6])

axs[0].scatter(my_v, v)
axs[0].set_xlabel('calculated_v (m/yr)')
axs[0].set_ylabel('itslive_v (m/yr)')

axs[1].hist(differ)
axs[1].set_title('distribution of differences\ncalculated_v - itslive_v')
axs[1].set_ylabel('count')
axs[1].set_xlabel('velocity difference (m/yr)')

axs[2].scatter(my_v_error, v_error)
axs[2].set_xlabel('calculated_v_error (m/yr)')
axs[2].set_ylabel('itslive_v_error (m/yr)')

for i in [0,2]:
    hf.xEqualsy(axs[i])
    
axs[0].set_xticks([0,200,400,600]);
axs[0].set_yticks([0,200,400,600]);

axs[1].set_xticks([-0.5,-0.25,0,0.25,0.5]);

axs[2].set_xticks([58,60,62,64,66,68]);
axs[2].set_yticks([58,60,62,64,66,68]);

hf.labelFig(fig, pos = 'topleft')
# fig.savefig('../figures/uncertainty_calculations.png', bbox_inches = 'tight', dpi = 300)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15566717361722565, 'fraction_numerical': 0.019671015768830544, 'mean_word_length': 2.4846245031727214, 'pattern_counts': {'":': 4, '<': 56, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 110, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '68055', 'license': 'MIT', 'n_tokens_mistral': 17922, 'n_tokens_neox': 16741, 'n_words': 3689, 'path': 'notebooks/cube_plotting.ipynb', 'repo_name': 'tlohde/its_live', 'size': '49972'}