added update to documentation for tide checker

thopri · thopri · commit 30b1faa31593 · 2020-05-14T11:59:43.000+01:00
diff --git a/docs/source/_static/comparision_fes.png b/docs/source/_static/comparision_fes.png
diff --git a/docs/source/tides.rst b/docs/source/tides.rst
@@ -40,6 +40,60 @@ not been tested with the global tide models.
 Other options include "ln_tide" a boolean that when set to true will generate tidal boundaries. "sn_tide_model" is a string that defines the model to use, currently only
 "fes" or "tpxo" are supported. "ln_trans" is a boolean that when set to true will interpolate transport rather than velocities.
 
+Harmonic Output Checker
+-----------------------
+
+There is an harmonic output checker that can be utilised to check the output of PyNEMO with a reference tide model. So far
+the only supported reference model is FES but TPXO will be added in the future. Any tidal output from PyNEMO can be checked
+(e.g. FES and TPXO). While using the same model used as input to check output doesn't improve accuracy, it does confirm that the
+output is within acceptable/expected limits of the nearest model reference point.
+
+There are differences as PyNEMO interpolates the harmonics and the tidal checker does not, so there can be some difference
+in the values particularly close to coastlines.
+
+The checker can be enabled by editing the following in the relevent bdy file::
+
+    ln_tide_checker = .true.                ! run tide checker on PyNEMO tide output
+    sn_ref_model    = 'fes'                 ! which model to check output against (FES only)
+
+The boolean determines if to run the checker or not, this takes place after creating the interpolated harmonics
+and writing them to disk. The string denotes which tide model to use as reference, so far only FES is supported.
+The string denoting model is not strictly needed, by default fes is used.
+
+The checker will output information regarding the checking to the NRCT log, and also write an spreadsheet to the output folder containing any
+exceedance values, the closest reference model value and their locations. Amplitude and phase are checked independently, so both have latitude and longitude
+associated with them. It is also useful to know the amplitude of a exceeded phase to see how much impact it will have so this
+is also written to the spreadsheet. An example output is shown below, as can be seen the majority of the amplitudes, both
+the two amplitudes exceedances and the ones associated with the phase exceedances are low (~0.01), so can most likely be ignored.
+There a few phase exceedances that have higher amplitudes (~0.2) which would potentially require further investigation. A common
+reason for such an exceedance is due to coastlines and the relevant point being further away from an FES data point.
+
+Tide Checker Example Output for M2 U currents
+---------------------------------------------
+
+.. figure:: _static/comparision_fes.png
+   :align:   center
+
+The actual thresholds for both amplitude and phase are based on the amplitude of the output or reference, this is due to
+different tolerances based on the amplitude. e.g. high amplitudes should have lower percentage differences to the FES reference,
+than lower ones simply due to the absolute amount of the ampltiude itself, e.g. a 0.1 m difference for a 1.0 m amplitude is
+acceptable but not for a 0.01 m amplitude. The smaller amplitudes contribute less to the overall tide height so larger percentage
+differences are acceptable. The same also applies to phases, where large amplitude phases have little room for differences but at
+lower amplitudes this is less critical so a higher threshold is tolerated.
+
+The following power functions are used to determine what threshold to apply based on the reference model amplitude.
+
+Amplitude Threshold
+-------------------
+
+.. important:: Percentage Exceedance = 26.933 * Reference Amplitude ^ -0.396'
+
+Phases Threshold
+----------------
+
+.. important:: Phase Exceedance = 5.052 * PyNEMO Amplitude ^ -0.60
+
+
 Future work
 -----------
 
diff --git a/inputs/namelist_cmems.bdy b/inputs/namelist_cmems.bdy
@@ -107,9 +107,9 @@
     sn_tide_model  = 'fes'                !  Name of tidal model (fes|tpxo)
     clname(1)      = 'M2'                 !  constituent name
     clname(2)      = 'S2'
-    clname(3)      = 'N2'
-    clname(4)      = 'O1'
-    clname(5)      = 'K1'
+    !clname(3)      = 'N2'
+    !clname(4)      = 'O1'
+    !clname(5)      = 'K1'
     !clname(6)      = 'K2'
     !clname(7)      = 'L2'
     !clname(8)      = 'NU2'
@@ -123,7 +123,6 @@
     ln_trans       = .false.                !  interpolate transport rather than velocities
     ln_tide_checker = .true.                ! run tide checker on PyNEMO tide output
     sn_ref_model    = 'fes'                 ! which model to check output against (FES only)
-    nn_amp_thres    = 0.10                  ! amplitude thresold to compare against (m)
 !------------------------------------------------------------------------------
 !  Time information
 !------------------------------------------------------------------------------
