Commit b45e6252 authored by Pierre Cazenave's avatar Pierre Cazenave

Add script to extract a variable of interest and plot it at a number of sites...

Add script to extract a variable of interest and plot it at a number of sites from a number of different model runs.
parent fdfbee92
"""
Compare the modelled temperature time series for n models against the CEFAS
time series.
The model grids must match exactly in the different model runs.
"""
import os
import glob
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.ticker import ScalarFormatter
from matplotlib.dates import DateFormatter
from PyFVCOM.read_FVCOM_results import ncread
from PyFVCOM.grid_tools import findNearestPoint
def timeseries(conf, sites, noisy=False):
"""
Child function which, given a dict of configuration parameters, loads the
FVCOM results and extracts the time series at the specified locations.
Results are returned as a dict whose keys match those in the conf input
dict.
Parameters
----------
conf : dict
Configuration dict of dicts. First key is the overall model run type,
second set are:
base : str
Where the files are stored
vars : list, tuple
Variable names to extract from the netCDF. Ignored if using
probes output.
dims : dict
Dict of the variable dimensions to use (e.g. {'time':':100'}).
Ignored if using probes output.
project : str
Main project name (e.g. pml-irish-sea)
domain : str
Model domain name (e.g. irish_sea)
version : str
Which grid version (e.g. v20)
casename : str
FVCOM casename (e.g. irish_sea_v20)
binname : str
Binary name used in model output file names
runname : str
Suffix appended to the year output directory
netCDF : bool
Whether to load from netCDF (True) or probes output (False)
sites : ndarray
Array of positions (lon/lat) for which to extract corresponding model
time series.
noisy : bool, optional
Set to True to enable verbose output (defaults to False).
Returns
-------
FVCOM : dict
Dict with the modelled time series for the buoy locations. Keys are the
same as in conf.
valid_indices : ndarray
Array of indices for the data in FVCOM at the sites specified.
"""
FVCOM, sites_clean = {}, {}
for p in conf:
base = conf[p]['base']
casename = conf[p]['casename']
runname = conf[p]['runname']
netCDF = conf[p]['netCDF']
dims = conf[p]['dims']
noisy = conf[p]['noisy']
if netCDF:
# Use the first model output to find the relevant indices and load those data.
fvcom = glob.glob(os.path.join(
conf[p]['base'], 'output', conf[p]['runname'],
'{0}_0001.nc'.format(conf[p]['domain']))
)
F = ncread(fvcom, vars=['lon', 'lat'])
_, _, _, timeseriesNodes = findNearestPoint(F['lon'], F['lat'], sites[:, 0], sites[:, 1])
timeseriesNodesClean = []
sites_clean = []
valid_indices = []
# Clear out the NaN values
for c, i in enumerate(timeseriesNodes):
if ~np.isnan(i):
timeseriesNodesClean.append(int(i))
sites_clean.append(sites[c])
valid_indices.append(c)
dims.update({'node':str(timeseriesNodesClean)})
F = ncread(fvcom, vars, dims=dims, noisy=noisy)
else:
F = {}
sites_clean = []
valid_indices = []
c = 0
for ii, s in enumerate(sites):
probes = glob.glob(
os.path.join(
base,
'output',
casename,
runname,
'probes',
'{}_t1.dat'.format(s.title().replace('Wavenet', 'WaveNet'))
)
)
if probes:
sites_clean.append(s)
valid_indices.append(c)
c += 1
if noisy:
print('Load probes output at {}'.format(
s.title().replace('Wavenet', 'WaveNet').replace('_', ' ')))
times, values = readProbes(probes, noisy=noisy)
if ii == 0:
F['time'] = times
F['temp'] = values[:, 0] # keep only surface data
else:
F['temp'] = np.column_stack((F['temp'], values[:, 0]))
F['temp'] = F['temp'].squeeze()
# Sort out the dates for plot_date().
F['times'] = F['time'].astype(np.float64) + 678575.0 # relative to 0001/01/01 00:00:00 + 1.
FVCOM[p] = F
return FVCOM, valid_indices
if __name__ == '__main__':
noisy = True
# Array of locations at which to extract time series (lon/lat).
sites = np.array((
(-78.1772, 43.9643),
(-78.2061, 43.9413),
(-78.2055, 43.9107),
(-78.0374, 43.9530),
(-77.8327, 43.9655),
(-77.8230, 44.0145),
(-77.8108, 43.9125)))
# What do we want to call the different model runs? These names are
# entirely arbitrary, but are used as keys for the dict which contains each
# model configuration and result.
names = ('typeIII', 'typeII', 'typeI', 'typeIA', 'type1')
#'kraus1972', 'paulsonandsimpson1977_composite4-5-6-9-10',
#'paulsonandsimpson1977_run1')
# Model inputs (define as many as you like)
dims = {'siglay':'0'}
vars = ['lat', 'lon', 'nv', 'h', 'time', 'Times', 'temp']
conf = {}
fvcom = []
for i, k in enumerate(names):
conf[k] = {}
conf[k]['project'] = 'Estuary'
conf[k]['domain'] = 'tst'
conf[k]['version'] = ''
conf[k]['casename'] = '{}_{}'.format(conf[k]['domain'], k)
conf[k]['binname'] = ''
conf[k]['runname'] = k
conf[k]['base'] = '/data/medusa/pica/models/FVCOM/fvcom-examples/{}/run/'.format(conf[k]['project'])
conf[k]['netCDF'] = True
conf[k]['dims'] = dims
conf[k]['vars'] = vars
conf[k]['noisy'] = noisy
# Get the output files.
fvcom.append(glob.glob(
os.path.join(
conf[k]['base'], 'output', conf[k]['runname'],
'{0}_0001.nc'.format(conf[k]['domain']))
))
# Load all the time series from the model output for the sites at which
# we've got observations.
F, valid_indices = timeseries(conf, sites)
# Fix number formatting.
date_fmt = DateFormatter('%m')
formatter = ScalarFormatter(useOffset=False)
# Plot a time series of the current temperature and salinity at the sites
# in timeseries.
colours = ['r', 'b', 'k', 'c', 'y']
styles = ['-', '--', '-.']
fig1 = plt.figure(figsize=(15, 10))
for i, c in enumerate(valid_indices):
ax1 = fig1.add_subplot(np.ceil(len(valid_indices) / 2.0), 2, i + 1)
# Add each set of model output.
for j, p in enumerate(names):
FVCOM = F[p]
# No need to use the timeseriesNodesClean value for the FVCOM index
# because we've only extracted the indices in timeseriesNodesClean i.e.
# FVCOM has only len(timeseriesNodesClean) positions in it.
if j > len(colours):
style = styles[j % len(styles)]
else:
style = styles[0]
T1 = ax1.plot_date(FVCOM['times'], FVCOM['temp'][:, i], '.-',
color=colours[j % len(colours)],
ls=style,
label=p.replace('type', 'Type '), zorder=200)
ax1.legend(frameon=False, loc=2,
ncol=len(names))
#ncol=np.ceil(len(names) / 2).astype(int))
ax1.set_xlim(np.min(FVCOM['times']) - 6, np.max(FVCOM['times']))
ax1.xaxis.set_major_formatter(date_fmt)
ax1.set_ylabel('Temperature ($^{\circ}C$)')
ax1.yaxis.set_major_formatter(formatter)
fig1.autofmt_xdate(rotation=0, ha='center')
fig1.tight_layout(pad=2)
fig1.show()
baseName = os.path.join('..',
'figures',
'fvcom_{}_timeseries_comparison'.format(
'_'.join(names)
)
)
#fig1.savefig('{}.png'.format(baseName), bbox_inches='tight', pad_inches=0.2, dpi=300)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment