Gallery 3: Sheet Plots

This page provides examples on:

• 2D and 3D plotting in plan view

• Extracting and plotting maximum statistics

• Plotting the difference between existing and developed conditions

• Exporting of ASCII Arc Grids

• Saving plots to disk

This notebook is used in combination with the TUFLOW FV Python Toolbox (tfv) package.

import numpy as np
import xarray as xr  # We utilise xarray to do all the heavy lifting 
import tfv.xarray
from pathlib import Path # We'll also make use of the `pathlib` module to assist with managing file-paths, although this is entirely optional! 
import matplotlib.pyplot as plt

Open TUFLOW FV Model Result

model_folder = Path('../../data')
model_file = 'HYD_002.nc'

fv = xr.open_dataset(model_folder / model_file, decode_times=False).tfv
# fv  # Uncomment this if you want to review the contents of the Dataframe

Single Sheet Plot

Quick plot - first time step of first variable

fv.plot()

<tfv.visual.SheetPatch at 0x28f1eef63c0>

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

We can also plot any specific time in the model

ts = '2011-05-02 00:00'
fv.plot('H', ts)

<tfv.visual.SheetPatch at 0x28f2092ac10>

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

Sheet Plot With Vectors and Save to Disk

# Plot a variable at a particular time with vectors
ts = '2011-05-02 00:00'

# Add figure and axes
fig, ax = plt.subplots()

# Common settings
ax.set_aspect('equal')

fv.plot('H',ts, ax=ax)
fv.plot_vector('V', ts, ax=ax, color='k')   # We can plot a vector overlay using the `plot_vector` method here. This plot uses depth averaged vectors however you can use the datum and limits variable to change this.

# Save to disk
output_folder = Path('./outputs')
output_folder.mkdir(exist_ok=True) # Make the folder, if it doesn't exist

plt.savefig(output_folder / 'Example_Sheet.jpg', dpi=300)

Sheet Plot at Different Depths

ts = '2011-05-02 00:00'
fig, axes = plt.subplots(ncols=3, figsize=(12, 6), constrained_layout=True)

sal_spec = dict(cmap='ocean', clim=(0, 36), shading='interp', colorbar=False)

ax = axes[0]
p1 = fv.plot('SAL', ts, ax=ax,  datum='depth', limits=(0,1), **sal_spec)  # By setting datum to depth, this is now the "surface 1m")
ax.set_title('Surface 1m')

ax = axes[1]
p2 = fv.plot('SAL', ts, ax=ax, datum='sigma', limits=(0,0.2), **sal_spec)
ax.set_title('Bottom 20% of Water Column')

ax = axes[2]
p3 = fv.plot('SAL', ts, ax=ax, **sal_spec)
plt.colorbar(p3.patch, ax=ax, orientation='vertical',label='Salinity (PSU)')
ax.set_title('Depth Averaged')

fig.suptitle('Example TUFLOW FV Stratified Result', fontweight='bold', x=0.55)
titles = {0: 'Velocity', 1: 'Salinity', 2: 'Temperature'}

Save Sheet To ASCII Arc Grid

ts = '2011-05-02 00:00'
cell_size =  5*10**-5
data = fv.xtr.get_sheet_cell('V', ts, datum='depth', limits=(0,1))

# Output the ASCII to disk
output_folder = Path('./outputs')
output_folder.mkdir(exist_ok=True) # Make the folder, if it doesn't exist

fv.xtr.write_data_to_ascii(output_folder / 'Surface_V_2011-05-02.asc', data, cell_size)

Extract Grid From Sheet And Output To GeoTiff

This example shows how to write a sheet to a custom grid resolution and extent using the get_sheet_grid method. A raster is then output to geotiff using the rioxarray package.

The get_sheet_method can be used to convert your TUFLOW FV mesh results to a regular grid. This may be useful if you need to compare TUFLOW FV results from different computational meshes, or if you need to compare a TUFLOW FV result with a data raster. For example, a DEM of surveyed bathmetry verses TUFLOW FV bed elevation outputs.

RioXarray is a wrapper around the rasterio package that allows for easy rasterio integration with xarray. It is a very powerful package that allows for the easy export of xarray data to a variety of formats.

Caution: The example below is based on a spherical model that uses latitude and longitude, thus the dx/dy parameters are very small representing degrees. If your model is in meters they should be changed to something more appropriate (e.g., 20m dx=20, dy=20, 200m dx=200,dy=200 and so on).

import rioxarray

# Define bounding box [xmin, ymin, xmax, ymax]
bbox = [159.08597, -31.39684, 159.09753, -31.3877]

# Extract the TUFLOW FV results to a regular grid. This example extracts a single timestep, however the default it for all timesteps. 
grd = fv.get_sheet_grid(time='2011-05-02 00:00',dx=0.00005, dy=0.00005, bbox=bbox)

