## Process EWR Hg Data into a surface # AJMR: September 2022 #%% Setup import pandas as pd import geopandas as gp import gpxpy import gpxpy.gpx import numpy as np import matplotlib.pyplot as plt import contextily as ctx import datetime import scipy as sp from scipy.interpolate import griddata from shapely import geometry, ops #%% Read in centerline shapefile river_centerline = gp.read_file('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/16_Waterline/Centerline_for_Modelling_UTMZ15.shp') river_centerlineExploded = river_centerline.explode(ignore_index=True) river_centerlineExploded.reset_index(inplace=True) tempMulti = river_centerlineExploded.iloc[[5,0,1,2,3,4,6,7,9], 4] # Put the sub-line coordinates into a list of sublists outcoords = [list(i.coords) for i in tempMulti] # Flatten the list of sublists and use it to make a new line river_centerline_merge = geometry.LineString([i for sublist in outcoords for i in sublist]) river_centerline_merge_gpd = gp.GeoSeries(river_centerline_merge) river_centerline_merge_gpd2 =\ gp.GeoDataFrame(geometry=gp.points_from_xy( river_centerline_merge.xy[0], river_centerline_merge.xy[1], crs="EPSG:32615")) # Add distance along centerline river_centerline_merge_gpd2['DistanceFromPrevious'] = river_centerline_merge_gpd2.distance(river_centerline_merge_gpd2.shift(1)) river_centerline_merge_gpd2['RiverKM'] = river_centerline_merge_gpd2['DistanceFromPrevious'].cumsum() river_centerline_merge_gpd2.iloc[0, 1] = 0 river_centerline_merge_gpd2.iloc[0, 2] = 0 #%% Read in combined dataset # obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/output/combined dataset-SGJ.csv") # obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/combined dataset.csv") # obs_IN = pd.read_csv("C:/Users/arey/Downloads/ExportDownloads_EwrrpDataExport_rharris_2023-05-12_104911/EwrrpDataExport_rharris_2023-05-12_104911.txt", sep='\t', encoding='UTF-16 LE') df = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Results.txt", sep = '\t', encoding='UTF-16 LE') locations = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Sites.txt", sep = '\t', encoding='UTF-16 LE') samples = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Samples.txt", sep = '\t', encoding='UTF-16 LE') #%% Merge locations and results obs_IN = pd.merge(left=df, right=locations, left_on= 'Samplesiteid', right_on= 'Id', how = 'left') # 1096 locations w/Hg obs_IN = pd.merge(left=obs_IN, right=samples, left_on= 'Sampleid', right_on= 'Id', how = 'left') # 1096 locations w/Hg # Add Datetime obs_IN['Sampledate'] = pd.to_datetime(obs_IN.Sampledate) #%% Process combined datsaset # Convert to geodataframe obs_gdf = gp.GeoDataFrame(obs_IN, geometry=gp.points_from_xy( obs_IN.loc[:, 'Longitude'], obs_IN.loc[:, 'Latitude'], crs="EPSG:4326")).to_crs(crs="EPSG:32615") # Add centerline and reset index obs_gdf = gp.sjoin_nearest(river_centerline_merge_gpd2, obs_gdf, how='right', max_distance=250).reset_index() # Sediment Hg GeoDataFrame where media is Sediment obsMask = (((obs_gdf.Media == 'SED') | (obs_gdf.Media == 'SOIL')) & ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna() & ((obs_gdf.TopDepth == 0) | (obs_gdf.TopDepth.isna()))) # obsMask = (((obs_gdf.Media_x == 'SED') | (obs_gdf.Media_x == 'SOIL')) & # ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') # | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & # (obs_gdf.DateTime > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna()) obs_HgSed = obs_gdf.loc[obsMask, :].copy() # Adjust Units obs_HgSed.loc[obs_HgSed.Unit == 'NG/G', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'NG/G', 'Samplevalue'] obs_HgSed.loc[obs_HgSed.Unit == 'ng/g', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'ng/g', 'Samplevalue'] obs_HgSed.loc[obs_HgSed.Unit == 'UG/G', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'UG/G', 'Samplevalue'] * 1000 obs_HgSed.loc[obs_HgSed.Unit == 'ug/g', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'ug/g', 'Samplevalue'] * 1000 obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Samplevalue'] * 0.001 obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Samplevalue'] * 0.001 # Save subset of columns as a shapefile obs_HgSed.loc[:, ['Sample_NG/G', 'Sampledate_x', 'RiverKM', 'Latitude', 'Longitude', 'geometry', 'Surveyname','Sitecode','Parameter']].to_file( 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgSed.shp') #%% Extract Hg data from surface water samples # Water Hg GeoDataFrame where media is SurfaceWater (SW) obsMask = (((obs_gdf.Media == 'SW')) & ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)-Total') | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna()) obs_HgWater = obs_gdf.loc[obsMask, :].copy() # Adjust Units obs_HgWater.loc[obs_HgWater.Unit == 'NG/L', 'Samplevalue'] = obs_HgWater.Samplevalue obs_HgWater.loc[obs_HgWater.Unit == 'ng/L', 'Samplevalue'] = obs_HgWater.Samplevalue obs_HgWater.loc[obs_HgWater.Unit == 'UG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 obs_HgWater.loc[obs_HgWater.Unit == 'ug/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 obs_HgWater.loc[obs_HgWater.Unit == 'MG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 obs_HgWater.loc[obs_HgWater.Unit == 'mg/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 # Drop where no result obs_HgWater = obs_HgWater.dropna(subset=['Samplevalue']) # Save subset of columns as a shapefile obs_HgWater.loc[:, ['Samplevalue', 'Sampledate', 'RiverKM', 'geometry', 'Contributor', 'Sitecode', 'Parameter']].to_file( 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgWater_RevB.shp') #%% River Profile Plotting fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 6)) axes.scatter(obs_HgSed.RiverKM, obs_HgSed.Samplevalue) axes.scatter(obs_HgSed.RiverKM, obs_HgSed.Samplevalue) axes.set_ylim(0, 20000) axes.set_xlim(0, 100000) axes.set_xlabel('Distance Along River [m]') axes.set_ylabel('Surface Sediment Hg [ng/g]') plt.show() fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Hg_SurfaceProfile.png', bbox_inches='tight', dpi=300) #%% Plot Vertical Profiles fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(7, 7)) obs_HgSed.plot.scatter('RiverKM', 'TopDepth_x', 12, 'Sample_NG/G', ax=axes, vmax=20000, cmap='viridis') axes.invert_yaxis() axes.set_xlabel('Distance Along River [m]') axes.set_ylabel('Sample Depth [cm]') axes.set_xlim([0, 100000]) # Get colorbar axis to set a legend cax = fig.get_axes()[1] #and we can modify it, i.e.: cax.set_ylabel('Sediment Hg (ng/g)') plt.show() fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Hg_Vertical.png', bbox_inches='tight', dpi=300) #%% Hg Map Plotting mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 7)) fig.patch.set_facecolor('white') # Full system limits # axes.set_xlim(350000, 520000) # axes.set_ylim(5510000, 5580000) # # Clay Lake Limits # axes.set_xlim(466418, 515846) # axes.set_ylim(5513437, 5545362) # Dryden Limits axes.set_xlim(497697, 515846) axes.set_ylim(5513437, 5525166) # Add Basemap ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615') # Plot Hg At surface obs_HgSed.plot(ax=axes, column='Samplevalue', legend=True, vmin=0, vmax=5000, markersize=2, cmap='viridis', legend_kwds={'label': 'Surface Sediment Hg (ng/g)'}) plt.show() fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceHgMap.png', bbox_inches='tight', dpi=300) #%% Interpolate Hg bathy_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Bed_Level.xyz', names=['x', 'y', 'z'], header=0, delim_whitespace=True) rough_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Roughness.xyz', names=['x', 'y', 'z'], header=0, delim_whitespace=True) HgInterp = griddata(np.array([dataOUTmax.geometry.x, dataOUTmax.geometry.y]).T, np.array(dataOUTmax.Samplevalue).T, np.array([bathy_xyz.x, bathy_xyz.y]).T, method='nearest') HgInterpOut = np.vstack((bathy_xyz.x, bathy_xyz.y, HgInterp)) HgInterpOut = np.transpose(HgInterpOut) # Set Wetlands to zero # Interpolate Roughness RoughInterp = sp.interpolate.griddata(np.transpose(np.array([rough_xyz.x, rough_xyz.y])), rough_xyz.z, np.transpose(np.array([bathy_xyz.x, bathy_xyz.y])), method='nearest') HgInterpOut[RoughInterp==0.05, 2] = 0 np.savetxt('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceInterFiltDB.xyz', HgInterpOut, delimiter=" ") #%% Hg Interp Plotting mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 7)) fig.patch.set_facecolor('white') # Clay Lake Limits axes.set_xlim(466418, 515846) axes.set_ylim(5513437, 5545362) # Add Basemap ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615') # Plot Hg axes.scatter(HgInterpOut[:, 0], HgInterpOut[:, 1], 5, HgInterpOut[:, 2], vmin=0, vmax=5000) plt.show()