<jupyter_start><jupyter_text>Setup<jupyter_code># import packages
%run ../../global_packages.py

# get the global parameters
%run ../../global_pars.py

# import your local functions
sys.path.insert(1, '../../')
from global_functions import *

# make sure the figures plot inline rather than at the end
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Parameters and Paths<jupyter_code>outfn = 'sla_processed.nc'

# get data from tigress
inpath = 'raw_data/concatenated/'
infn = 'all_dt_global_twosat_phy_l4__vDT2018_19930101_to_20200603.nc'<jupyter_output><empty_output><jupyter_text>Get Data and Subset<jupyter_code>ds = xr.open_dataset(inpath + infn)

# Subset ------------------------------------------------#
lat_slice = slice(lat_bounds[0], lat_bounds[1])
lon_slice = slice(lon_bounds[0], lon_bounds[1])
time_slice = slice(ts,te)

ds = ds.sel(lat=lat_slice,lon=lon_slice, time = time_slice)

# Get daily sla
daily_sla = ds['sla']
lat = daily_sla.lat.values
lon = daily_sla.lon.values
ds<jupyter_output><empty_output><jupyter_text>Detrend SLAResources: [link](https://ecco-v4-python-tutorial.readthedocs.io/ECCO_v4_Example_calculations_with_scalar_quantities.html), [link](https://groups.google.com/g/xarray/c/fz7HHgpgwk0), [link](https://stackoverflow.com/questions/28647172/numpy-polyfit-doesnt-handle-nan-values), [link](https://stackoverflow.com/questions/28647172/numpy-polyfit-doesnt-handle-nan-values), [link](https://stackoverflow.com/questions/17638137/curve-fitting-to-a-time-series-in-the-format-datetime)<jupyter_code>%%time

# stack lat and lon into a single dimension called allpoints
stacked = daily_sla.stack(allpoints=['lat','lon'])

# set places where there are nans to zero since polyfit can't deal with them
stacked_nonan = stacked.fillna(0)

# convert date to a number to polyfit can handle it
datenum = dates.date2num(stacked_nonan.time)
daily_sla_slope, daily_sla_intercept = np.polyfit(datenum, stacked_nonan, 1)

#reshape the data
daily_sla_slope = np.reshape(daily_sla_slope, daily_sla.shape[1:3])
daily_sla_intercept = np.reshape(daily_sla_intercept, daily_sla.shape[1:3])

# define a function to compute a linear trend of a timeseries
def linear_detrend(y):
    x = dates.date2num(y.time)
    m, b = np.polyfit(x, y, 1)
    # we need to return a dataarray or else xarray's groupby won't be happy
    return xr.DataArray(y - (m*x + b))

# apply the function over allpoints to calculate the trend at each point
daily_sla_dtrnd = stacked_nonan.groupby('allpoints').apply(linear_detrend)
# unstack back to lat lon coordinates
daily_sla_dtrnd = daily_sla_dtrnd.unstack('allpoints')

# fill all points we set originally to zero back to nan
daily_sla_dtrnd = daily_sla_dtrnd.where(~np.isnan(daily_sla))
# find global means for comparison
stacked = daily_sla.stack(allpoints=['lat','lon'])
global_mean = stacked.mean(dim='allpoints',skipna=True)
# global_mean.plot()

# find global means for comparison 
stacked = daily_sla_dtrnd.stack(allpoints=['lat','lon'])
global_mean_dtrnd = stacked.mean(dim='allpoints',skipna=True)
# global_mean_dtrnd.plot()
# plt.legend


fig = plt.figure(figsize=(6,3),dpi=200)
plt.plot(global_mean.time,global_mean)
plt.plot(global_mean_dtrnd.time,global_mean_dtrnd, label = 'Detrended')
plt.legend()
plt.title('Indian Ocean Sea Level Anomaly')
plt.ylabel('Sea Level Anomaly ($m$)')
plt.xlabel('Time')<jupyter_output>/home/jennap/anaconda3/lib/python3.7/site-packages/pandas/plotting/_converter.py:129: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.

To register the converters:
	>>> from pandas.plotting import register_matplotlib_converters
	>>> register_matplotlib_converters()
  warnings.warn(msg, FutureWarning)<jupyter_text>Downsample to Monthly and Seasonal Temporal Resolution<jupyter_code># monthly
mon_sla_dtrnd = daily_sla_dtrnd.resample(time='1MS').mean(dim="time")<jupyter_output>/home/jennap/anaconda3/lib/python3.7/site-packages/xarray/core/nanops.py:161: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Find Monthly Climatology and AnomalyResources: [link](http://xarray.pydata.org/en/stable/examples/monthly-means.html)<jupyter_code># monthly clim
mon_sla_mon_clim_dtrnd = daily_sla_dtrnd.groupby('time.month').mean('time') 

# monthly anom = monthly avg data - monthly climatology
mon_sla_mon_anom_dtrnd = mon_sla_dtrnd.groupby('time.month') - mon_sla_mon_clim_dtrnd<jupyter_output><empty_output><jupyter_text>Find Coastlines<jupyter_code># ------------------------------------------------------------#
# block out a few places that you don't want to include
# ------------------------------------------------------------#

mask = mon_sla_dtrnd[0,:,:]
lat = daily_sla.lat.values
lon = daily_sla.lon.values
xx,yy = np.meshgrid(lon,lat)
# set Gulfs to Nans
mask = mask.where(~((yy>22) & (xx<57.5)))
mask = mask.where(~((yy>10) & (xx<51)))

# AS inlets
mask = mask.where(~((yy>21.5) & (yy<23) & (xx>72) & (xx<73)))
mask = mask.where(~((yy>22.25) & (yy<23) & (xx>69) & (xx<70.5)))
mask = mask.where(~((yy>23.25) & (yy<24) & (xx>68.5) & (xx<70)))

# sri lanka are
mask = mask.where(~((yy>8) & (yy<10.5) & (xx>78) & (xx<80)))

# Sumatra area
mask = mask.where(~((yy>4.5) & (yy<5.5) & (xx>95.5) & (xx<100)))
mask = mask.where(~((yy>3) & (yy<5.5) & (xx>97.5) & (xx<100)))
mask = mask.where(~((yy>0) & (yy<10) & (xx>99) & (xx<100)))
mask = mask.where(~((yy>0) & (yy<2.5) & (xx>97) & (xx<98)),1)

# equator
mask = mask.where(~(yy<0))
# ------------------------------------------------------------#
# locate points along the BoB and AS
# ------------------------------------------------------------#

def find_coast(arr):
    import itertools
    
    # create empty lists to add row and col info to
    rowind = []
    colind = []

    # create an array of indices
    rows = np.arange(arr.shape[0])
    cols = np.arange(arr.shape[1])

    # find if the sum of a block around a point is a nan (meaning one of the values at least must have been a nan)
    for row,col in itertools.product(rows,cols):
#         cond = (~np.isnan(arr[row,col])) & (np.isnan(np.sum(arr[max(0,row-1):min(arr.shape[0],row+2),max(0,col-1):min(arr.shape[1],col+2)])))
        rowcond = (np.isnan(np.sum(arr[max(0,row-1):min(arr.shape[0],row+2),col])))
        colcond = (np.isnan(np.sum(arr[row,max(0,col-1):min(arr.shape[1],col+2)])))

#         if  (~np.isnan(arr[row,col])) & cond):
        if  (~np.isnan(arr[row,col])) & (rowcond | colcond):
            rowind.append(rows[row].tolist())
            colind.append(cols[col].tolist())
    
    return np.array(rowind), np.array(colind)

sta_latinds_unord, sta_loninds_unord = find_coast(np.array(mask))
# ------------------------------------------------------------#
# Remove Manually some other points
# ------------------------------------------------------------#

# sta_latinds_unord = sta_latinds_unord[46:]
# sta_loninds_unord = sta_loninds_unord[46:]

sta_latinds_unord = sta_latinds_unord[125:]
sta_loninds_unord = sta_loninds_unord[125:]

# AS
lons_from_sta = lon[sta_loninds_unord]
lats_from_sta = lat[sta_latinds_unord]

ind = ~((lons_from_sta <52) & (lats_from_sta < 15))
sta_latinds_unord = sta_latinds_unord[ind]
sta_loninds_unord = sta_loninds_unord[ind]

#BoB
lons_from_sta = lon[sta_loninds_unord]
lats_from_sta = lat[sta_latinds_unord]

ind = ~(lons_from_sta >99)
sta_latinds_unord = sta_latinds_unord[ind]
sta_loninds_unord = sta_loninds_unord[ind]

# ------------------------------------------------------------#
# Plot
# ------------------------------------------------------------#

bounds = [40,100,-5,30]

def add_land(ax,bounds):
    res = '50m'
    from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    land = cfeature.NaturalEarthFeature('physical', 'land', res,
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])
    ax.add_feature(land,color='gray',zorder = 0)
    ax.background_patch.set_facecolor('k')
    ax.coastlines(resolution=res,zorder = 0)
    g = ax.gridlines(draw_labels=True,alpha=0)
    g.xlabels_top = False
    g.ylabels_right = False
    g.xlabel_style = {'size': 15}
    g.ylabel_style = {'size': 15}
    g.xformatter = LONGITUDE_FORMATTER
    g.yformatter = LATITUDE_FORMATTER
    ax.axes.axis('tight')
    ax.set_extent(bounds, crs=ccrs.PlateCarree())
    
cmin = -0.2
cmax = 0.2
levels = np.linspace(cmin, cmax, 10)

# Start figure
fig = plt.figure(figsize=(16, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
add_land(ax,bounds)
# cf = ax.contourf(lon, lat,sla_mon_anom[0,:,:],levels = levels,cmap=plt.cm.PuOr_r, extend="both")
plt.scatter(lon[sta_loninds_unord],lat[sta_latinds_unord],
            s = 160,c=np.arange(len(sta_latinds_unord)),marker = '.',cmap =plt.cm.PiYG)
plt.colorbar(label = 'Station No.')
plt.title('Unordered Stations')
# ------------------------------------------------------------#
# put in the right order
# ------------------------------------------------------------#


def order_coast(loninds,latinds,sta_zero):
    # find based on radius
    zipped_lists = zip(loninds, latinds)
    sorted_pairs = sorted(zipped_lists, reverse=True)

    # sort by lon
    tuples = zip(*sorted_pairs)
    lon_list,lat_list  = [ list(tuple) for tuple in  tuples]

    pos = []
    curr_sta = []
    rem_sta = sorted_pairs
    for i in range(len(lon_list)):
#         print(i)
        if i == 0:
            curr_sta.append(sta_zero)
            rem_sta.remove(sta_zero) 
            prev_sta = sta_zero
        else:
            prev_sta = curr_sta[i-1]

        start_len = len(curr_sta)
        for j,(lo, la) in enumerate(rem_sta):
            next_sta = (rem_sta[j])
            
            diff = tuple(map(lambda l, k: l - k, curr_sta[i], next_sta))
            
            # check uplr first
            if(next_sta != prev_sta) & (all(np.abs(diff) == [0,1])) | (all(np.abs(diff) == [1,0])):
                curr_sta.append(next_sta)
                rem_sta.remove(next_sta) 
                break

            # then check diagonals
            elif (next_sta != prev_sta) & (all(np.abs(diff) == [1,1])):
                curr_sta.append(next_sta)
                rem_sta.remove(next_sta) 
                break
                
        if len(curr_sta) == start_len:
            print('No Next Station Found. Returning Previous Stations Only.')
            print(curr_sta[i],rem_sta)
            break
        
    sta_lonind, sta_latind  = map(np.array, zip(*curr_sta))
    
    return sta_lonind, sta_latind

#initial station
sta_zero = (sta_loninds_unord[0],sta_latinds_unord[0]) # 98.625, 10.375

#calculate
sta_loninds, sta_latinds =order_coast(sta_loninds_unord,sta_latinds_unord,sta_zero)

# remove the one weird nan - use next code block to help identify it
sta_loninds = np.delete(sta_loninds,85)
sta_latinds = np.delete(sta_latinds,85)

# cut off extra arabian sea
sta_latinds = sta_latinds
sta_loninds = sta_loninds


# ------------------------------------------------------------#
# identify first stations for EQ, BOB, and AS
# ------------------------------------------------------------#
# old country based ones
# JA = 0
# sBoB = 120
# mBoB = 210
# SL = 300
# mAS = 410
# eAS = 494
# loc_list = [JA,sBoB,mBoB,SL,mAS,eAS]

JA = 0
sEBoB = 75
sWBoB = 225
sEAS = 294
sWAS = 415
eWAS = 504 # 494

loc_list = [JA,sEBoB,sWBoB,sEAS,sWAS,eWAS]


# ------------------------------------------------------------#
# plot
# ------------------------------------------------------------#
bounds = [40,100,-5,30]

def add_land(ax,bounds):
    res = '50m'
    from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    land = cfeature.NaturalEarthFeature('physical', 'land', res,
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])
    ax.add_feature(land,color='gray',zorder = 0)
    ax.background_patch.set_facecolor('k')
    ax.coastlines(resolution=res,zorder = 0)
    g = ax.gridlines(draw_labels=True,alpha=0)
    g.xlabels_top = False
    g.ylabels_right = False
    g.xlabel_style = {'size': 15}
    g.ylabel_style = {'size': 15}
    g.xformatter = LONGITUDE_FORMATTER
    g.yformatter = LATITUDE_FORMATTER
    ax.axes.axis('tight')
    ax.set_extent(bounds, crs=ccrs.PlateCarree())
    
cmin = -0.2
cmax = 0.2
levels = np.linspace(cmin, cmax, 10)

# Start figure
fig = plt.figure(figsize=(16, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
add_land(ax,bounds)

p =plt.scatter(lon[sta_loninds],lat[sta_latinds],
            s = 100,c=np.arange(len(sta_latinds)),marker = '.',cmap =plt.cm.PiYG) 
plt.scatter(lon[sta_loninds[loc_list]],lat[sta_latinds[loc_list]],
            s = 200,c='darkblue',marker = '.')
plt.colorbar(p,label = 'Station No.')
plt.title('Ordered Stations')

print('lons: ',lon[sta_loninds[loc_list]])
print('lats: ',lat[sta_latinds[loc_list]])
# use this to find the random nan value that appeasr someotimes

data = np.zeros([12,sta_loninds.shape[0]])

for jj in range(sta_loninds.shape[0]):
        data[:,jj] = mon_sla_dtrnd[0,sta_latinds[jj],sta_loninds[jj]]
        
np.argwhere(np.isnan(data))<jupyter_output><empty_output><jupyter_text>Save Data<jupyter_code># convert to xarray dataset
ds_out = xr.Dataset(coords={'lon': mon_sla_dtrnd.lon,
                    'lat': mon_sla_dtrnd.lat,
                    'time': mon_sla_dtrnd.time})

# add variables to dataset
ds_out["daily_sla"]=xr.DataArray(daily_sla_dtrnd,dims = ['day_time','lat', 'lon'],
                     coords =[daily_sla_dtrnd.time,daily_sla_dtrnd.lat,daily_sla_dtrnd.lon])

ds_out["mon_sla"]=xr.DataArray(mon_sla_dtrnd,dims = ['time','lat', 'lon'],
                     coords =[mon_sla_dtrnd.time,daily_sla_dtrnd.lat,daily_sla_dtrnd.lon])

# clim
ds_out["mon_sla_mon_clim"]=xr.DataArray(mon_sla_mon_clim_dtrnd,dims = ['month','lat', 'lon'],
                     coords =[mon_sla_mon_clim_dtrnd.month,daily_sla_dtrnd.lat,daily_sla_dtrnd.lon])

# anom
ds_out["mon_sla_mon_anom"]=xr.DataArray(mon_sla_mon_anom_dtrnd,dims = ['time','lat', 'lon'],
                   coords =[mon_sla_mon_anom_dtrnd.time,daily_sla_dtrnd.lat,daily_sla_dtrnd.lon])

no_stations =np.linspace(1,len(sta_loninds),len(sta_loninds))

ds_out['sta_loninds'] = xr.DataArray(sta_loninds,dims = ['no_stations'],coords =[no_stations])
ds_out['sta_latinds'] = xr.DataArray(sta_latinds,dims = ['no_stations'],coords =[no_stations])
ds_out['loc_list'] = xr.DataArray(loc_list,dims = ['loc_list'],coords =[loc_list])

ds_out['sta_lon'] = xr.DataArray(lon[sta_loninds],dims = ['no_stations'],coords =[no_stations])
ds_out['sta_lat'] = xr.DataArray(lat[sta_latinds],dims = ['no_stations'],coords =[no_stations])
ds_out['loc_list'] = xr.DataArray(loc_list,dims = ['loc_list'],coords =[loc_list])



# delete if already present
if os.path.isfile(outfn):
    os.remove(outfn)

ds_out.to_netcdf(outfn,mode='w',format = "NETCDF4")
ds_out<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1949286846275753, 'fraction_numerical': 0.02159270998415214, 'mean_word_length': 3.1987801497088992, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 56, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '143971', 'license': 'MIT', 'n_tokens_mistral': 5938, 'n_tokens_neox': 5647, 'n_words': 950, 'path': 'data_processing/2_SLA/0_process_sla.ipynb', 'repo_name': 'jlpearso/coastlines-at-risk-of-hypoxia-in-the-northern-indian-ocean', 'size': '15144'}