Case study: Citizen science phosphate data comparisons

# match samples that occur within 24 hours of one another
merged_timeseries_eawq_ti = pd.DataFrame()
# process one site at a time
for idx in merged_sites_wq_ti.index:
    
    # download the thames initiative time series data
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_sites_wq_ti.loc[idx,'properties.id_right']}&datasetID=thames_initiative&determinand=Soluble%20reactive%20phosphorus%20(%CE%BCg/L)")
    ti_timeseries=pd.json_normalize(response.json(), record_path = 'data')

    #format the timeseries as a dataframe
    if len(ti_timeseries) > 0:
        ti_timeseries['dti'] = pd.DatetimeIndex(ti_timeseries['datetime'])
        ti_timeseries.loc[ti_timeseries['value'] == '','value'] = np.nan #set missing values to nan
        ti_timeseries.loc[ti_timeseries['value'] == 'No Sample','value'] = np.nan #set missing values to nan
        ti_timeseries['value'] = ti_timeseries['value'].astype('float')
        ti_timeseries['value'] = ti_timeseries['value']/1000 #ug to mg
        ti_timeseries['date'] = ti_timeseries['dti'].dt.date
        ti_timeseries['jul2'] = pd.DatetimeIndex(ti_timeseries['dti']).to_julian_date()
        ti_timeseries['site'] = merged_sites_wq_ti.loc[idx,'properties.id_right']
        ti_timeseries['distance'] = merged_sites_wq_ti.loc[idx,'distance']

        # download the ea wq time series data
        response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_sites_wq_ti.loc[idx,'properties.id_left']}&datasetID=ea_water_quality&determinand=0180")
        eawq_timeseries=pd.json_normalize(response.json(), record_path = 'data')

        # format the timeseries as a dataframe
        if len(eawq_timeseries)>0:
            eawq_timeseries['dti'] = pd.DatetimeIndex(eawq_timeseries['datetime'])
            eawq_timeseries['date'] = eawq_timeseries['dti'].dt.date
            eawq_timeseries['jul2'] = pd.DatetimeIndex(eawq_timeseries['dti']).to_julian_date()
            eawq_timeseries['site'] = merged_sites_wq_ti.loc[idx,'properties.id_left']
            # merge the ea wq and ti timeseries based on sampling time
            # tolerance 1 = 1 day = 24 hours
            merged_timeseries_eawq_ti = pd.concat([merged_timeseries_eawq_ti,pd.merge_asof(ti_timeseries.sort_values('jul2'),eawq_timeseries.sort_values('jul2'), left_on='jul2', right_on='jul2', direction='nearest', tolerance=1)])
            
# calculate how far apart the samples were in time
merged_timeseries_eawq_ti['timedelta'] = merged_timeseries_eawq_ti['dti_x'] - merged_timeseries_eawq_ti['dti_y']
merged_timeseries_eawq_ti=merged_timeseries_eawq_ti.filter(['value_x','value_y','dti_x','dti_y','site_x','site_y', 'distance', 'timedelta'])
merged_timeseries_eawq_ti=merged_timeseries_eawq_ti.rename(columns={'value_x':'ti_value', 'value_y':'eawq_value','dti_x':'ti_dti','dti_y':'eawq_dti','site_x':'ti_site','site_y':'eawq_site'})
merged_timeseries_eawq_ti = merged_timeseries_eawq_ti[merged_timeseries_eawq_ti['eawq_value'].notnull()]
merged_timeseries_eawq_ti['hours_diff'] = abs(merged_timeseries_eawq_ti['timedelta'].dt.total_seconds()/3600)