diff --git a/pynemo/profile.py b/pynemo/profile.py
@@ -443,7 +443,7 @@ def process_bdy(setup_filepath=0, mask_gui=False):
 
             if settings['tide_checker'] == True:
                 logger.info('tide checker starting now.....')
-                tt_test = tt.main(setup_filepath,settings['amp_thres'],settings['ref_model'])
+                tt_test = tt.main(setup_filepath,settings['ref_model'])
                 if tt_test == 0:
                     logger.info('tide checker ran successfully, check spreadsheet in output folder')
                 if tt_test !=0:
@@ -457,7 +457,7 @@ def process_bdy(setup_filepath=0, mask_gui=False):
 
             if settings['tide_checker'] == True:
                 logger.info('tide checker starting now.....')
-                tt_test = tt.main(setup_filepath,settings['amp_thres'],settings['ref_model'])
+                tt_test = tt.main(setup_filepath,settings['ref_model'])
                 if tt_test == 0:
                     logger.info('tide checker ran successfully, check spreadsheet in output folder')
                 if tt_test !=0:
diff --git a/pynemo/tests/nemo_tide_test.py b/pynemo/tests/nemo_tide_test.py
@@ -32,6 +32,15 @@
 a separate sheet in the spreadsheet. The name of the spreadsheet contains meta data showing thresholds and reference
 model used. Units for threshold are meters and degrees.
 
