335 lines
12 KiB
Python
335 lines
12 KiB
Python
"""
|
||
@author: AJMR
|
||
|
||
Plotting and animation script for NTC
|
||
|
||
June 2, 2022
|
||
"""
|
||
|
||
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
|
||
|
||
from dfm_tools.get_nc import get_netdata, get_ncmodeldata, plot_netmapdata
|
||
from dfm_tools.get_nc_helpers import get_timesfromnc, get_ncvardimlist, get_hisstationlist
|
||
|
||
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', 'Propwash')
|
||
|
||
dataPath = "FlowFM_map.nc"
|
||
histPath = "FlowFM_his.nc"
|
||
|
||
#%% Import Model results for static images
|
||
|
||
# Define which models to import
|
||
#modelPlot = [3, 5, 6, 7, 8]
|
||
#modelPlot = [2, 4, 5, 6, 7, 8]
|
||
# modelPlot = [18, 20, 21, 22]
|
||
modelPlot = [25, 26]
|
||
|
||
ugrid_all = [None] * (max(modelPlot)+1)
|
||
tSteps = [None] * (max(modelPlot)+1)
|
||
modelX = [None] * (max(modelPlot)+1)
|
||
modelY = [None] * (max(modelPlot)+1)
|
||
modelZ = [None] * (max(modelPlot)+1)
|
||
modelMaxShear = [None] * (max(modelPlot)+1)
|
||
|
||
|
||
for i in modelPlot:
|
||
file_nc_map = Path(runLog['Run Location'][i]) / dataPath
|
||
|
||
tSteps[i] = get_timesfromnc(file_nc=file_nc_map.as_posix(), varname='time')
|
||
|
||
# Use final timestep
|
||
tStep = len(tSteps[i])-1
|
||
|
||
# Get Var info
|
||
vars_pd, dims_pd = get_ncvardimlist(file_nc=file_nc_map.as_posix())
|
||
|
||
# Import
|
||
ugrid_all[i] = get_netdata(file_nc=file_nc_map.as_posix())
|
||
modelX[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_face_x')
|
||
modelY[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_face_y')
|
||
modelMaxShear[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_tausmax',
|
||
timestep=tStep)
|
||
|
||
#%% 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()
|
||
|