fig, axes=plt.subplots(nrows=2)
fig.set_size_inches(10,10)
#scatter plot
sns.scatterplot(x='ti_value',y='eawq_value',data=merged_timeseries_eawq_ti, ax=axes[0])
axes[0].set_ylabel('EA phosphate concentration (mg/L)')
axes[0].set_xlabel('TI soluble reactive phosphorus concentration (mg/L)')
#ax.set_ylim(0,0.4)
x0, x1 = axes[0].get_xlim()
y0, y1 = axes[0].get_ylim()
lims = [max(x0, y0), min(x1, y1)]
axes[0].plot(lims, lims, '-r')
r, p = pearsonr(merged_timeseries_eawq_ti['ti_value'], merged_timeseries_eawq_ti['eawq_value'])
r2 = r**2
axes[0].text(x=0,y=1.5,s=r"$r^{2}$"+f"={np.round(r2,3)}",c='red')
axes[0].text(x=-.1,y=2.4,s="A")
# make a timeline
# Sort by ti_dti
df =merged_timeseries_eawq_ti.sort_values(by='ti_dti')
# stagger the stems 
levels = [1, 2, 3, -1, -2, -3] * len(df)
levels = levels[:len(df)]
# Draw horizontal line
axes[1].axhline(0, color='black', linewidth=1)
# Draw stems and points
sns.scatterplot(ax=axes[1],x='ti_dti', y=[0]*len(df), color='blue', s=50, zorder=3,data=df, legend=False, edgecolor='k')
ci = 0 
for i in df.index:
    if pd.notna(df.loc[i,'ti_dti']):
        axes[1].vlines(df.loc[i,'ti_dti'], 0, levels[ci], linewidth=1)
    ci += 1 
# Ensure annotations aren't too long
df['annot'] = df["ti_site"].str.slice(0,10)
# Annotate each point
for d, l, sname, di, hd in zip(df['ti_dti'], levels, df['annot'],df['distance'],df['hours_diff']):
        di = np.round(di,2)
        hd = np.round(hd,2)
        axes[1].text(d, l, f"{sname}", ha='center', va='bottom' if l > 0 else 'top', fontsize=8, rotation=30)
# Format x-axis
axes[1].xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
# Remove y-axis 
axes[1].get_yaxis().set_visible(False)
#remove some spines
axes[1].spines[["left", "top", "right"]].set_visible(False)
axes[1].set_ylim([-7,5])
axes[1].set_xticks(ticks=axes[1].get_xticks(),rotation=45,labels=axes[1].get_xticklabels())
axes[1].set_xlabel("")
axes[1].text(x=axes[1].get_xticks()[0],y=-5,s=f"μ distance (m): { np.round(df['distance'].mean(),2)} \nμ time difference (hrs): { np.round(df['hours_diff'].mean(),2)}")
axes[1].text(x=axes[1].get_xticks()[0],y=5,s=f"    B")

fig.text(.5, -.05, "Figure 1. A: A scatter plot comparison of EA and TI phosphate measurements. \n The red line is equivalence and the r^2 value represents the correlation between the two measures. \n B: A timeline showing when and where the samples were matched, labelled using the site name and date corresponding to the TI data in each pair. \n The average distance and time between the samples is given in text.", ha='center')

merged_timeseries_fww_ti = pd.DataFrame()
for idx in merged_sites_fww_ti.index:
    
    # download the thames initiative (ti) timeseries data
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_sites_fww_ti.loc[idx,'properties.id_right']}&datasetID=thames_initiative&determinand=Soluble%20reactive%20phosphorus%20(%CE%BCg/L)")
    ti_timeseries=pd.json_normalize(response.json(), record_path = 'data')

    # format the ti data as a pandas dataframe
    if len(ti_timeseries) > 0:
        ti_timeseries['dti'] = pd.DatetimeIndex(ti_timeseries['datetime'])
        ti_timeseries.loc[ti_timeseries['value'] == '','value'] = np.nan #set missing values to nan
        ti_timeseries.loc[ti_timeseries['value'] == 'No Sample','value'] = np.nan #set missing values to nan
        ti_timeseries['value'] = ti_timeseries['value'].astype('float')
        ti_timeseries['value'] = ti_timeseries['value']/1000 #ug to mg
        ti_timeseries['date'] = ti_timeseries['dti'].dt.date
        ti_timeseries['jul2'] = pd.DatetimeIndex(ti_timeseries['dti']).to_julian_date()
        ti_timeseries['site'] = merged_sites_fww_ti.loc[idx,'properties.id_right']
        ti_timeseries['distance'] = merged_sites_fww_ti.loc[idx,'distance']

        # download the FreshWater Watch (fww) timeseries data
        response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_sites_fww_ti.loc[idx,'properties.id_left']}&datasetID=freshwater_watch&determinand=chemical_phosphate")
        fww_timeseries=pd.json_normalize(response.json(), record_path = 'data')

        # format the fww data as a pandas dataframe
        if len(fww_timeseries) > 0:
            fww_timeseries['dti'] = pd.DatetimeIndex(fww_timeseries['datetime']).tz_convert(None)
            fww_timeseries['date'] = fww_timeseries['dti'].dt.date
            fww_timeseries['jul2'] = pd.DatetimeIndex(fww_timeseries['dti']).to_julian_date()
            fww_timeseries['site'] = merged_sites_fww_ti.loc[idx,'properties.id_right']
            # merge the ti and fww data based on sampling time
            # tolerance = 1 = 24 hours
            merged_timeseries_fww_ti = pd.concat([merged_timeseries_fww_ti,pd.merge_asof(ti_timeseries.sort_values('jul2'),fww_timeseries.sort_values('jul2'), left_on='jul2', right_on='jul2', direction='nearest', tolerance=1)])