+Update: fill values for FES are commonly returned at coastlines, this is due to the nearest FES cell being land but PyNEMO
+will have interpolated data from the water. In instance the code checks the cells aroud the fill value and averages both
+amplitude and phase (using HsinG,HcosG) to act as a reference.
+
+Phase threshold is not longer required as it is applied using an function that references amplitude, the idea is that the
+threshold is low for high amplitudes, e.g. 5 degrees for 1.0m and high for low amplitudes 80 degrees for 0.01 m.
+
+Amplitudes at phase exceedance locataions are also returned to allow assessment of the impact, e.g. low amplitude low impact
+
 """
 from netCDF4 import Dataset
 import numpy as np
@@ -47,7 +56,7 @@
 
 # TODO: add TPXO read and subset functionality currently only uses FES as "truth"
 
-def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.1,model='fes',model_res=1/16):
+def main(bdy_file='inputs/namelist_cmems.bdy',model='fes'):
     logger.info('============================================')
     logger.info('Start Tide Test Logging: ' + time.asctime())
     logger.info('============================================')
@@ -60,7 +69,7 @@ def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.1,model='f
     if model == 'fes':
         logger.info('using FES as reference.......')
         # open writer object to write pandas dataframes to spreadsheet
-        writer = pd.ExcelWriter(settings['dst_dir'] + 'exceed_values_amp_thres-'+str(amplitude_threshold)+'_reference_model-'+str(model)+'.xlsx', engine='xlsxwriter')
+        writer = pd.ExcelWriter(settings['dst_dir'] + 'comparision_with_'+str(model)+'.xlsx', engine='xlsxwriter')
         for key in constituents:
             for j in range(len(grids)):
                 out_fname = settings['dst_dir']+settings['fn']+'_bdytide_'+constituents[key].strip("',/\n")+'_grd_'+grids[j]+'.nc'
@@ -74,7 +83,7 @@ def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.1,model='f
                 # subset FES to match PyNEMO list of lat lons
                 subset_fes = subset_reference(pynemo_out, fes)
                 # compare the two lists (or dicts really)
-                error_log = compare_tides(pynemo_out, subset_fes, amplitude_threshold, model_res)
+                error_log = compare_tides(pynemo_out, subset_fes)
                 # return differences above threshold as a Pandas Dataframe and name using HC and Grid
                 error_log.name = constituents[key].strip("',/\n") + grids[j]
                 # if the dataframe is empty (no exceedances) then discard dataframe and log the good news
@@ -150,6 +159,7 @@ def read_fes(fes_fname,grid):
     # subset FES dict from read_FES, this uses find_nearest to find nearest FES point using PyNEMO dict from extract_PyNEMO
     # It also converts FES amplitude from cm to m.
 def subset_reference(pynemo_out, reference):
+    model_res = np.abs(reference['lon'][0]-reference['lon'][1])
     idx_lat = np.zeros(np.shape(pynemo_out['lat']))
     for i in range(np.shape(pynemo_out['lat'])[1]):
         idx_lat[0, i] = find_nearest(reference['lat'], pynemo_out['lat'][0, i])
@@ -167,7 +177,7 @@ def subset_reference(pynemo_out, reference):
         for i in range(np.shape(amp_sub)[1]):
             # if a fill value in FES subset is found
             if amp_sub[0, i] == 184467436613926912.0000:
-                logger.warning('found fill value in FES subset, taking nanmean from surrounding points')
+                logger.warning('found fill value in FES subset, taking nanmean from surrounding amplitude points')
                 # if there are fill values surrounding subset fill value change these to NaN
                 if reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]]== 184467436613926912.0000:
                     reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]] = np.nan
@@ -187,7 +197,7 @@ def subset_reference(pynemo_out, reference):
                     reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]-1] = np.nan
                 # nan mean surrounding points to replace fill value subset point
                 amp_sub[0,i] = np.nanmean([reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]], \
-                               reference['amp'][idx_lat[0,i], idx_lon[0,i]+1], \
+                              reference['amp'][idx_lat[0,i], idx_lon[0,i]+1], \
                               reference['amp'][idx_lat[0,i]-1, idx_lon[0,i]], \
                               reference['amp'][idx_lat[0,i], idx_lon[0,i]-1], \
                               reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]]+1, \
@@ -197,8 +207,10 @@ def subset_reference(pynemo_out, reference):
                                            ])
         phase_sub = reference['phase'][idx_lat, idx_lon]
         for i in range(np.shape(phase_sub)[1]):
+            # if a fill value in FES subset is found
             if phase_sub[0, i] == 18446744073709551616.0000:
-                logger.warning('found fill value in FES subset, taking nanmean value from surrounding points')
+                logger.warning('found fill value in FES subset, taking nanmean from surrounding phase points')
+                # if there are fill values surrounding subset fill value change these to NaN
                 if reference['phase'][idx_lat[0, i] + 1, idx_lon[0, i]] == 18446744073709551616.0000:
                     reference['phase'][idx_lat[0, i] + 1, idx_lon[0, i]] = np.nan
                 if reference['phase'][idx_lat[0, i], idx_lon[0, i] + 1] == 18446744073709551616.0000:
@@ -215,7 +227,7 @@ def subset_reference(pynemo_out, reference):
                     reference['phase'][idx_lat[0, i] - 1, idx_lon[0, i] + 1] = np.nan
                 if reference['phase'][idx_lat[0, i] + 1, idx_lon[0, i] - 1] == 18446744073709551616.0000:
                     reference['phase'][idx_lat[0, i] + 1, idx_lon[0, i] - 1] = np.nan
-
+                # calculate HcosG and then average
                 HcosG = np.nanmean([reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]]*np.cos(
                              reference['phase'][idx_lat[0, i]+1, idx_lon[0, i]]*np.pi/180),
                          reference['amp'][idx_lat[0, i], idx_lon[0, i]+1] * np.cos(
@@ -233,7 +245,7 @@ def subset_reference(pynemo_out, reference):
                          reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]-1] * np.cos(
                              reference['phase'][idx_lat[0, i]+1, idx_lon[0, i]-1] * np.pi / 180),
                          ])
-
+                # calculate HsinG and then average
                 HsinG = np.nanmean([reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]]*np.sin(
                              reference['phase'][idx_lat[0, i]+1, idx_lon[0, i]]*np.pi/180),
                          reference['amp'][idx_lat[0, i], idx_lon[0, i]+1] * np.sin(
@@ -251,33 +263,39 @@ def subset_reference(pynemo_out, reference):
                          reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]-1] * np.sin(
                              reference['phase'][idx_lat[0, i]+1, idx_lon[0, i]-1] * np.pi / 180),
                          ])
-
+                # convert back to phase
                 phase_sub[0,i] = np.arctan2(HsinG,HcosG)
 
     lat_sub = reference['lat'][idx_lat]
     lon_sub = reference['lon'][idx_lon]
-    subset = {'lat':lat_sub,'lon':lon_sub,'amp':amp_sub,'phase':phase_sub}
+    subset = {'lat':lat_sub,'lon':lon_sub,'amp':amp_sub,'phase':phase_sub,'model_res':model_res}
     return subset
 
     # takes pynemo extract dict, subset fes dict, and the thresholds and model res passed to main function.
     # returns a Pandas Dataframe with any PyNEMO values that exceed the nearest FES point by defined threshold
     # It also checks lats and lons are within the model reference resolution
     # i.e. ensure closest model reference point is used.
-def compare_tides(pynemo_out,subset,amp_thres,model_res):
+def compare_tides(pynemo_out,subset):
     # compare lat and lons
     diff_lat = np.abs(pynemo_out['lat']-subset['lat'])
     diff_lon = np.abs(pynemo_out['lon'] - subset['lon'])
-    exceed_lat = diff_lat > model_res
-    exceed_lon = diff_lon > model_res
+    exceed_lat = diff_lat > subset['model_res']
+    exceed_lon = diff_lon > subset['model_res']
     exceed_sum = np.sum(exceed_lat+exceed_lon)
     if exceed_sum > 0:
         raise Exception('Dont Panic: Lat and/or Lon further away from model point than model resolution')
     # surpress warnings as NaNs from averaging surrounding pixels can cause issues
     with warnings.catch_warnings():
         warnings.simplefilter("ignore", category=RuntimeWarning)
         # compare amp
-        abs_amp = np.abs(pynemo_out['amp']-subset['amp'])
-        abs_amp_thres = abs_amp > amp_thres
+        abs_amp_diff = np.abs(pynemo_out['amp']-subset['amp'])
+        # calculate threshold in percentage terms
+        logger.info('percentage amplitude exceedance calculated using the following.....')
+        amp_percentage_exceed = 26.933 * subset['amp'] ** -0.396
+        logger.info('Percentage Exceedance = 26.933 * Reference Amplitude ^ -0.396')
+        # work out difference based on percentage and reference amplitude
+        percent_diff = (abs_amp_diff / pynemo_out['amp']) * 100
+        abs_amp_thres = percent_diff > amp_percentage_exceed
         err_amp = pynemo_out['amp'][abs_amp_thres].tolist()
         err_amp_lats = pynemo_out['lat'][abs_amp_thres].tolist()
         err_amp_lons = pynemo_out['lon'][abs_amp_thres].tolist()
@@ -297,7 +315,7 @@ def compare_tides(pynemo_out,subset,amp_thres,model_res):
         abs_ph[abs_ph < 0.0 ] = abs_ph[abs_ph < 0.0] *-1
         # calculate phase threshold based on amplitude and power relationship
         # as amplitude decreases the phase exceedance allowed increases.
-        logger.info('phase exccedance calculates using the following.....')
+        logger.info('phase exceedance calculated using the following.....')
         phase_thres = 5.052 * pynemo_out['amp'] ** -0.60
         logger.info('Exceedance = 5.052 * Amplitude ^ -0.60')
         abs_ph_thres = abs_ph > phase_thres
diff --git a/pynemo/tide/fes_extract_HC.py b/pynemo/tide/fes_extract_HC.py
@@ -48,8 +48,8 @@ def __init__(self, settings, lat, lon, grid_type):
             self.mask_dataset = {}
 
             # extract lon and lat z data
-            lon_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lon'])
-            lat_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lat'])
+            lon_z = np.array(Dataset(settings['tide_fes']+constituents[i]+'_Z.nc').variables['lon'])
+            lat_z = np.array(Dataset(settings['tide_fes']+constituents[i]+'_Z.nc').variables['lat'])
             lon_resolution = lon_z[1] - lon_z[0]
             data_in_km = 0 # added to maintain the reference to matlab tmd code
 
@@ -63,8 +63,8 @@ def __init__(self, settings, lat, lon, grid_type):
                 #read and convert the height_dataset file to complex and store in dicts
                 hRe = []
                 hIm = []
-                lat_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lat'][:])
-                lon_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lon'][:])
+                lat_z = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_Z.nc').variables['lat'][:])
+                lon_z = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_Z.nc').variables['lon'][:])
                 for ncon in range(len(constituents)):
                     amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:])))
                     # set fill values to zero
@@ -90,14 +90,12 @@ def __init__(self, settings, lat, lon, grid_type):
 
                 URe = []
                 UIm = []
-                lat_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lat'][:])
-                lon_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lon'][:])
+                lat_u = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_U.nc').variables['lat'][:])
+                lon_u = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_U.nc').variables['lon'][:])
                 for ncon in range(len(constituents)):
                     amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
                     # set fill values to zero
                     amp[amp == 18446744073709551616.00000] = 0
-                    # convert amp units to m/s
-                    amp = amp/100.00
                     phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
                     phase[phase == 18446744073709551616.00000] = 0
                     URe.append(amp*np.cos(phase*(np.pi/180)))
@@ -114,14 +112,12 @@ def __init__(self, settings, lat, lon, grid_type):
 
                 VRe = []
                 VIm = []
-                lat_v = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['lat'][:])
-                lon_v = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['lon'][:])
+                lat_v = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_V.nc').variables['lat'][:])
+                lon_v = np.array(Dataset(settings['tide_fes'] + constituents[i] + '_V.nc').variables['lon'][:])
                 for ncon in range(len(constituents)):
                     amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
                     # set fill value to zero
                     amp[amp == 18446744073709551616.00000] = 0
-                    # convert amp units to m/s
-                    amp = amp/100.00
                     phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
                     phase[phase == 18446744073709551616.00000] = 0
                     VRe.append(amp*np.cos(phase*(np.pi/180)))
