added tide checker function to BDY file

inputs/namelist_cmems.bdy

@@ -110,8 +110,11 @@
     !clname(3)      = 'O1'
     !clname(4)      = 'K1'
     !clname(5)      = 'N2'
-    ln_trans       = .false.               !  interpolate transport rather than
-                                          !  velocities
+    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.25                  ! amplitude thresold to compare against (m)
+    nn_phase_thres   = 20.0                 ! phase threshold to compare against (degrees)
 !------------------------------------------------------------------------------
 !  Time information
 !------------------------------------------------------------------------------

pynemo/profile.py

@@ -54,6 +54,7 @@ from pynemo.reader.factory import GetFile
 from pynemo.reader import factory
 from pynemo.tide import nemo_bdy_tide3 as tide
 from pynemo.tide import nemo_bdy_tide_ncgen
+from pynemo.tests import nemo_tide_test as tt
 from pynemo.utils import Constants
 from pynemo.gui.nemo_bdy_mask import Mask as Mask_File
 from pynemo import nemo_bdy_dl_cmems as dl_cmems
@@ -437,13 +438,26 @@ def process_bdy(setup_filepath=0, mask_gui=False):
         if settings['tide_model']=='tpxo':
             cons = tide.nemo_bdy_tide_rot(
                 Setup, DstCoord, bdy_ind['t'], bdy_ind['u'], bdy_ind['v'],
-                                                            settings['clname'], settings['tide_model'])
+                                                               settings['clname'], settings['tide_model'])
+            if settings['tide_checker'] == True:
+                logger.info('tide checker starting now.....')
+                tt_test = tt.main(setup_filepath,settings['amp_thres'],settings['phase_thres'],settings['ref_model'])
+                if tt_test == 0:
+                    logger.info('tide checker ran successfully, check spreadsheet in output folder')
+                if tt_test !=0:
+                    logger.warning('error running tide checker')
+
         elif settings['tide_model']=='fes':
             cons = tide.nemo_bdy_tide_rot(
                 Setup, DstCoord, bdy_ind['t'], bdy_ind['u'], bdy_ind['v'],
                                                             settings['clname'],settings['tide_model'])
-            logger.warning('Tidal model: %s, not yet properly implimented',
-                         settings['tide_model'])
+            if settings['tide_checker'] == True:
+                logger.info('tide checker starting now.....')
+                tt_test = tt.main(setup_filepath,settings['amp_thres'],settings['phase_thres'],settings['ref_model'])
+                if tt_test == 0:
+                    logger.info('tide checker ran successfully, check spreadsheet in output folder')
+                if tt_test !=0:
+                    logger.warning('error running tide checker')
         else:
             logger.error('Tidal model: %s, not recognised', 
                          settings['tide_model'])

pynemo/tests/nemo_tide_test.py

@@ -46,7 +46,7 @@ logging.basicConfig(filename='nrct.log', level=logging.INFO)
 
 # TODO: add TPXO read and subset functionality currently only uses FES as "truth"
 
-def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_threshold=15.00,model_res=1/16,model='fes'):
+def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_threshold=15.00,model='fes',model_res=1/16):
     logger.info('============================================')
     logger.info('Start Tide Test Logging: ' + time.asctime())
     logger.info('============================================')
@@ -85,7 +85,7 @@ def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_t
                     logger.warning('Exceedance in thesholds detected, check spreadsheet in dst_dir')
                     error_log.to_excel(writer,sheet_name=error_log.name)
                 # close writer object and save excel spreadsheet
-                writer.save()
+        writer.save()
     # code runs here if TPXO is requested as reference this hasn't been written yet so raises exception
     elif model == 'tpxo':
         logger.info('using TPXO as reference.......')
@@ -129,11 +129,10 @@ def read_fes(fes_fname,grid):
         fes_amp = np.array(fes_tide.variables['amplitude'][:])
         fes_amp = fes_amp / 100
         fes_phase = np.array(fes_tide.variables['phase'][:])