# calculate how far apart the samples were in time
merged_timeseries_fww_ti['timedelta'] = merged_timeseries_fww_ti['dti_x'] - merged_timeseries_fww_ti['dti_y']
merged_timeseries_fww_ti=merged_timeseries_fww_ti.filter(['value_x','value_y','dti_x','dti_y','site_x','site_y', 'distance', 'timedelta'])
merged_timeseries_fww_ti=merged_timeseries_fww_ti.rename(columns={'value_x':'ti_value', 'value_y':'fww_value','dti_x':'ti_dti','dti_y':'fww_dti','site_x':'ti_site','site_y':'fww_site'})
merged_timeseries_fww_ti = merged_timeseries_fww_ti[merged_timeseries_fww_ti['fww_value'].notnull()]
merged_timeseries_fww_ti_dd = merged_timeseries_fww_ti.drop_duplicates()
merged_timeseries_fww_ti_dd = merged_timeseries_fww_ti_dd.reset_index()
merged_timeseries_fww_ti_dd['hours_diff'] = abs(merged_timeseries_fww_ti_dd['timedelta'].dt.total_seconds()/3600)

# assess whether the results agree or not
fwwDictMin = {"<0.02":0,"0.02-0.05": 0.02,"0.05-0.1":0.05,"0.1-0.2":0.1,"0.2-0.5":0.2,"0.5-1":0.5}
fwwDictMax = {"<0.02":.02,"0.02-0.05": 0.05,"0.05-0.1":0.1,"0.1-0.2":0.2,"0.2-0.5":0.5,"0.5-1":1}
merged_timeseries_fww_ti_dd['agree'] = False
for idx in merged_timeseries_fww_ti_dd.index:
    fwwcat = merged_timeseries_fww_ti_dd.loc[idx,'fww_value']
    if merged_timeseries_fww_ti_dd.loc[idx,'ti_value'] >= fwwDictMin[fwwcat] and merged_timeseries_fww_ti_dd.loc[idx,'ti_value'] <= fwwDictMax[fwwcat]:
        merged_timeseries_fww_ti_dd.loc[idx,'agree'] = True

# plot the results
fig, axes = plt.subplots(nrows=2)
# strip plot
sns.stripplot(x='fww_value',y='ti_value',data=merged_timeseries_fww_ti_dd, order=["<0.02","0.02-0.05","0.05-0.1","0.1-0.2","0.2-0.5","0.5-1"], hue='agree', ax=axes[0])
axes[0].set_xlabel('FWW phosphate concentration (mg/L)')
axes[0].set_ylabel('TI soluble reactive \n phosphorus concentration (mg/L)')
axes[0].set_yticks([0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.6])
axes[0].axhline(y=0.02,ls='--')
axes[0].axhline(y=0.05,ls='--')
axes[0].axhline(y=0.1,ls='--')
axes[0].axhline(y=0.2,ls='--')
axes[0].axhline(y=0.5,ls='--')
axes[0].text(x=-.4,y=.85,s="A")
# timeline
# Sort by ti_dti
df = merged_timeseries_fww_ti_dd.sort_values(by='ti_dti')
# stagger the stems 
levels = [1, 2, 3, -1, -2, -3] * len(df)
levels = levels[:len(df)]
# Draw horizontal line
axes[1].axhline(0, color='black', linewidth=1)
# Draw stems and points
sns.scatterplot(ax=axes[1],x='ti_dti', y=[0]*len(df), hue='agree', s=50, zorder=3,data=df, legend=False, edgecolor='k')
ci = 0 
for i in df.index:
    if pd.notna(df.loc[i,'ti_dti']):
        axes[1].vlines(df.loc[i,'ti_dti'], 0, levels[ci], linewidth=1)
    ci += 1 
