From 0e3ad301365c5a4bb5d60f3ccc01d8e8f9575cfa Mon Sep 17 00:00:00 2001
From: thopri <thopri@noc.ac.uk>
Date: Thu, 7 May 2020 15:30:19 +0100
Subject: [PATCH] added tide checker function to BDY file
---
inputs/namelist_cmems.bdy | 7 +-
pynemo/profile.py | 20 ++++-
pynemo/tests/nemo_tide_test.py | 18 ++--
pynemo/tide/fes_extract_HC.py | 155 +++++++++++++++++----------------
4 files changed, 115 insertions(+), 85 deletions(-)
diff --git a/inputs/namelist_cmems.bdy b/inputs/namelist_cmems.bdy
index af97def..353da6d 100644
--- a/inputs/namelist_cmems.bdy
+++ b/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
!------------------------------------------------------------------------------
diff --git a/pynemo/profile.py b/pynemo/profile.py
index 0f126a8..2d1e8f8 100644
--- a/pynemo/profile.py
+++ b/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'])
diff --git a/pynemo/tests/nemo_tide_test.py b/pynemo/tests/nemo_tide_test.py
index 0a8afa7..0f39575 100644
--- a/pynemo/tests/nemo_tide_test.py
+++ b/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()
diff --git a/pynemo/tide/fes_extract_HC.py b/pynemo/tide/fes_extract_HC.py
index 04f50de..8822f54 100644
--- a/pynemo/tide/fes_extract_HC.py
+++ b/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]
--
GitLab