-        fes_phase[fes_phase > 180.0] = fes_phase[fes_phase > 180.0] - 360.0
+
     if grid != 'Z':
         fes_amp = np.array(fes_tide.variables[grid+'a'][:])
         fes_phase = np.array(fes_tide.variables[grid+'g'][:])
-        fes_phase[fes_phase > 180.0] = fes_phase[fes_phase > 180.0] - 360.0
 
     fes_lat = fes_tide.variables['lat'][:]
     fes_lon = fes_tide.variables['lon'][:]
@@ -175,6 +174,7 @@ def compare_tides(pynemo_out,subset,amp_thres,phase_thres,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')
+
     # compare amp
     abs_amp = np.abs(pynemo_out['amp']-subset['amp'])
     abs_amp_thres = abs_amp > amp_thres
@@ -187,8 +187,16 @@ def compare_tides(pynemo_out,subset,amp_thres,phase_thres,model_res):
     err_ref_lons_amp = subset['lon'][abs_amp_thres].tolist()
 
     # compare phase
-    abs_ph = np.abs(pynemo_out['phase']-subset['phase'])
+    # change from -180-180 to 0to 360 for both pynemo and subset.
+    pynemo_out['phase'][pynemo_out['phase'] < 0.0] = pynemo_out['phase'][pynemo_out['phase'] < 0.0] + 360.0
+    subset['phase'][subset['phase'] < 0.0] = subset['phase'][subset['phase'] < 0.0] + 360.0
+    # compare phase angles between 0 and 360.
+    abs_ph = 180 - abs(abs(pynemo_out['phase'] - subset['phase']) - 180)
+    # values outside of 0 to 360  (such as erroneous fill values) end up negative
+    # so multiply by -1 to ensure they are identified as exceeding threshold
+    abs_ph[abs_ph < 0.0 ] = abs_ph[abs_ph < 0.0] *-1
     abs_ph_thres = abs_ph > phase_thres
+    
     err_pha = pynemo_out['phase'][abs_ph_thres[0,:]].tolist()
     err_pha_lats = pynemo_out['lat'][abs_ph_thres].tolist()
     err_pha_lons = pynemo_out['lon'][abs_ph_thres].tolist()

pynemo/tide/fes_extract_HC.py

@@ -53,77 +53,82 @@ class HcExtract(object):
             lon_resolution = lon_z[1] - lon_z[0]
             data_in_km = 0 # added to maintain the reference to matlab tmd code
 
-            # extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
-            mask_z = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:])))
-            mask_z[mask_z != 18446744073709551616.00000] = 1
-            mask_z[mask_z == 18446744073709551616.00000] = 0
-            self.mask_dataset[mz_name] = mask_z
-
-            mask_u = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:])))
-            mask_z[mask_z != 18446744073709551616.00000] = 1
-            mask_u[mask_u == 18446744073709551616.00000] = 0
-            self.mask_dataset[mu_name] = mask_u
-
-            mask_v = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:])))
-            mask_z[mask_z != 18446744073709551616.00000] = 1
-            mask_v[mask_v == 18446744073709551616.00000] = 0
-            self.mask_dataset[mv_name] = mask_v
-
-            #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'][:])
-            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
-                amp[amp == 18446744073709551616.00000] = 0
-                # convert to m from cm
-                amp = amp/100.00
-                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
-                phase[phase == 18446744073709551616.00000] = 0
-                hRe.append(amp*np.cos(phase*(np.pi/180)))
-                hIm.append(-amp*np.sin(phase*(np.pi/180)))
-            hRe = np.stack(hRe)
-            hIm = np.stack(hIm)
-            self.height_dataset = [lon_z,lat_z,hRe,hIm]
-
-            #read and convert the velocity_dataset files to complex
-            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'][:])
-            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 to cm from m
-                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)))
-                UIm.append(-amp*np.sin(phase*(np.pi/180)))
-            URe = np.stack(URe)
-            UIm = np.stack(UIm)
-            self.Uvelocity_dataset = [lon_u,lat_u,URe,UIm]
-
-            VRe = []
-            VIm = []
-            lat_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lat'][:])
-            lon_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_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 m to cm
-                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)))
-                VIm.append(-amp*np.sin(phase*(np.pi/180)))
-            VRe = np.stack(VRe)
-            VIm = np.stack(VIm)
-            self.Vvelocity_dataset = [lon_v,lat_v,VRe,VIm]
+            if grid_type == 'z' or grid_type == 't':
+                # extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
+                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['amplitude'][:])))
+                mask[mask != 18446744073709551616.00000] = 1
+                mask[mask == 18446744073709551616.00000] = 0
+                self.mask_dataset[mz_name] = mask
+
+                #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'][:])
+                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
+                    amp[amp == 18446744073709551616.00000] = 0
+                    # convert to m from cm
+                    amp = amp/100.00
+                    phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
+                    # set fill values to 0
+                    phase[phase == 18446744073709551616.00000] = 0
+                    # convert to real and imaginary conjugates and also convert to radians.
+                    hRe.append(amp*np.cos(phase*(np.pi/180)))
+                    hIm.append(-amp*np.sin(phase*(np.pi/180)))
+                hRe = np.stack(hRe)
+                hIm = np.stack(hIm)
+                self.height_dataset = [lon_z,lat_z,hRe,hIm]
+
+            elif grid_type == 'u':
+                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['Ua'][:])))
+                mask[mask != 18446744073709551616.00000] = 1
+                mask[mask == 18446744073709551616.00000] = 0
+                self.mask_dataset[mu_name] = mask
+                #read and convert the velocity_dataset files to complex
+
+                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'][:])
+                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 to cm from m
+                    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)))
+                    UIm.append(-amp*np.sin(phase*(np.pi/180)))
+                URe = np.stack(URe)
+                UIm = np.stack(UIm)
+                self.Uvelocity_dataset = [lon_u,lat_u,URe,UIm]
+
+            elif grid_type == 'v':
+                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['Va'][:])))
+                mask[mask != 18446744073709551616.00000] = 1
+                mask[mask == 18446744073709551616.00000] = 0
+                self.mask_dataset[mv_name] = mask
+
+                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'][:])
+                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 m to cm
+                    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)))
+                    VIm.append(-amp*np.sin(phase*(np.pi/180)))
+                VRe = np.stack(VRe)
+                VIm = np.stack(VIm)
+                self.Vvelocity_dataset = [lon_v,lat_v,VRe,VIm]
 
             # open grid variables these are resampled TPXO parameters so may not work correctly.
             self.grid = Dataset(settings['tide_fes']+'grid_fes.nc')