# Ensure annotations aren't too long
df['annot'] = df["ti_site"].str.slice(0,10)
# Annotate each point
for d, l, sname, di, hd in zip(df['ti_dti'], levels, df['annot'],df['distance'],df['hours_diff']):
        di = np.round(di,2)
        hd = np.round(hd,2)
        axes[1].text(d, l, f"{sname}", ha='center', va='bottom' if l > 0 else 'top', fontsize=8, rotation=30)
# Format x-axis
axes[1].xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
# Remove y-axis 
axes[1].get_yaxis().set_visible(False)
#remove some spines
axes[1].spines[["left", "top", "right"]].set_visible(False)
axes[1].set_ylim([-7,5])
axes[1].set_xticks(ticks=axes[1].get_xticks(),rotation=45,labels=axes[1].get_xticklabels())
axes[1].set_xlabel("")
axes[1].text(x=axes[1].get_xticks()[0],y=-5,s=f"μ distance (m): { np.round(df['distance'].mean(),2)} \nμ time difference (hrs): { np.round(df['hours_diff'].mean(),2)}")
axes[1].text(x=axes[1].get_xticks()[0],y=5,s=f"    B")

fig.tight_layout()
fig.set_size_inches(10,10)
fig.text(.5, -.05, "Figure 2. A: A scatter plot comparison of FWW and TI phosphate measurements. \n Scatter points marked orange agree between the two datasets. \n B: A timeline showing when and where the samples were matched, labelled using the site name and date corresponding to the TI data in each pair. \n Timepoints coloured orange showed agreement between the two datasets. \n The average distance and time between the samples is given in text.", ha='center')

# make a table to show accuracy of different FWW classifications
fig, ax = plt.subplots()
# table
summary = (
    merged_timeseries_fww_ti_dd
    .groupby('fww_value')['agree']
    .agg(
        n='count',
        proportion_agree='mean'
    )
)
order = ["<0.02","0.02-0.05","0.05-0.1","0.1-0.2","0.2-0.5","0.5-1"]
summary = summary.reindex(order)
table_df = summary.copy()
table_df['proportion_agree'] = table_df['proportion_agree'].map('{:.2f}'.format)
table_df['n'] = table_df['n'].astype(int)
ax.axis('off')
ax.set_title("FWW Kit Accuracy Table")
tbl = ax.table(
    cellText=table_df.values,
    colLabels=['Count (n)', 'Proportion correct'],
    rowLabels=table_df.index,
    loc='center',
    cellLoc='center'
)
# Formatting
tbl.auto_set_font_size(False)
tbl.set_fontsize(10)
tbl.scale(1, 1.4)

# Header styling
for (row, col), cell in tbl.get_celld().items():
    if row == 0:
        cell.set_text_props(weight='bold')
        cell.set_facecolor('#E6E6E6')
    if col == -1:  # row labels
        cell.set_text_props(weight='bold')
        cell.set_facecolor('#F5F5F5')

fig.text(.5, .7, "Table 1. A: A comparison of FWW and TI phosphate measurements.\nThe FWW citizen return a result as a range of possible phosphate values.\nThis table presents the count of each type of result after matching to the TI data (Count (n) column)\nand what proportion of these results were correct based on the correpsonding TI result (Proportion correct column).", ha='center')
   

