From 8a78279eaab79b01d71bf849704bd20632890391 Mon Sep 17 00:00:00 2001 From: Alexander Rey Date: Thu, 29 Jun 2023 16:55:43 -0400 Subject: [PATCH] Final Baird Commit --- Azure/secret.py | 0 EWR_D3D/WSP_DB.py | 0 EWR_Data/EWR_Hg_Syn.py | 238 ++++++++++++++++ EWR_Data/createSeg.py | 0 NTC_DFM/Delft3DBoundsCompare.py | 0 NTC_DFM/NTC_VelAnimation.py | 464 ++++++++++++++++++++++++++++++++ SatBathy/.idea/.gitignore | 0 SatBathy/.idea/SatBathy.iml | 8 + 8 files changed, 710 insertions(+) create mode 100644 Azure/secret.py create mode 100644 EWR_D3D/WSP_DB.py create mode 100644 EWR_Data/EWR_Hg_Syn.py create mode 100644 EWR_Data/createSeg.py create mode 100644 NTC_DFM/Delft3DBoundsCompare.py create mode 100644 NTC_DFM/NTC_VelAnimation.py create mode 100644 SatBathy/.idea/.gitignore create mode 100644 SatBathy/.idea/SatBathy.iml diff --git a/Azure/secret.py b/Azure/secret.py new file mode 100644 index 0000000..e69de29 diff --git a/EWR_D3D/WSP_DB.py b/EWR_D3D/WSP_DB.py new file mode 100644 index 0000000..e69de29 diff --git a/EWR_Data/EWR_Hg_Syn.py b/EWR_Data/EWR_Hg_Syn.py new file mode 100644 index 0000000..11faf1d --- /dev/null +++ b/EWR_Data/EWR_Hg_Syn.py @@ -0,0 +1,238 @@ +## 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() + + + + diff --git a/EWR_Data/createSeg.py b/EWR_Data/createSeg.py new file mode 100644 index 0000000..e69de29 diff --git a/NTC_DFM/Delft3DBoundsCompare.py b/NTC_DFM/Delft3DBoundsCompare.py new file mode 100644 index 0000000..e69de29 diff --git a/NTC_DFM/NTC_VelAnimation.py b/NTC_DFM/NTC_VelAnimation.py new file mode 100644 index 0000000..4efd0b6 --- /dev/null +++ b/NTC_DFM/NTC_VelAnimation.py @@ -0,0 +1,464 @@ +""" +@author: AJMR + +Plotting and animation script for NTC sediment tests + +April 17, 2023 +""" + +import os +import pandas as pd +import geopandas as gpd +import netCDF4 as nc +import math +import numpy as np + +import matplotlib as mpl +import matplotlib.pyplot as plt +from scipy.interpolate import LinearNDInterpolator, interp1d +import matplotlib.animation as animation + +import dfm_tools as dfmt + +# Geotiff +import cartopy.crs as ccrs +import contextily as ctx +import pathlib as pl + +from shapely.geometry import Point, MultiPoint +import alphashape + +import rioxarray +import xarray as xr + + +import cartopy.crs as ccrs +import contextily as ctx +from dfm_tools.regulargrid import scatter_to_regulargrid + +from pathlib import Path + +# Colormaps +import cm_xml_to_matplotlib as cm + +#%% Read Model Log + +runLog = pd.read_excel('C:/Users/arey/files/Projects/Newtown/Model Log NTC.xlsx', 'ModelLog') + +dataPath = "FlowFM/dflowfm/output/FlowFM_map.nc" +histPath = "FlowFM/dflowfm/output/FlowFM_his.nc" + +#%% Import Model results for static images + +# Define which models to import +modelPlot = [144] + +d3dfm_DataArray = [None] * (max(modelPlot)+1) +sedMass = [None] * (max(modelPlot)+1) # Available Mass of Sediment +erodeSed = [None] * (max(modelPlot)+1) # Erosion and deposition +initialBed = [None] * (max(modelPlot)+1) # Erosion and deposition +modelCRS = [None] * (max(modelPlot)+1) # Model CRS +sedConc = [None] * (max(modelPlot)+1) # Model CRS + +for i in modelPlot: + file_nc_map = Path(runLog['Run Location'][i]) / dataPath + + # Open Map file as xarrray Dataset + d3dfm_DataArray[i] = dfmt.open_partitioned_dataset(file_nc_map.as_posix()) + + # Get Var info + vars_pd = dfmt.get_ncvarproperties(d3dfm_DataArray[i]) + + # Get last timestep + tStep = d3dfm_DataArray[i].time[-1].values + + # Import sand sediment class + sedMass[i] = d3dfm_DataArray[i].mesh2d_bodsed[-1, :, 2].compute() + initialBed[i]= d3dfm_DataArray[i].mesh2d_s1[0, :] - d3dfm_DataArray[i].mesh2d_waterdepth[0, :] + finalBed = d3dfm_DataArray[i].mesh2d_s1[-1, :] - d3dfm_DataArray[i].mesh2d_waterdepth[-1, :] + + erodeSed[i] = finalBed - initialBed[i] + erodeSed[i] = erodeSed[i].compute() + initialBed[i] = initialBed[i].compute() + + # Import concentration of sediment 2 at timestep 95 and layer 10 + sedConc[i] = d3dfm_DataArray[i].mesh2d_sedfrac_concentration[94, 2, :, 9].compute() + + modelCRS[i] = vars_pd['EPSG_code']['projected_coordinate_system'] + +#%% Define axis limits +axLim_Xmin = [] +axLim_Xmax = [] + +axLim_Ymin = [] +axLim_Ymax = [] + +# Set axis limits +# Full domain +axLim_Xmin.append(298500) +axLim_Xmax.append(306800) + +axLim_Ymin.append(58500) +axLim_Ymax.append(68000) + +# Zoom North River +axLim_Xmin.append(303500) +axLim_Xmax.append(305500) + +axLim_Ymin.append(66000) +axLim_Ymax.append(68000) + +# Zoom Tribs +axLim_Xmin.append(305500) +axLim_Xmax.append(306800) + +axLim_Ymin.append(60000) +axLim_Ymax.append(62500) + +# Zoom far trib +axLim_Xmin.append(305800) +axLim_Xmax.append(305900) + +axLim_Ymin.append(60150) +axLim_Ymax.append(60400) + + +#%% Plot using DFM functions + +for i in modelPlot: + # Create figure + fig, axesFig = plt.subplots(1, 4, figsize=(8, 2)) + + axes = axesFig.flatten() + + # Loop through axis and plot with different limits + for subIDX in range(0, 4): + + # # Plot Sediment Mass + # smp = sedMass[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='turbo', + # vmin=0, vmax=2000, add_colorbar=False) + + # # Plot Erosion + # smp = erodeSed[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='coolwarm', + # vmin=-1, vmax=1, add_colorbar=False) + + # Plot Bed + smp = initialBed[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='nipy_spectral', + vmin=-25, vmax=0, add_colorbar=False) + + + # Set axis labels + axes[subIDX].set_xlabel('') + axes[subIDX].set_ylabel('') + axes[subIDX].set_title('') + axes[subIDX].set_xticks([]) + axes[subIDX].set_yticks([]) + + axes[subIDX].set_xlim(axLim_Xmin[subIDX], axLim_Xmax[subIDX]) + axes[subIDX].set_ylim(axLim_Ymin[subIDX], axLim_Ymax[subIDX]) + + # Add basemap + # ctx.add_basemap(ax=axes, source=ctx.providers.Esri.WorldImagery, crs=modelCRS[i], attribution=False) + mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' + ctx.add_basemap(axes[subIDX], source=mapbox, crs=modelCRS[i]) + + + # cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Sand Mass [kg/m2]') + # cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Erosion/ deposition [m]') + cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Initial Bed [mNAVD88]') + + plt.show() + + # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/SedMass_' + runLog['Run'][i] + '.png', + # bbox_inches='tight', dpi=300) + + # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/ErodeDep_' + runLog['Run'][i] + '.png', + # bbox_inches='tight', dpi=300) + + # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/InitialBed_' + runLog['Run'][i] + '.png', + # bbox_inches='tight', dpi=300) + +#%% Plot using DFM functions +#Cmap Path +cmap_path = 'C:/Users/arey/Repo/MATLAB_Q/Downloads/KeyColormaps/' + +for i in modelPlot: + fig, axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': ccrs.epsg(32118)}, figsize=(6, 6)) + fig.patch.set_facecolor('white') + + fig.tight_layout(pad=3.0) + + # Load cmaps + cmap = 'AJMR_Sed5_RevE.xml' + plotcmap = cm.make_cmap(cmap_path + cmap) + + pc = plot_netmapdata(ugrid_all[i].verts, values=modelMaxShear[i][0, :], + ax=axes, linewidth=0.5, cmap=plotcmap) + + axes.set_xlim(305900, 306400) + axes.set_ylim(61500, 62200) + + # axes.quiver(modelX[i], modelX[i], U[i], V[i], color='w', scale=1.00) + + mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' + ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118') + + extent = axes.get_extent() + axes.set_xticks(np.linspace(extent[0], extent[1], 4)) + axes.set_yticks(np.linspace(extent[2], extent[3], 4)) + + pc.set_clim([0, 0.145]) + cbarTicks = [0, 0.0378, 0.063, 0.0826, 0.11, 0.145] + cbar = fig.colorbar(pc, ax=axes, shrink=0.95, + ticks=cbarTicks) + cbar.set_label('Total Max Shear Stress [N/m^2]') + + # Add label on colorbar + cbar.ax.text(0.08, (cbarTicks[0]+cbarTicks[1])/2, 'CLAY', ha='center', va='center', + rotation=90, fontweight='bold') + cbar.ax.text(0.08, (cbarTicks[1]+cbarTicks[2])/2, 'FSILT', ha='center', va='center', + rotation=90, fontweight='bold') + cbar.ax.text(0.08, (cbarTicks[2]+cbarTicks[3])/2, 'MSILT', ha='center', va='center', + rotation=90, fontweight='bold') + cbar.ax.text(0.08, (cbarTicks[3]+cbarTicks[4])/2, 'CSILT', ha='center', va='center', + rotation=90, fontweight='bold') + cbar.ax.text(0.08, (cbarTicks[4]+cbarTicks[5])/2, 'FSAND', ha='center', va='center', + rotation=90, fontweight='bold') + + plt.show() + + fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs/Shear_Class_' + runLog['Run'][i] + '.png', + bbox_inches='tight', dpi=300) + +#%% Import Model results for animations + +# modelPlot = [20, 22] +# modelPlot = [18, 20, 21, 22] +# modelPlot = [18, 21] +modelPlot = [25, 26] + +modelSedConc = [None] * (max(modelPlot)+1) +modelMaxShear_animate = [None] * (max(modelPlot)+1) +ugrid_all = [None] * (max(modelPlot)+1) +tSteps = [None] * (max(modelPlot)+1) + +for i in modelPlot: + file_nc_map = Path(runLog['Run Location'][i]) / dataPath + + # Get Var info + vars_pd, dims_pd = get_ncvardimlist(file_nc=file_nc_map.as_posix()) + + # Import + tSteps[i] = get_timesfromnc(file_nc=file_nc_map.as_posix(), varname='time') + + ugrid_all[i] = get_netdata(file_nc=file_nc_map.as_posix()) + modelSedConc[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_sedfrac_concentration', + timestep='all', station=0, layer=0) #timestep='all' OR timestep=range(0, 30) + modelMaxShear_animate[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_tausmax', + timestep='all') + + +# %% Import water levels from hist +modelHistWL = [None] * (max(modelPlot)+1) +modelHistTime = [None] * (max(modelPlot)+1) +# tSteps = [None] * (max(modelPlot)+1) + +for i in modelPlot: + file_nc_hist = Path(runLog['Run Location'][i]) / histPath + + vars_pd, dims_pd = get_ncvardimlist(file_nc=file_nc_hist.as_posix()) + + # Get station names + histStations = get_hisstationlist(file_nc=file_nc_hist.as_posix(), varname='waterlevel') + + # Import + modelHistTime[i] = get_timesfromnc(file_nc=file_nc_hist.as_posix(), varname='time') + modelHistWL[i] = get_ncmodeldata(file_nc=file_nc_hist.as_posix(), varname='waterlevel', + timestep='all', station=1) + + +#%% Sediment animations +plt.rcParams['animation.ffmpeg_path'] = \ + 'C:/Users/arey/Local/ffmpeg-2022-02-14-git-59c647bcf3-full_build/bin/ffmpeg.exe' + +mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckex9vtri0o6619p55sl5qiyv/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' +x, y, arrow_length = 0.93, 0.95, 0.12 +# fontprops = fm.FontProperties(size=12) + +# Set model to plot +# modelPlot = [20, 22] +# modelPlot = [18, 20, 21, 22] +modelPlot = [25, 26] +cmap_path = 'C:/Users/arey/Repo/MATLAB_Q/Downloads/KeyColormaps/' + +for m in modelPlot: + vmax = 0.75 + vmin = 0 + # vmax = 0.145 + # vmin = 0 + + cbarTicks = [0, 0.0378, 0.063, 0.0826, 0.11, 0.145] + + cmapName = 'turbo' + cmap = mpl.cm.turbo + # cmapName = 'AJMR_Sed5_RevE.xml' + # cmap = cm.make_cmap(cmap_path + cmapName) + + + # Zoom limits + # xLimits = [305900, 306500] + # yLimits = [61400, 62200] + + # Wide Limits + xLimits = [302719.4, 306934.2] + yLimits = [60263.7, 64194.8] + + # Setup Video + metadata = dict(title='NTC Sediment Animation', artist='Matplotlib', + comment='AJMR June 28, 2022') + writer = animation.FFMpegWriter(fps=2, metadata=metadata, codec='h264', bitrate=5000) + fig, axes = plt.subplots(figsize=(6, 6)) + + writer.setup(fig, '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + + '06_Models/04_Delft3D/Figures/SedConc/HQ3_Wide_RevC_Long_Sed_Cond_' + + runLog['Run'][m] + '.mp4') + + mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' + + for i in np.arange(1, len(modelSedConc[m])-1, 1): #289 + # Clear axes for video + axes.cla() + + # Create new figure for stills + # fig, axes = plt.subplots(figsize=(6, 6)) + + # Clear inset WL Plot + # if i > 1: + # axin1.cla() + + # Plot sediment + pc = plot_netmapdata(ugrid_all[m].verts, values=modelSedConc[m][i, 0, :, 0], + ax=axes, cmap=cmapName, + antialiaseds=False) + + # Plot shear + # pc = plot_netmapdata(ugrid_all[m].verts, values=modelMaxShear_animate[m][i,:], + # ax=axes, cmap=cmap, + # antialiaseds=False) + + # Set map limits + axes.set_xlim(xLimits[0], xLimits[1]) + axes.set_ylim(yLimits[0], yLimits[1]) + + # Add basemap + ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118') + + # Set map ticks + axes.set_xticks(np.linspace(xLimits[0], xLimits[1], 4)) + axes.set_yticks(np.linspace(yLimits[0], yLimits[1], 4)) + + # Color Limits + pc.set_clim([vmin, vmax]) + + # Add title + axes.title.set_text('Bottom Sediment Concentration at: ' + tSteps[m][i].strftime("%Y/%m/%d %H:%M")) + # axes.title.set_text('Total Max Shear Stress at: ' + str(tSteps[m][i])) + + # Add inset wl plot + axin1 = axes.inset_axes([0.1, 0.1, 0.4, 0.15]) + axin1.plot(modelHistTime[m], modelHistWL[m]) + # Add current point + # Find closest time + abs_deltas_from_target_date = np.absolute(modelHistTime[m] - tSteps[m][i]) + axin1.scatter(modelHistTime[m][np.argmin(abs_deltas_from_target_date)], modelHistWL[m][np.argmin(abs_deltas_from_target_date)], 10, 'r') + + # axin1.set_xticks([modelHistTime[m][0], modelHistTime[m][len(modelHistTime[m])-1]]) + axin1.set_xticks([]) + axin1.set_yticks([]) + + plt.setp(axin1.get_xticklabels(), backgroundcolor="white") + plt.setp(axin1.get_yticklabels(), backgroundcolor="white") + + if i == 1: #i < 1000 + fig.subplots_adjust(right=0.8) + norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) + + # Create colorbar axis + # cax = plt.axes([0.82, 0.25, 0.04, 0.5]) + cax = plt.axes([0.85, 0.15, 0.04, 0.7]) + + # Add colorbar with designated tick marks + # cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap, + # norm=norm, + # orientation='vertical', + # ticks=cbarTicks) + + # Add colorbar without designated tick marks + cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap, + norm=norm, + orientation='vertical') + + # Colorbar lavel + cb1.set_label('Sediment Concentration [$kg/m^3$]') + # cb1.set_label('Total Max Shear Stress [N/m^2]') + + # Add label on colorbar + # cb1.ax.text(0.08, (cbarTicks[0] + cbarTicks[1]) / 2, 'CLAY', ha='center', va='center', + # rotation=90, fontweight='bold') + # cb1.ax.text(0.08, (cbarTicks[1] + cbarTicks[2]) / 2, 'FSILT', ha='center', va='center', + # rotation=90, fontweight='bold') + # cb1.ax.text(0.08, (cbarTicks[2] + cbarTicks[3]) / 2, 'MSILT', ha='center', va='center', + # rotation=90, fontweight='bold') + # cb1.ax.text(0.08, (cbarTicks[3] + cbarTicks[4]) / 2, 'CSILT', ha='center', va='center', + # rotation=90, fontweight='bold') + # cb1.ax.text(0.08, (cbarTicks[4] + cbarTicks[5]) / 2, 'FSAND', ha='center', va='center', + # rotation=90, fontweight='bold') + + + # Change label font size for only the primary axis + for item in ([axes.xaxis.label, axes.yaxis.label] + + axes.get_xticklabels() + axes.get_yticklabels()): + item.set_fontsize(4) + + + # Save video frame + writer.grab_frame() + plt.show() + print(i) + + # Save stills + fig.savefig('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + + '06_Models/04_Delft3D/Figures/SedConc/Wide_RevC_SedConc_' + + runLog['Run'][m] + '_Frame_' + str(i) + '.png', + bbox_inches='tight', dpi=300) + # fig.savefig('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + + # '06_Models/04_Delft3D/Figures/ShearStress/Wide_ShearStress_' + + # runLog['Run'][m] + '_Frame_' + str(i) + '.png', + # bbox_inches='tight', dpi=300) + + # Save video + writer.finish() + + + +#%% Save Sediment Geotiff + +# Create mask within bounds +boundSelect = 2 +croppedGridMask = (d3dfm_DataArray[i].mesh2d_face_x >= axLim_Xmin[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_x <= axLim_Xmax[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_y >= axLim_Ymin[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_y <= axLim_Ymax[boundSelect]) +croppedGridMask = croppedGridMask.compute() + +# Step through time steps +for sedTstep in range(105, 120): + # Rasterize sediment concentration within bounds + sedConcTif = dfmt.rasterize_ugrid(d3dfm_DataArray[i]['mesh2d_sedfrac_concentration'].isel( + time=sedTstep, mesh2d_nLayers=4, nSedSus=0).where(croppedGridMask, drop=True).to_dataset(), + resolution=2) + + # Add CRS + sedConcTif.rio.write_crs("epsg:32118", inplace=True) + sedConcTif.rio.to_raster('C:/Users/arey/files/Projects/Newtown/SedFigs2/GeoTiff_RevB_' + + 'Sed_Conc4_Time_' + str(sedTstep) + '.tif') diff --git a/SatBathy/.idea/.gitignore b/SatBathy/.idea/.gitignore new file mode 100644 index 0000000..e69de29 diff --git a/SatBathy/.idea/SatBathy.iml b/SatBathy/.idea/SatBathy.iml new file mode 100644 index 0000000..d0876a7 --- /dev/null +++ b/SatBathy/.idea/SatBathy.iml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file