Interactive application¶
Note
In this notebook, we explore some interactive plotting of gospl
outputs that might be useful for sharing the results within a jupyter environment.
We will create a Dashboard
with two HoloViews
objects:
a
panel
pn.widgets.Select
object that contains a list ofXarray
variables, anda
hvPlot
object that takes the selected variable on input.
See also
An in-depth description of the approach quickly presented here is well discussed in a recent paper by Signell & Pothina (2019)1.
Load the required Python libraries¶
First of all, load the necessary libraries. These are the ones we discussed previously:
numpy
matplotlib
cartopy
panel
xarray
holoviews
geoviews
import os
import pooch
import numpy as np
import xarray as xr
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))
import cmocean
import hvplot.xarray
import panel as pn
import holoviews as hv
from holoviews import opts, dim
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs
import geoviews as gv
from geoviews import tile_sources as gvts
import geoviews.feature as gf
from cartopy import crs as ccrs
import pandas as pd
from script import readOutput as rout
gv.extension('bokeh')
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
We first define a folder where exported files will be stored:
out_path = 'export'
if not os.path.exists(out_path):
os.makedirs(out_path)
Netcdf outputs¶
We start with the netcdf
outputs as it is the most common one used in our field.
Building netcdf file¶
Netcdf
exports are done by using the readOutput.py
script presented in the previous notebook. Here we export all the time steps at once by looping through the number of outputs (5 in this case).
Arguments for readOutput.py
The readOutput.py
script main function requires several arguments:
path
: the path to the input filefilename
: the name of the input filestep
: the step you wish to output (here set to 5 corresponding to the last output based on the input parameters: start time 0 year, end time 50 thousand years with an output every 10 thousand years)nbstep
: the number of time steps to plot (useful if one want to output anetdcf
file containing all time steps (done in the following section).uplift_forcing
: set to False as we are not considering any tectonic forcing
Then the buildLonLatMesh
function is used to interpolate (using a kd-tree
approach) the gospl
variables on a regular mesh. It also provides a way to limit the created netcdf
file by defining a bounding box
:
uplift_forcing = False
# Specifying the grid resolution in degrees
reso = 0.1
# Total number of outputs
nbstep = 5
# Bounding box
bb = [-134,20,-46,80]
def build_ncgrid(ncout, nbstep, bbox, reso, fname, uplift_forcing):
# Looping through the output time steps
for k in range(nbstep+1):
if k == 0:
# Calling the initialisation function for our class
ncgrid = rout.readOutput(path='./', filename=fname,
step=k, nbstep=nbstep+1,
uplift=uplift_forcing)
else:
# Update the variables after the first time steps
ncgrid.update(step=k, uplift=uplift_forcing)
# Build the regular grid defining the bounding box
ncgrid.buildLonLatMesh(res=reso, nghb=3, box=bbox)
# Export as netcdf
ncgrid.exportNetCDF(ncfile = ncout)
return
ncout = os.path.join(out_path, "GoMresult.nc")
# We commented the next line and loaded the dataset from figshare (see below)
#build_ncgrid(ncout, nbstep, bb, reso, fname='input.yml', uplift_forcing=False)
Here we will download the file directly from figshare:
ncout = pooch.retrieve(
url="https://ndownloader.figshare.com/files/27579284",
known_hash="sha256:fdc4a4ac70bca66ed20c3a126fdbf8fd2b33a2b47170e45c53eaf96f20f42a32",
downloader=pooch.HTTPDownloader(progressbar=True),
fname="GoMresult.nc",
)
Open the netcdf
file with xarray
file:
ds = xr.open_dataset(ncout, decode_times=False)
ds
<xarray.Dataset> Dimensions: (latitude: 681, longitude: 881, time: 6) Coordinates: * time (time) float64 0.0 1e+04 2e+04 3e+04 4e+04 5e+04 * latitude (latitude) float64 12.0 12.1 12.2 12.3 ... 79.8 79.9 80.0 * longitude (longitude) float64 -134.0 -133.9 -133.8 ... -46.1 -46.0 Data variables: elevation (time, latitude, longitude) float64 ... erodep (time, latitude, longitude) float64 ... precipitation (time, latitude, longitude) float64 ... drainageArea (time, latitude, longitude) float64 ... basinID (time, latitude, longitude) int32 ... Attributes: description: gospl outputs history: Created Thu Apr 8 17:39:07 2021
- latitude: 681
- longitude: 881
- time: 6
- time(time)float640.0 1e+04 2e+04 3e+04 4e+04 5e+04
- units :
- days since 0000-01-01
- calendar :
- 365_day
array([ 0., 10000., 20000., 30000., 40000., 50000.])
- latitude(latitude)float6412.0 12.1 12.2 ... 79.8 79.9 80.0
- units :
- degrees_north
array([12. , 12.1, 12.2, ..., 79.8, 79.9, 80. ])
- longitude(longitude)float64-134.0 -133.9 ... -46.1 -46.0
- units :
- degrees_east
array([-134. , -133.9, -133.8, ..., -46.2, -46.1, -46. ])
- elevation(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- erodep(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- precipitation(time, latitude, longitude)float64...
- units :
- m/yr
[3599766 values with dtype=float64]
- drainageArea(time, latitude, longitude)float64...
- units :
- m2
[3599766 values with dtype=float64]
- basinID(time, latitude, longitude)int32...
- units :
- int
[3599766 values with dtype=int32]
- description :
- gospl outputs
- history :
- Created Thu Apr 8 17:39:07 2021
Plotting a specific variable¶
We will plot the basin index for a specific region.
Let’s start by clipping the area to reduce the Dataset size. We will clip the spatial extent based on longitudinal and latitudinal values.
This is done using the sel
function with the slice
method.
min_lon = -130 # lower left longitude
min_lat = 21 # lower left latitude
max_lon = -55 # upper right longitude
max_lat = 70 # upper right latitude
# Defining the boundaries
lon_bnds = [min_lon, max_lon]
lat_bnds = [min_lat, max_lat]
# Time interval
tsteps = [1.e04, 5.e04]
# Performing the reduction
ds_clip = ds.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds), time=slice(*tsteps))
ds_clip
<xarray.Dataset> Dimensions: (latitude: 491, longitude: 751, time: 5) Coordinates: * time (time) float64 1e+04 2e+04 3e+04 4e+04 5e+04 * latitude (latitude) float64 21.0 21.1 21.2 21.3 ... 69.8 69.9 70.0 * longitude (longitude) float64 -130.0 -129.9 -129.8 ... -55.1 -55.0 Data variables: elevation (time, latitude, longitude) float64 ... erodep (time, latitude, longitude) float64 ... precipitation (time, latitude, longitude) float64 ... drainageArea (time, latitude, longitude) float64 ... basinID (time, latitude, longitude) int32 ... Attributes: description: gospl outputs history: Created Thu Apr 8 17:39:07 2021
- latitude: 491
- longitude: 751
- time: 5
- time(time)float641e+04 2e+04 3e+04 4e+04 5e+04
- units :
- days since 0000-01-01
- calendar :
- 365_day
array([10000., 20000., 30000., 40000., 50000.])
- latitude(latitude)float6421.0 21.1 21.2 ... 69.8 69.9 70.0
- units :
- degrees_north
array([21. , 21.1, 21.2, ..., 69.8, 69.9, 70. ])
- longitude(longitude)float64-130.0 -129.9 ... -55.1 -55.0
- units :
- degrees_east
array([-130. , -129.9, -129.8, ..., -55.2, -55.1, -55. ])
- elevation(time, latitude, longitude)float64...
- units :
- metres
[1843705 values with dtype=float64]
- erodep(time, latitude, longitude)float64...
- units :
- metres
[1843705 values with dtype=float64]
- precipitation(time, latitude, longitude)float64...
- units :
- m/yr
[1843705 values with dtype=float64]
- drainageArea(time, latitude, longitude)float64...
- units :
- m2
[1843705 values with dtype=float64]
- basinID(time, latitude, longitude)int32...
- units :
- int
[1843705 values with dtype=int32]
- description :
- gospl outputs
- history :
- Created Thu Apr 8 17:39:07 2021
Using GeoViews
we visualise the basin index over time:
# Specify the dataset, its coordinates and requested variable
dataset = gv.Dataset(ds_clip, ['longitude', 'latitude', 'time'],
'basinID', crs=crs.PlateCarree())
images = dataset.to(gv.Image,dynamic=True)
# Loading coastlines from Cartopy
coastline = gf.coastline(line_width=1,line_color='k').opts(projection=ccrs.PlateCarree(),
scale='10m')
# Slider location
hv.output(widget_location='bottom')
# Create stack of images grouped by time
images.opts(active_tools=['wheel_zoom', 'pan'], cmap='jet',
colorbar=True, width=800, height=500,
) * coastline
/usr/share/miniconda/envs/stellar/lib/python3.8/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/110m/physical/ne_110m_coastline.zip
warnings.warn('Downloading: {}'.format(url), DownloadWarning)
/usr/share/miniconda/envs/stellar/lib/python3.8/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/10m/physical/ne_10m_coastline.zip
warnings.warn('Downloading: {}'.format(url), DownloadWarning)