merged_timeseries = pd.DataFrame()
for idx in merged_site_wq_fww.index:
    
    #download the fww timeseries data
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_site_wq_fww.loc[idx,'properties.id_right']}&datasetID=freshwater_watch&determinand=chemical_phosphate")
    fww_timeseries=pd.json_normalize(response.json(), record_path = 'data')

    # format into a dataframe
    if len(fww_timeseries) > 0:
        fww_timeseries['dti'] = pd.DatetimeIndex(fww_timeseries['datetime']).tz_convert(None)
        fww_timeseries['date'] = fww_timeseries['dti'].dt.date
        fww_timeseries['jul2'] = pd.DatetimeIndex(fww_timeseries['dti']).to_julian_date()
        fww_timeseries['site'] = merged_site_wq_fww.loc[idx,'properties.id_right']
        fww_timeseries['distance'] = merged_site_wq_fww.loc[idx,'distance']

        # download the ea timeseries data
        response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={merged_site_wq_fww.loc[idx,'properties.id_left']}&datasetID=ea_water_quality&determinand=0180")
        eawq_timeseries=pd.json_normalize(response.json(), record_path = 'data')

        # format into a dataframe
        if len(eawq_timeseries)>0:
            eawq_timeseries['dti'] = pd.DatetimeIndex(eawq_timeseries['datetime'])
            eawq_timeseries['date'] = eawq_timeseries['dti'].dt.date
            eawq_timeseries['jul2'] = pd.DatetimeIndex(eawq_timeseries['dti']).to_julian_date()
            eawq_timeseries['site'] = merged_site_wq_fww.loc[idx,'properties.id_left']
            #the line below merges the fww and ti data based on time, and ensures dates are no more than one day apart. 
            merged_timeseries = pd.concat([merged_timeseries,pd.merge_asof(fww_timeseries.sort_values('jul2'),eawq_timeseries.sort_values('jul2'), left_on='jul2', right_on='jul2', direction='nearest', tolerance=1)])
            
# calculate the time difference between merged samples 
merged_timeseries['timedelta'] = merged_timeseries['dti_x'] - merged_timeseries['dti_y']
merged_timeseries=merged_timeseries.filter(['value_x','value_y','dti_x','dti_y','site_x','site_y', 'distance', 'timedelta'])
merged_timeseries=merged_timeseries.rename(columns={'value_x':'fww_value', 'value_y':'eawq_value','dti_x':'fww_dti','dti_y':'eawq_dti','site_x':'fww_site','site_y':'eawq_site'})
merged_timeseries = merged_timeseries[merged_timeseries['eawq_value'].notnull()]
merged_timeseries['hours_diff'] = abs(merged_timeseries['timedelta'].dt.total_seconds()/3600)
merged_timeseries.reset_index(inplace=True)

fwwDictMin = {"<0.02":0,"0.02-0.05": 0.02,"0.05-0.1":0.05,"0.1-0.2":0.1,"0.2-0.5":0.2,"0.5-1":0.5}
fwwDictMax = {"<0.02":.02,"0.02-0.05": 0.05,"0.05-0.1":0.1,"0.1-0.2":0.2,"0.2-0.5":0.5,"0.5-1":1}
merged_timeseries['agree'] = False
for idx in merged_timeseries.index:
    fwwcat = merged_timeseries.loc[idx,'fww_value']
    # deal with string ea values
    if merged_timeseries.loc[idx,'eawq_value'] == '<0.01':
        merged_timeseries.loc[idx,'eawq_value'] = 0.01
    if merged_timeseries.loc[idx,'eawq_value'] >= fwwDictMin[fwwcat] and merged_timeseries.loc[idx,'eawq_value'] <= fwwDictMax[fwwcat]:
        merged_timeseries.loc[idx,'agree'] = True

# plot the results
fig, axes = plt.subplots(nrows=2)
# strip plot
sns.stripplot(x='fww_value',y='eawq_value',data=merged_timeseries, order=["<0.02","0.02-0.05","0.05-0.1","0.1-0.2","0.2-0.5","0.5-1"], hue='agree', ax=axes[0])
axes[0].set_xlabel('FWW phosphate concentration (mg/L)')
axes[0].set_ylabel('EA phosphate concentration (mg/L)')
axes[0].set_yticks([0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.6])
axes[0].axhline(y=0.02,ls='--')
axes[0].axhline(y=0.05,ls='--')
axes[0].axhline(y=0.1,ls='--')
axes[0].axhline(y=0.2,ls='--')
axes[0].axhline(y=0.5,ls='--')
axes[0].text(x=-.4,y=.55,s="A")