# Ouptut salinity to geotiff
output_folder = Path('./outputs')
output_folder.mkdir(exist_ok=True) # Make the folder, if it doesn't exist

grd['SAL'].rio.to_raster(output_folder / 'raster_geotiff.tif')

Multiple Sheet Plots

%matplotlib inline
# With datum='depth', this is the surface 1m. 
# Experiment with datum='height' to see how the results change. 
datum = 'depth'
limits = (0,1)

fig, axes = plt.subplots(ncols=3, figsize=(12, 6), constrained_layout=True)

ts = '2011-05-02 00:00'

ax = axes[0]
p1 = fv.plot('V', ts, ax=ax,  shading='interp', cmap='turbo', clim=(0, 0.25), datum=datum, limits=limits, colorbar=False)  # By setting datum to depth, this is now the "surface 1m")
v1 = fv.plot_vector('V', ts, ax=ax, color='k')   # We can plot a vector overlay using the `plot_vector` method here
ax.set_title('Velocity')
plt.colorbar(p1.patch, ax=ax, orientation='horizontal', label='Speed (m/s)')

ax = axes[1]
p2 = fv.plot('SAL', ts, ax=ax,  shading='interp', cmap='ocean', clim=(0, 36), datum=datum, limits=limits, colorbar=False)
ax.set_title('Salinity')
plt.colorbar(p2.patch, ax=ax, orientation='horizontal', label='Salinity (PSU)')

ax = axes[2]
p3 = fv.plot('TEMP', ts, ax=ax, shading='interp', cmap='plasma', clim=(5, 25), datum=datum, limits=limits, colorbar=False)
ax.set_title('Temperature')
plt.colorbar(p3.patch, ax=ax, orientation='horizontal', label='Temp (°C)')

plt.show()

Impact Sheet Multiple Models

# Open second set of results. This model has a 40 degree celcius thermal discharge in the surface 1m. 
model_file = 'HYD_002'
model_file_therm = 'HYD_002_Thermal_Discharge.nc'  

# Open file
fv_therm = xr.open_dataset(model_folder / model_file_therm, decode_times=False).tfv

# Subtracting based on the 3D RAW DATA
fv_therm['IMPACT'] = fv_therm['TEMP'] - fv['TEMP']

ts = '2011-05-02 00:00'
fig, axes = plt.subplots(ncols=3, figsize=(12, 6), constrained_layout=True)

temp_spec = dict(cmap='jet', clim=(15, 25), shading='interp', colorbar=False)
diff_spec = dict(cmap='jet', clim=(-2, 2), shading='interp', colorbar=False)

ax = axes[0]
p1 = fv.plot('TEMP', ts, ax=ax,  datum='depth', limits=(0,1), **temp_spec)  # By setting datum to depth, this is now the "surface 1m")
ax.set_title('Existing Conditions')

ax = axes[1]
p2 = fv_therm.plot('TEMP', ts, ax=ax, datum='depth', limits=(0,1.0), **temp_spec)
ax.set_title('Developed Conditions')

ax = axes[2]
p3 = fv_therm.plot('IMPACT', ts, ax=ax, datum='depth', limits=(0,1.0), **diff_spec)
plt.colorbar(p3.patch, ax=ax, orientation='vertical',label='Temperature Difference (°C)')
ax.set_title('Difference - Developed minus Existing')

fig.suptitle('Impact of Thermal Discharge at ' + ts, fontweight='bold', x=0.45)
titles = {0: 'Existing Conditions', 1: 'Developed Conditions', 2: 'Impact'}

plt.show()

Extract Contours From Sheet

This example shows how you can extract contours from a TUFLOW FV result as an array or a vector dataset. This is can be useful for highlighting thresholds in your model, or finding the area of a plume size, and so on.

The default way of using this function will return a dictionary containing keys corresponding to the requested levels, and values that are a list of numpy arrays containing the contour coordinates (Nx2 with X and Y columns).

We’ll follow on from the impact plot sample shown above. We will extract contours of the temperature difference above 0.5, 0.75 and 1.0°C.

levels = [0.5, 0.75, 1.0]  # temperature difference in degrees

gdf = fv_therm.get_contours(levels, 'IMPACT', ts, datum='depth', limits=(0, 1.0))

# contours is the default output - a dictionary 
# gdf is the optional output, a geodataframe. 
gdf

	level	geometry