@@ -139,7 +144,7 @@ class HcExtract(object):
         if glob == 1:
             lon_z = np.concatenate(([lon_z[0]-lon_resolution, ], lon_z,[lon_z[-1]+lon_resolution, ]))
             height_z = np.concatenate(([height_z[-1, :], ], height_z, [height_z[0, :],]), axis=0)
-            mask_z = np.concatenate(([mask_z[-1, :], ], mask_z, [mask_z[0, :], ]), axis=0)
+            mask_z = np.concatenate(([mask[-1, :], ], mask, [mask[0, :], ]), axis=0)
 
         #adjust lon convention
         xmin = np.min(lon)
@@ -250,9 +255,9 @@ class HcExtract(object):
                                                                 maskedpoints, lonlat)
 
             #for velocity_dataset values
-            if height_data is not None:
-                data_temp[cons_index, :, 0] = data_temp[cons_index, :, 0]/height_data*100
-                data_temp[cons_index, :, 1] = data_temp[cons_index, :, 1]/height_data*100
+            #if height_data is not None:
+            #    data_temp[cons_index, :, 0] = data_temp[cons_index, :, 0]/height_data*100
+            #    data_temp[cons_index, :, 1] = data_temp[cons_index, :, 1]/height_data*100
 
             zcomplex = np.array(data_temp[cons_index, :, 0], dtype=complex)
             zcomplex.imag = data_temp[cons_index, :, 1]