# timeline
# Sort by ti_dti
df = merged_timeseries.sort_values(by='eawq_dti')
# stagger the stems 
levels = [1, 2, 3, -1, -2, -3] * len(df)
levels = levels[:len(df)]
# Draw horizontal line
axes[1].axhline(0, color='black', linewidth=1)
# Draw stems and points
sns.scatterplot(ax=axes[1],x='eawq_dti', y=[0]*len(df), hue='agree', s=50, zorder=3,data=df, legend=False, edgecolor='k')
ci = 0 
for i in df.index:
    if pd.notna(df.loc[i,'eawq_dti']):
        axes[1].vlines(df.loc[i,'eawq_dti'], 0, levels[ci], linewidth=1)
    ci += 1 
# Ensure annotations aren't too long
df['annot'] = df["fww_site"].str.slice(0,10)
# Annotate each point
for d, l, sname, di, hd in zip(df['eawq_dti'], levels, df['annot'],df['distance'],df['hours_diff']):
        di = np.round(di,2)
        hd = np.round(hd,2)
        axes[1].text(d, l, f"{sname}", ha='center', va='bottom' if l > 0 else 'top', fontsize=8, rotation=30)
# Format x-axis
axes[1].xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
# Remove y-axis 
axes[1].get_yaxis().set_visible(False)
#remove some spines
axes[1].spines[["left", "top", "right"]].set_visible(False)
axes[1].set_ylim([-7,5])
axes[1].set_xticks(ticks=axes[1].get_xticks(),rotation=45,labels=axes[1].get_xticklabels())
axes[1].set_xlabel("")
axes[1].text(x=axes[1].get_xticks()[0],y=-5,s=f"μ distance (m): { np.round(df['distance'].mean(),2)} \nμ time difference (hrs): { np.round(df['hours_diff'].mean(),2)}")
fig.text(x=0.17,y=.55,s="B")

fig.tight_layout()
fig.set_size_inches(10,10)
fig.text(.5, -.05, "Figure 2. A: A scatter plot comparison of FWW and EA phosphate measurements. \n Scatter points marked orange agree between the two datasets. \n B: A timeline showing when and where the samples were matched, labelled using the site name and date corresponding to the EA data in each pair. \n Timepoints coloured orange showed agreement between the two datasets. \n The average distance and time between the samples is given in text.", ha='center')

# make a table to show accuracy of different FWW classifications
fig, ax = plt.subplots()
# table
summary = (
    merged_timeseries
    .groupby('fww_value')['agree']
    .agg(
        n='count',
        proportion_agree='mean'
    )
)
order = ["<0.02","0.02-0.05","0.05-0.1","0.1-0.2","0.2-0.5","0.5-1"]
summary = summary.reindex(order)
table_df = summary.copy()
table_df['proportion_agree'] = table_df['proportion_agree'].map('{:.2f}'.format)
table_df['n'] = table_df['n'].astype(float)
ax.axis('off')
ax.set_title("FWW Kit Accuracy Table")
tbl = ax.table(
    cellText=table_df.values,
    colLabels=['Count (n)', 'Proportion correct'],
    rowLabels=table_df.index,
    loc='center',
    cellLoc='center'
)
# Formatting
tbl.auto_set_font_size(False)
tbl.set_fontsize(10)
tbl.scale(1, 1.4)

# Header styling
for (row, col), cell in tbl.get_celld().items():
    if row == 0:
        cell.set_text_props(weight='bold')
        cell.set_facecolor('#E6E6E6')
    if col == -1:  # row labels
        cell.set_text_props(weight='bold')
        cell.set_facecolor('#F5F5F5')
fig.text(.5, .7, "Table 2. A: A comparison of FWW and EA phosphate measurements.\nThe FWW citizen return a result as a range of possible phosphate values.\nThis table presents the count of each type of result after matching to the EA data (Count (n) column)\nand what proportion of these results were correct based on the correpsonding EA result (Proportion correct column).", ha='center')