0	0.50	POLYGON ((159.08624 -31.37594, 159.08642 -31.3...
1	0.75	POLYGON ((159.08819 -31.38019, 159.08831 -31.3...
2	1.00	POLYGON ((159.09403 -31.39216, 159.09396 -31.3...

fig, axes = plt.subplots(ncols=2, figsize=(10, 6), constrained_layout=True, sharey=True)
for ax in axes:
    ax.set_aspect('equal')

# Plot 1: The impacts from TUFLOW FV
fv_therm.plot('IMPACT', time=ts, ax=axes[0], datum='depth', limits=(0, 1), boundary=True, **diff_spec)
axes[0].set_title('Impacts From TUFLOW FV')

# Plot 2: The contoured up impact zones above our arbitary thresholds
fv_therm.plot('IMPACT', ax=axes[1],alpha=0.0, boundary=True, **diff_spec) # We're using this for the boundary line only :) 
gdf.plot(ax=axes[1], column='level', cmap='Reds', lw=1)
axes[1].set_title('Contoured Impacts Above Our Thresholds')

plt.show()

We can also export a GeoDataFrame to shapefile, which can be useful for further analysis, or for plotting in external programs like QGIS.

By default, the geodataframe we’ve made doesn’t know what coordinate system it is in. You can add on a coordinate system (to make it load in QGIS/ArcGIS easier) by using gdf = gdf.set_crs(<epsg_code>) as below.

gdf = gdf.set_crs(4326)
gdf.to_file(output_folder / 'Impact_Contours.shp')

Plot Sheet Max

fv.get_statistics('max', 'SAL', datum='sigma').plot()
plt.title('Maximum Salinity from Map Output')

Text(0.5, 1.0, 'Maximum Salinity from Map Output')

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

Save Max to asc

import numpy as np

stat_max = fv.get_statistics('max', 'SAL', datum='sigma')
cell_size =  5*10**-5

# Output the ASCII to disk
output_folder = Path('./outputs')
output_folder.mkdir(exist_ok=True) # Make the folder, if it doesn't exist

fv.xtr.write_data_to_ascii(output_folder / 'Max_Depth_Ave_Salinity.asc', np.squeeze(stat_max['SAL_max']), cell_size)

Save Sheet Plot Animation

To save an animation it is required to have the FFmpeg program on your computer’s Path environment variable.

# Output the ASCII to disk
output_folder = Path('./outputs')
output_folder.mkdir(exist_ok=True) # Make the folder, if it doesn't exist

out_file = output_folder / 'Salinty_Animation.mp4'
from matplotlib import animation

time_start = '2011-05-02 00:00'
time_end = '2011-05-02 12:00'

fvx = fv.sel(Time=slice(time_start, time_end))

# Colour settings
col_spec = dict(cmap='ocean', clim=(10, 36), shading='interp', edgecolor='face')

# Vector settings
vec_spec = dict(angles='uv', scale_units='dots', scale=1/200, width=0.002,
                headwidth=2.5, headlength=4, pivot='middle', minlength=2)

variable = 'SAL'
vector_variable = 'V'

fig, ax = plt.subplots()
ax.set_aspect('equal', adjustable='datalim')

sheet = fvx.plot(variable, 1, ax=ax, **col_spec) 
vec = fvx.plot_vector(vector_variable, 1, ax=ax, **vec_spec) 

def animate(i):
    sheet.set_time_current(i)
    vec.set_time_current(i)
    date = sheet.get_time_current()
    ax.set_title(date.strftime('%Y-%m-%d %H:%M'))

nframes = fvx.sizes['Time']
anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat=False)
anim.save(
    out_file, 
    writer=animation.FFMpegWriter(),
    dpi=300)

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

Interactive Sheet Plot

It is a common task to view your results interactively so that you can get a feel for how a model is responding and behaving.

While this can be handled very well in QGIS with the TUFLOW Viewer plugin, it can also be done directly through this notebook interface.

The syntax for the interactive viewer is generally same as using fv.plot. Any keyword arguments that can be used on fv.plot(..) can also be used here, such as datum, cmap, and so on.

Vectors can also be added onto this figure. The vectors keyword argument will enable velocity if True is passed through (and V_x and V_y are present), however other variables can be plotted by specifying them directly as a list (e.g., vectors=['SEDLOAD_TOTAL_x', 'SEDLOAD_TOTAL_y']). Keyword arguments for the plot_vector command can also be used here, passed through as a dictionary to vector_kwargs - for example, vector_kwargs=dict(scale=10, color='red')

%matplotlib widget

Note that interactive plots don’t work well in the documentation so we’ve copied in a screenshot instead

fv.plot_interactive('V', shading='interp', cmap='turbo', clim=(0,0.5), datum='depth', vectors=True)

To turn off interactive plotting (which can often make your sessions run faster, a good idea when you don’t need it), use %matplotlib inline

plt.close('all')
%matplotlib inline

This concludes the examples on sheet plotting.