%pylab inline
Populating the interactive namespace from numpy and matplotlib
"Validation class" is a hybrid class which requires both a numerical-model-based object (e.g. "FVCOM" or "Station" object) and at least one measurement-based object (e.g. "ADCP", "TideGauge" or "Drifter" object or a list of their combination). Its aim is to provide validation benchmarks in order to access the accuracy of a FVCOM simulation compared to matching measurements.
As any other library in Python, PySeidon-dvt has to be first imported before to be used. Here we will use an alternative import statement compared to the one previoulsy presented:
from pyseidon_dvt import *
/home/grumpynounours/anaconda/lib/python2.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`. "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)
*** Star *** here means all. Usually this form of statements would import the entire library. In the case of PySeidon-dvt, this statement will import the following object classes: FVCOM, Station, Validation, ADCP, Tidegauge and Drifter. Only the Validation class will be tackle in this tutorial.
Python is by definition an object oriented language...and so is matlab. PySeidon-dvt is based on this notion of object, so let us define our first "Validation" object.
*Exercise 1: *
*Answer: *
Validation?
According to the documentation, in order to define a Validation object, the required inputs are:
Note:
*Exercise 2: *
*Answer: *
drift = Drifter('./data4tutorial/drifter_GP_01aug2013.mat')
adcp = ADCP('./data4tutorial/adcp_GP_01aug2013.mat')
tg = TideGauge('./data4tutorial/tidegauge_GP_01aug2013.mat')
station = Station('http://ecoii.acadiau.ca/thredds/dodsC/ecoii/test/Station3D_dngrid_BF_20130730_20130809.nc')
fvcom = FVCOM('http://ecoii.acadiau.ca/thredds/dodsC/ecoii/test/FVCOM3D_dngrid_GP_20130801_20130802.nc')
val1 = Validation(adcp, station, flow=10.0)
val2 = Validation([adcp, drift, tg], fvcom, flow='daf')
val3 = Validation(adcp, fvcom, flow=-5.0)
Retrieving data from /media/grumpynounours/Seagate Backup Plus Drive/FVCOM_simulations/sim_tuto/Station3D_dngrid_BF_20130730_20130809.nc ... Initialisation... Retrieving data from /media/grumpynounours/Seagate Backup Plus Drive/FVCOM_simulations/sim_tuto/FVCOM3D_dngrid_GP_20130801_20130802.nc ... Initialisation... Station site: GP_130730_TA
The Validation object possesses 2 attributes and 5 methods. They would appear by typing *val1. Tab* for instance.
An attribute is a quantity intrinsic to its object. A method is an intrinsic function which changes an attribute of its object. Contrarily a function will generate its own output:
The Validation attributes are:
The Validation methods & functions are:
*Exercise 3: *
*Answer: *
val1.validate_data(plot=True, save_csv=True, filename='test1')
val1.taylor_diagram()
val1.validate_harmonics(filename='test2')
/home/grumpynounours/anaconda/lib/python2.7/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice warnings.warn("Mean of empty slice", RuntimeWarning) /home/grumpynounours/anaconda/lib/python2.7/site-packages/matplotlib/figure.py:397: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure "matplotlib is currently using a non-GUI backend, "
---Validation benchmarks--- MDPO pbias gear r2 RMSE NOF CF NSE \ adcp_GP_01aug2013 0.0 -311.17 ADCP 1.00 0.16 0.00 99.29 0.98 adcp_GP_01aug2013 100.0 85950.52 ADCP 0.96 2.42 20.57 56.03 0.61 adcp_GP_01aug2013 340.0 162.73 ADCP 0.89 0.24 12.06 17.02 -0.89 adcp_GP_01aug2013 0.0 -109.28 ADCP 0.99 0.26 0.00 90.78 0.97 adcp_GP_01aug2013 0.0 -84.77 ADCP 0.99 0.35 0.00 86.52 0.96 adcp_GP_01aug2013 30.0 -14.39 ADCP 0.91 0.27 14.18 36.17 0.72 adcp_GP_01aug2013 70.0 351.18 ADCP 0.89 6.20 2.84 87.94 0.86 phase SI bias NRMSE corr POF NMSE skill \ adcp_GP_01aug2013 0.0 -3.54 0.14 12.47 1.00 0.00 1.56 1.00 adcp_GP_01aug2013 0.0 -2744.23 -0.76 62.15 0.98 1.42 38.63 0.94 adcp_GP_01aug2013 0.0 4.72 0.08 132.11 0.94 39.01 174.53 0.81 adcp_GP_01aug2013 0.0 3.10 -0.09 18.64 0.99 0.00 3.48 0.99 adcp_GP_01aug2013 0.0 3.91 -0.07 24.33 0.99 0.00 3.95 0.99 adcp_GP_01aug2013 0.0 0.21 -0.19 19.36 0.96 0.71 3.75 0.94 adcp_GP_01aug2013 0.0 14.10 7.93 6.89 0.94 3.55 13.53 0.97 Type MSE MDNO SD adcp_GP_01aug2013 el 0.02 0.0 0.35 adcp_GP_01aug2013 cubic_speed 5.85 210.0 1.40 adcp_GP_01aug2013 u 0.06 250.0 0.45 adcp_GP_01aug2013 v 0.07 40.0 0.48 adcp_GP_01aug2013 vel 0.08 70.0 0.50 adcp_GP_01aug2013 speed 0.07 290.0 0.52 adcp_GP_01aug2013 dir 1013.62 30.0 7.71
/home/grumpynounours/anaconda/lib/python2.7/site-packages/matplotlib/lines.py:1106: UnicodeWarning: Unicode unequal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if self._markerfacecolor != fc: /home/grumpynounours/anaconda/lib/python2.7/site-packages/matplotlib/lines.py:1082: UnicodeWarning: Unicode unequal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if self._markeredgecolor != ec:
solve:
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-5-65c6c7a24794> in <module>() 1 val1.validate_data(plot=True, save_csv=True, filename='test1') 2 val1.taylor_diagram() ----> 3 val1.validate_harmonics(filename='test2') /home/grumpynounours/Desktop/Github/PySeidon_dvt/pyseidon_dvt/validationClass/validationClass.pyc in validate_harmonics(self, filename, save_csv, debug, debug_plot) 605 # Make directory 606 self._make_save_folder() --> 607 self._validate_harmonics(filename, save_csv, debug, debug_plot) 608 self.HarmonicBenchmarks = self._HarmonicBenchmarks 609 else: /home/grumpynounours/Desktop/Github/PySeidon_dvt/pyseidon_dvt/validationClass/validationClass.pyc in _validate_harmonics(self, filename, save_csv, debug, debug_plot) 290 self.Variables.obs.elCoef = solve(obs_time, obs_el, None, obs_lat, 291 constit='auto', trend=False, Rayleigh_min=0.95, --> 292 method='ols', conf_int='linear') 293 self.Variables.sim.elCoef = solve(mod_time, mod_el, None, mod_lat, 294 constit='auto', trend=False, Rayleigh_min=0.95, /home/grumpynounours/anaconda/lib/python2.7/site-packages/UTide-0.1b0.dev0-py2.7.egg/utide/_solve.pyc in solve(tin, uin, vin, lat, **opts) 44 ''' 45 ---> 46 coef = _solv1(tin, uin, vin, lat, **opts) 47 48 return coef /home/grumpynounours/anaconda/lib/python2.7/site-packages/UTide-0.1b0.dev0-py2.7.egg/utide/_solve.pyc in _solv1(tin, uin, vin, lat, **opts) 52 53 print('solve: ') ---> 54 packed = _slvinit(tin, uin, vin, lat, **opts) 55 t, u, v, tref, lor, elor, opt, tgd, uvgd = packed 56 nt = len(t) /home/grumpynounours/anaconda/lib/python2.7/site-packages/UTide-0.1b0.dev0-py2.7.egg/utide/_solve.pyc in _slvinit(tin, uin, vin, lat, **opts) 260 opt['equi'] = False 261 lor = np.ptp(t) --> 262 elor = lor*nt / (nt-1) 263 tref = 0.5*(t[0]+t[-1]) 264 NameError: global name 'nt' is not defined
*Exercise 4: *
*Answer: *
val2.validate_data(save_csv=True, filename='test3')
val2.taylor_diagram()
---Validation benchmarks--- MDPO pbias gear r2 RMSE NOF CF NSE phase \ adcp_GP_01aug2013 0.0 -227.44 ADCP 1.00 0.12 0.00 99.12 0.99 0.0 adcp_GP_01aug2013 50.0 -39.91 ADCP 0.99 1.09 0.00 71.93 0.92 0.0 adcp_GP_01aug2013 340.0 -92.39 ADCP 0.95 0.21 47.37 10.53 -0.87 0.0 adcp_GP_01aug2013 0.0 -5.49 ADCP 1.00 0.10 0.00 98.25 0.99 0.0 adcp_GP_01aug2013 0.0 -7.63 ADCP 1.00 0.24 0.00 98.25 0.99 0.0 adcp_GP_01aug2013 10.0 -6.04 ADCP 0.98 0.12 0.00 85.96 0.94 0.0 adcp_GP_01aug2013 10.0 29.89 ADCP 0.99 6.54 0.88 95.61 0.98 0.0 SI bias NRMSE corr POF NMSE skill \ adcp_GP_01aug2013 -2.63 0.11 9.26 1.00 0.00 0.86 1.00 adcp_GP_01aug2013 1.19 -0.37 27.76 0.99 0.00 7.71 0.98 adcp_GP_01aug2013 3.74 -0.05 128.56 0.97 26.32 165.28 0.83 adcp_GP_01aug2013 0.23 -0.02 6.99 1.00 0.00 0.49 1.00 adcp_GP_01aug2013 0.54 -0.03 16.80 1.00 0.00 0.77 1.00 adcp_GP_01aug2013 0.09 -0.08 8.82 0.99 0.00 0.78 0.99 adcp_GP_01aug2013 0.50 7.06 7.27 0.99 0.00 1.81 1.00 Type MSE MDNO SD adcp_GP_01aug2013 el 0.01 0.0 0.31 adcp_GP_01aug2013 cubic_speed 1.19 150.0 0.93 adcp_GP_01aug2013 u 0.04 330.0 0.44 adcp_GP_01aug2013 v 0.01 0.0 0.29 adcp_GP_01aug2013 vel 0.02 0.0 0.33 adcp_GP_01aug2013 speed 0.02 30.0 0.33 adcp_GP_01aug2013 dir 138.00 10.0 6.54
/home/grumpynounours/anaconda/lib/python2.7/site-packages/numpy/core/_methods.py:59: RuntimeWarning: Mean of empty slice. warnings.warn("Mean of empty slice.", RuntimeWarning)
---Validation benchmarks--- MDPO pbias gear r2 RMSE NOF CF \ drifter_GP_01aug2013 2.000000 69.90 Drifter 0.71 3.56 27.27 39.39 drifter_GP_01aug2013 1.000000 13.66 Drifter 0.82 0.12 6.06 69.70 drifter_GP_01aug2013 10.999999 11.53 Drifter 0.62 0.25 0.00 30.30 drifter_GP_01aug2013 10.999999 11.23 Drifter 0.67 0.28 0.00 30.30 drifter_GP_01aug2013 10.999999 11.23 Drifter 0.67 0.25 0.00 30.30 drifter_GP_01aug2013 3.000000 -0.07 Drifter 0.78 11.14 6.06 51.52 NSE phase SI bias NRMSE corr POF NMSE \ drifter_GP_01aug2013 0.63 0.0 -1.64 -1.52 57.13 0.84 3.03 32.64 drifter_GP_01aug2013 0.81 0.0 -1.05 -0.02 40.36 0.91 6.06 16.29 drifter_GP_01aug2013 0.03 0.0 0.15 0.20 14.57 0.79 54.55 2.12 drifter_GP_01aug2013 0.15 0.0 0.16 0.20 15.98 0.82 54.55 2.02 drifter_GP_01aug2013 0.15 0.0 0.14 0.20 14.20 0.82 54.55 2.02 drifter_GP_01aug2013 0.71 0.0 0.05 -0.06 12.38 0.88 6.06 0.24 skill Type MSE MDNO SD drifter_GP_01aug2013 0.87 cubic_speed 12.70 9.999999 2.00 drifter_GP_01aug2013 0.95 u 0.02 3.000000 0.31 drifter_GP_01aug2013 0.75 v 0.06 1.000000 0.43 drifter_GP_01aug2013 0.77 vel 0.06 1.000000 0.43 drifter_GP_01aug2013 0.77 speed 0.06 1.000000 0.43 drifter_GP_01aug2013 0.93 dir 21.41 1.000000 1.93 ---Validation benchmarks--- MDPO pbias gear r2 RMSE NOF CF NSE \ tidegauge_GP_01aug2013 0 973.35 TideGauge 1.0 0.12 0.0 100.0 0.99 phase SI bias NRMSE corr POF NMSE skill Type \ tidegauge_GP_01aug2013 0.0 10.3 0.12 9.24 1.0 0.0 0.85 1.0 el MSE MDNO SD tidegauge_GP_01aug2013 0.02 0 0.32
*Exercise 5: *
val3.validate_data(save_csv=True, phase_shift=True, filename='test4')
val3.benchmarks_map()
val3.write_validation_report(report_title="val3_report.pdf")
---Validation benchmarks--- MDPO pbias gear r2 RMSE NOF CF NSE \ adcp_GP_01aug2013 0.0 -227.44 ADCP 1.00 0.12 0.00 99.12 0.99 adcp_GP_01aug2013 250.0 69.58 ADCP 0.98 1.59 0.00 54.39 0.88 adcp_GP_01aug2013 130.0 9517.82 ADCP 0.82 0.31 60.53 14.04 -4.91 adcp_GP_01aug2013 100.0 41.05 ADCP 1.00 0.28 0.00 71.05 0.96 adcp_GP_01aug2013 0.0 32.16 ADCP 1.00 0.43 0.00 98.25 0.98 adcp_GP_01aug2013 300.0 12.07 ADCP 0.93 0.22 0.00 42.98 0.84 adcp_GP_01aug2013 40.0 57.73 ADCP 0.97 9.18 0.00 85.09 0.93 phase SI bias NRMSE corr POF NMSE skill \ adcp_GP_01aug2013 0.0 -2.63 0.11 9.26 1.00 0.00 0.86 1.00 adcp_GP_01aug2013 0.0 1.09 1.02 33.17 0.99 5.26 11.00 0.96 adcp_GP_01aug2013 0.0 -140.94 -0.21 243.03 0.91 7.02 590.64 0.62 adcp_GP_01aug2013 0.0 0.55 0.21 18.76 1.00 0.00 3.52 0.99 adcp_GP_01aug2013 0.0 0.83 0.17 28.14 1.00 0.00 2.14 0.99 adcp_GP_01aug2013 0.0 0.16 0.17 14.66 0.97 0.00 2.15 0.96 adcp_GP_01aug2013 0.0 0.81 15.90 10.20 0.99 3.51 6.09 0.98 Type MSE MDNO SD adcp_GP_01aug2013 el 0.01 0.0 0.31 adcp_GP_01aug2013 cubic_speed 2.54 20.0 1.05 adcp_GP_01aug2013 u 0.10 390.0 0.54 adcp_GP_01aug2013 v 0.08 0.0 0.47 adcp_GP_01aug2013 vel 0.05 0.0 0.42 adcp_GP_01aug2013 speed 0.05 10.0 0.42 adcp_GP_01aug2013 dir 498.76 0.0 15.39
/home/grumpynounours/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py:1215: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if aspect == 'normal': /home/grumpynounours/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py:1220: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal elif aspect in ('equal', 'auto'):
...plotting adcp_GP_01aug2013...
<matplotlib.figure.Figure at 0x7f209e38fd10>
This collection and analysis of a set of statistics mostly adhere to the benchmarks defined as standards for hydrodynamic model validation by NOAA [1]. Additional statistics have been added to provide additional clarity on the skill of the model [2, 3, 4].
The present validation set is usually performed the following variables yet could be extended to other hydrodynamic quantities:
*Benchmarks*
Following is a list of the statistics used to evaluate model skil
The value of X is dependent on the variable and the statistic it's being used for. As defined here, X equals 10% of the data range. NOAA additionally defines limits on the values for determining the skill of the model. Accordingly, CF can be no lower than 90%, NOF/POF can be no higher than 1%, and MDPO/MDNO can be no longer than 24hr/1440min. ideally, the Willmott Skill and R2 should be 1, and the phase should be 0.
*References*
[1] K. W. Hess, T. F. Gross, R. A. Schmalz, J. G. Kelley, F. Aikman and E. Wei, "NOS standards for evaluating operational nowcast and forecst hydrodynamic model systems," National Oceanic and Atmospheric Administration, Silver Srping, Maryland, 2003.
[2] K. Gunn and C. Stock-Williams, "On validating numerical hydrodynamic models of complex tidal flow," International Journal of Marine Energy, Vols. 3-4, no. Special, pp. 82-97, 2013.
[3] N. Georgas and A. F. Blumberg, "Establishing Confidence in Marine Forecast Systems: The Design and Skill Assessment of the New York Harbor Observation and Prediction System, Version 3 (NYHOPS v3)," in 11th International Conference on Estuarine and Coastal Modeling, Seattle, Washington, United States, 2010.
[4] Y. Liu, P. MacCready, H. M. Barbara, E. P. Dever, M. Kosro and N. S. Banas, "Evaluation of a coastal ocean circulation model for the Columbia River plume in summer 2004," Journal of Geophysical Research, vol. 114, no. C2, p. 1978–2012, 2009.
As beta tester, your first assignement is to report bugs...yet not everything is a bug. The first thing to check before to report a bug is to verify that your version of PySeidon-dvt is up-to-date. The best way to keep up with the package evolution is to *git* to *clone* the repository, use *pull* to update it and *re-install* it if needed.
The second thing to check before to report a bug is to verify that the bug is *reproducible. When running into a bug, double check that your inputs fit the description of the documentation then turn the debug flag on* (e.g. output = tidegaugeobject.function(inputs, debug=True)) and submit the command again. If the error re-occurs then report it (i.e. copy entire error message + command and send it to package administrator)
Your second role as beta-tester is to submit suggestions and critics to the developpers regarding the functioning and functionality of the package. Beta testing phase is the best opportunity to steer a project towards the applications you would like to be tackled...