# Download all TI measurements
ti_all = []
for idx in thamesinit_sites.index:
    # download the thames initiative time series data for each site
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={thamesinit_sites.loc[idx,'properties.id']}&datasetID=thames_initiative&determinand=Soluble%20reactive%20phosphorus%20(%CE%BCg/L)")
    ti_all.append(pd.json_normalize(response.json(), record_path = 'data'))
# concat it all together
ti_df = pd.concat(ti_all).reset_index()
# convert non numerical data to nan
ti_df['value'] = pd.to_numeric(ti_df['value'], errors='coerce')
ti_df['value_mg'] = ti_df['value']/1000

# Download all EA measurements
ea_all = []
for idx in eawq_sites.index:
    # download the thames initiative time series data for each site
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={eawq_sites.loc[idx,'properties.id']}&datasetID=ea_water_quality&determinand=0180")
    ea_all.append(pd.json_normalize(response.json(), record_path = 'data'))
# concat it all together
ea_df = pd.concat(ea_all).reset_index()
# convert non numerical data to nan
ea_df['value'] = pd.to_numeric(ea_df['value'], errors='coerce')

# Download all FWW measurements.
fww_all = []
for idx in fww_sites.index:
    # download the thames initiative time series data for each site
    response = requests.get(f"https://oxfordrivers.ceh.ac.uk/getTimeseries?siteID={fww_sites.loc[idx,'properties.id']}&datasetID=freshwater_watch&determinand=chemical_phosphate")
    fww_all.append(pd.json_normalize(response.json(), record_path = 'data'))
# concat it all together
fww_df = pd.concat(fww_all).reset_index()
fww_df['value']=pd.Categorical(fww_df['value'], categories=['<0.02','0.02-0.05','0.05-0.1','0.1-0.2','0.2-0.5','0.5-1','>1','No_water_available',''])
fww_df.drop(fww_df[fww_df['value'] == 'No_water_available'].index, inplace=True)
fww_df.drop(fww_df[fww_df['value'] == ''].index, inplace=True)
fww_df['value']=pd.Categorical(fww_df['value'], categories=['<0.02','0.02-0.05','0.05-0.1','0.1-0.2','0.2-0.5','0.5-1','>1'])

# Bin EA and TI data so it matches FWW
bins = [0, 0.02, 0.05, 0.1, 0.2, 0.5, 1, np.inf]
labels = [
    "0–0.02", "0.02–0.05", "0.05–0.1",
    "0.1–0.2", "0.2–0.5", "0.5–1", ">1"
]
ea_df["value_bin"] = pd.cut(ea_df["value"], bins=bins, labels=labels)
ti_df["value_bin"] = pd.cut(ti_df["value_mg"], bins=bins, labels=labels)

# Now make comparing the counts of the relative bins 
fig, axes = plt.subplots(nrows=3, figsize=(8,16))

sns.countplot(data=fww_df,x='value',ax=axes[0])
ax02 = axes[0].twinx()
sns.countplot(data=fww_df,x='value',ax=ax02, stat='percent')
axes[0].set_xlabel("")
axes[0].text(s="FWW",x=.75,y=.8,transform=axes[0].transAxes,color='tab:blue', size=18)

sns.countplot(data=ti_df,x='value_bin',ax=axes[1],color='tab:orange')
ax12 = axes[1].twinx()
sns.countplot(data=ti_df,x='value_bin',ax=ax12, stat='percent',color='tab:orange')
axes[1].set_xlabel("")
axes[1].text(s="TI",x=.75,y=.8,transform=axes[1].transAxes,color='tab:orange', size=18)

sns.countplot(data=ea_df,x='value_bin',ax=axes[2],color='tab:green')
ax22 = axes[2].twinx()
sns.countplot(data=ea_df,x='value_bin',ax=ax22, stat='percent',color='tab:green')
axes[2].text(s="EA",x=.75,y=.8,transform=axes[2].transAxes,color='tab:green', size=18)

axes[2].set_xlabel("Phosphate concentration (mg/L)")
fig.set_tight_layout(tight=True)
fig.text(.5, -.05, "Figure 3. Distributions of phosphate concentrations in the entire FWW (top), TI (middle) and EA (bottom) datasets on the Oxford Rivers Portal. \n The continuous TI and EA datasets have been discretised to match the categories on the FWW test strips.", ha='center', size =14)