<jupyter_start><jupyter_text>okay - running the code with the low res data works; results are all 0 (zeros) for high res data though...so now need to compare the formats of the datasets.  ncdump in terminal~~~~~~~~~~~~~~~~~~~~~~~~~high res~~~~~~~~~~~~~~~~~~~~~~~~~HR data - 188 profiles557 N_LEVELS = depths?46 parametersmissing data recorded as -1e10~~~~~~~~~~~~~~~~~~~~~~~~~low res~~~~~~~~~~~~~~~~~~~~~~~~~LR data - 188 profiles67 N_LEVELS = depths?46 parametersmissing data recorded as -1e10N_LEVELS DIFFER - high res has much greater pressure/depth resolution than low res(Nitrate is low res anyway) see below:<jupyter_code># fd
print('High Res')
HRfd.Pressure.isel(N_PROF=0).values
print('Low Res')
LRfd.Pressure.isel(N_PROF=0).values
LRfd.Pressure.isel(N_PROF=0).shape
help(np.convolve)
zi = np.arange(0, 1600, 5)
print(zi)
zi.shape
np.interp?
print(len(LRfd.Nitrate))
np.asarray?
print(LRfd.Pressure[:,::-1])
#help(np.interp)
#help(np.asarray)
#interpolate the data onto the standard depth grid given by x_int
def interpolate(x_int, xvals, yvals):
    yvals_int = []
    for n in range(0, len(yvals)): #len(yvals) being the number of profiles
        yvals_int.append(np.interp(x_int, xvals[n, :], yvals[n, :]))
    #convert the interpolated data from a list to numpy array
    return np.asarray(yvals_int)
abs.argmin?
#.argmin
print(LRfd.Nitrate.shape)
print(len(LRfd.Nitrate))
LRfd.Nitrate.isel(N_PROF=0).values # show values of nitrate for first profile - in rev because deepest measurement first
# Interpolate nitrate and poc
zi = np.arange(0, 1600, 5)
nitr_int = interpolate(zi, LRfd.Pressure[:, ::-1], LRfd.Nitrate[:, ::-1])

print(nitr_int[0]) # values for profile 1
print(len(zi))
print(len(nitr_int))

print(zi.shape)
print(nitr_int.shape)
#help(np.array)
# calculate the vertically integrated data column inventory using the composite trapezoidal rule
def integrate(zi, data, depth_range):
    n_profs = len(data)
    zi_start = abs(zi - depth_range[0]).argmin()
    zi_end =  abs(zi - depth_range[1]).argmin()
    zi_struct = np.ones((n_profs, 1)) * zi[zi_start : zi_end]
    data = data[:, zi_start : zi_end]
    col_inv = []
    for n in range(0, len(data)):
        col_inv.append(np.trapz(data[n,:][~np.isnan(data[n,:])], zi_struct[n,:][~np.isnan(data[n,:])]))
    return col_inv
depth_range=[25,200]
zi_start = abs(zi - depth_range[0]).argmin()
print(zi)
print(depth_range[0])
print(zi_start)
depth_range=[25,200]
zi_end = abs(zi - depth_range[1]).argmin()
print(zi)
print(depth_range[1])
print(zi_end)
depth_range=[25,200]
zi_start = abs(zi - depth_range[0]).argmin()
print(depth_range[0])
print(zi_start)
zi_end = abs(zi - depth_range[1]).argmin()
print(depth_range[1])
print(zi_end)

n_profs = len(nitr_int)
zi_struct = np.ones((n_profs, 1)) * zi[zi_start : zi_end+1] # could add +1 to get all the way to 200?
print(zi_struct)
zi_struct.shape
#np.trapz?
np.where?
zi = np.arange(0, 1600, 2)
print(zi)
# Integrate nitrate and poc
nitr = np.array(integrate(zi, nitr_int, [25, 200]))
print(nitr.shape)
print(nitr.size)
print(nitr)
help(abs)

def find_extrema(data, date_range, find_func):
    # Find indices of float profiles in the date range
    date_mask = (JULDLR > date_range[0]) & (JULDLR < date_range[1])
    
    # Get the index where the data is closest to the find_func
    index = np.where(data[date_mask] == find_func(data[date_mask]))[0][0]
    
    # Get the average data for the month of the extrema
    month_start = JULDLR[date_mask][index].replace(day = 1)
    month_dates = (JULDLR > month_start) & (JULDLR < month_start + pd.Timedelta(days = 30))
    month_avg = np.mean(data[date_mask])
    
    return month_avg, JULDLR[date_mask][index], data[date_mask][index]

years = [2015, 2016, 2017, 2018]
nitr_extrema = []
nitr_ancp = []
for y in years:
    winter_range = [pd.datetime(y, 8, 1), pd.datetime(y, 12, 1)]
    summer_range = [pd.datetime(y, 12, 1), pd.datetime(y + 1, 4, 1)]

    # Find maximum winter and minimum summer nitrate
    avg_max_nitr, max_nitr_date, max_nitr = find_extrema(nitr, winter_range, np.max)
    avg_min_nitr, min_nitr_date, min_nitr = find_extrema(nitr, summer_range, np.min)
    
    # Convert to annual nitrate drawdown
    redfield_ratio = 106.0/16.0 #106C:16NO3-
    
    # Nitrate units: umol/kg
    
    nitr_drawdown = (avg_max_nitr - avg_min_nitr)/1000.0 * redfield_ratio
    
    nitr_ancp.append(nitr_drawdown)
    
    nitr_extrema.append(((max_nitr, max_nitr_date), (min_nitr, min_nitr_date)))
    
    print(y, max_nitr_date, max_nitr, avg_max_nitr)
    print(y, min_nitr_date, min_nitr, avg_min_nitr)
    
print('And ANCP are:', years)
print(nitr_ancp)
JULDLR.month.shape
winter_range = [pd.datetime(2015, 8, 1), pd.datetime(2015, 12, 1)]
date_mask = (JULDLR > winter_range[0]) & (JULDLR < winter_range[1])
index = np.where(nitr[date_mask] == np.max(nitr[date_mask]))[0][0]
print(index)
month_start2 = JULDLR[date_mask][index]
month_start = JULDLR[date_mask][index].replace(day = 1)
print(month_start2)
print(month_start)
month_dates = (JULDLR > month_start) & (JULDLR < month_start + pd.Timedelta(days = 30)) 
print(month_dates)
month_dates.shape
#pd.Timedelta?
month_avg = np.mean(nitr[date_mask])
print(month_avg)
print(month_start.month)
print(month_start)
# for 2015: winter max in month 9
    # index - number 5 within datemask

month_avg = np.mean(nitr[date_mask])     
print(month_avg)
#month_ave = np.mean(nitr[date_mask] where month = month_start.month)

print(date_mask)
month_mask = (JULDLR.month == month_start.month) & (JULDLR.year == month_start.year)

month_mask2 = (JULDLR.month[date_mask] == month_start.month)

#month_mask = (JULDLR.month == month_start.month)
#month_mask = (JULDLR == month_start.year)
print(month_mask)
print(month_mask2)

#print(month_start.year)
month_avg = np.mean(nitr[month_mask])
print(month_avg)
help(plt.axvspan)
#try to integrate using high res data
print(HRfd.Pressure.shape)
print(HRfd.Pressure.size)
print('High Res')
HRfd.Pressure.isel(N_PROF=1).values
# Integrate nitrate and poc
nitr = np.array(integrate(zi, nitr_int, [25, 200])) # integrate interpolated nitrate values between 25m-200m



# my 'years' range
timeseries = np.arange(2002, 2021, 1)
print(timeseries)
#HRfd.Lat.isel(N_PROF=0).values #for just one profile
HRfd.Lat.values
zi = np.arange(0, 1600, 5) # 5 = 320 depth intervals between 0m to 1595m
nitr_int = interpolate(zi, HRfd.Pressure[:, ::-1], HRfd.Nitrate[:, ::-1]) # interpolate nitrate values across zi depth intervals for all 188 profiles

# Integrate nitrate and poc - total nitrate in upper 200m
upperlim=25
lowerlim=200
nitr = np.array(integrate(zi, nitr_int, [upperlim, lowerlim])) # integrate interpolated nitrate values between 25m-200m 

print(nitr_int[0])
#HRfd.nitr_int.isel(N_PROF=0).values
HRfd.Pressure[:, ::-1].values
zi = np.arange(0, 1600, 5) # 5 = 320 depth intervals between 0m to 1595m
nitr_int = interpolate(zi, LRfd.Pressure[:, :], LRfd.Nitrate[:, :]) # interpolate nitrate values across zi depth intervals for all 188 profiles

# Integrate nitrate and poc - total nitrate in upper 200m
upperlim=25
lowerlim=200
nitr = np.array(integrate(zi, nitr_int, [upperlim, lowerlim])) # integrate interpolated nitrate values between 25m-200m 

print(nitr_int[10,9])

#print(nitr)
print(nitr_int.shape)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1456164383561644, 'fraction_numerical': 0.03575342465753425, 'mean_word_length': 5.176818950930626, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, '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': '240167', 'license': 'MIT', 'n_tokens_mistral': 3017, 'n_tokens_neox': 2789, 'n_words': 648, 'path': 'notebooks/nitrate_drawdown_workingitout.ipynb', 'repo_name': 'ClaraCDouglas/net_community_production', 'size': '7300'}