Skip to content

GitLab

Commit d5d9f5c9 authored by thopri's avatar thopri
Browse files

updated testing function to cope wiht fill values in FES

parent 602269e8
Show whitespace changes
Inline Side-by-side
Showing with 172 additions and 73 deletions
inputs/namelist_cmems.bdy
View file @ d5d9f5c9
......@@ -113,8 +113,7 @@
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)
nn_amp_thres = 0.30 ! amplitude thresold to compare against (m)
!------------------------------------------------------------------------------
! Time information
!------------------------------------------------------------------------------
......
pynemo/profile.py
View file @ d5d9f5c9
......@@ -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['phase_thres'],settings['ref_model'])
tt_test = tt.main(setup_filepath,settings['amp_thres'],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['phase_thres'],settings['ref_model'])
tt_test = tt.main(setup_filepath,settings['amp_thres'],settings['ref_model'])
if tt_test == 0:
logger.info('tide checker ran successfully, check spreadsheet in output folder')
if tt_test !=0:
......
pynemo/tests/nemo_tide_test.py
View file @ d5d9f5c9
......@@ -38,6 +38,7 @@ import numpy as np
import logging
import time
import pandas as pd
import warnings
from pynemo import nemo_bdy_setup as setup
# log to PyNEMO log file
......@@ -46,7 +47,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='fes',model_res=1/16):
def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,model='fes',model_res=1/16):
logger.info('============================================')
logger.info('Start Tide Test Logging: ' + time.asctime())
logger.info('============================================')
......@@ -59,7 +60,7 @@ def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_t
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)+'_phase_thres-'+str(phase_threshold)+'_reference_model-'+str(model)+'.xlsx', engine='xlsxwriter')
writer = pd.ExcelWriter(settings['dst_dir'] + 'exceed_values_amp_thres-'+str(amplitude_threshold)+'_reference_model-'+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'
......@@ -73,7 +74,7 @@ def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_t
# 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, phase_threshold, model_res)
error_log = compare_tides(pynemo_out, subset_fes, amplitude_threshold, model_res)
# 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
......@@ -119,6 +120,7 @@ def extract_PyNEMO_output(out_fname,grid):
lat = np.array(nav_lat[nbjdta, nbidta])
lon = np.array(nav_lon[nbjdta, nbidta])
pynemo_out = {'lat':lat,'lon':lon,'amp':amp,'phase':phase}
tide_out.close()
return pynemo_out
# read FES netcdf file, convert lon to -180 to 180, rather than 0-360 it also converts amplitude from CM to M
......@@ -142,6 +144,7 @@ def read_fes(fes_fname,grid):
# change to -180 to 180 lonitude convention
fes_lon[fes_lon > 180.0] = fes_lon[fes_lon > 180.0] - 360.0
fes_dict = {'lat':fes_lat,'lon':fes_lon,'amp':fes_amp,'phase':fes_phase}
fes_tide.close()
return fes_dict
# subset FES dict from read_FES, this uses find_nearest to find nearest FES point using PyNEMO dict from extract_PyNEMO
......@@ -158,7 +161,98 @@ def subset_reference(pynemo_out, reference):
idx_lon = idx_lon.astype(np.int64)
amp_sub = reference['amp'][idx_lat, idx_lon]
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
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')
# 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
if reference['amp'][idx_lat[0,i], idx_lon[0,i]+1]== 184467436613926912.0000:
reference['amp'][idx_lat[0, i], idx_lon[0,i]+1] = np.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
if reference['amp'][idx_lat[0,i], idx_lon[0,i]-1]== 184467436613926912.0000:
reference['amp'][idx_lat[0, i], idx_lon[0, i]-1] = np.nan
if reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]+1]== 184467436613926912.0000:
reference['amp'][idx_lat[0, i]+1, idx_lon[0, i]+1] = np.nan
if reference['amp'][idx_lat[0,i]-1, idx_lon[0,i]-1]== 184467436613926912.0000:
reference['amp'][idx_lat[0, i]-1, idx_lon[0, i]-1] = np.nan
if reference['amp'][idx_lat[0,i]-1, idx_lon[0,i]+1]== 184467436613926912.0000:
reference['amp'][idx_lat[0, i]-1, idx_lon[0, i]+1] = np.nan
if reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]-1]== 184467436613926912.0000:
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]-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, \
reference['amp'][idx_lat[0,i]-1, idx_lon[0,i]-1], \
reference['amp'][idx_lat[0,i]-1, idx_lon[0,i]+1], \
reference['amp'][idx_lat[0,i]+1, idx_lon[0,i]-1] \
])
phase_sub = reference['phase'][idx_lat, idx_lon]
for i in range(np.shape(phase_sub)[1]):
if phase_sub[0, i] == 18446744073709551616.0000:
logger.warning('found fill value in FES subset, taking nanmean value from surrounding points')
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:
reference['phase'][idx_lat[0, i], idx_lon[0, i] + 1] = np.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:
reference['phase'][idx_lat[0, i], 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
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
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
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
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(
reference['phase'][idx_lat[0, i], idx_lon[0, i]+1] * np.pi / 180),
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(
reference['phase'][idx_lat[0, i], idx_lon[0, i]-1] * np.pi / 180),
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),
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),
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),
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),
])
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(
reference['phase'][idx_lat[0, i], idx_lon[0, i]+1] * np.pi / 180),
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(
reference['phase'][idx_lat[0, i], idx_lon[0, i]-1] * np.pi / 180),
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),
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),
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),
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),
])
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}
......@@ -168,7 +262,7 @@ def subset_reference(pynemo_out, reference):
# 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,phase_thres,model_res):
def compare_tides(pynemo_out,subset,amp_thres,model_res):
# compare lat and lons
diff_lat = np.abs(pynemo_out['lat']-subset['lat'])
diff_lon = np.abs(pynemo_out['lon'] - subset['lon'])
......@@ -177,7 +271,8 @@ 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')
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
......@@ -198,6 +293,11 @@ def compare_tides(pynemo_out,subset,amp_thres,phase_thres,model_res):
# 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
# 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.....')
phase_thres = 5.052 * pynemo_out['amp'] ** -0.60
logger.info('Exceedance = 5.052 * Amplitude ^ -0.60')
abs_ph_thres = abs_ph > phase_thres
err_pha = pynemo_out['phase'][abs_ph_thres[0,:]].tolist()
......
pynemo/tide/fes_extract_HC.py
View file @ d5d9f5c9
......@@ -133,6 +133,8 @@ class HcExtract(object):
# open grid variables these are resampled TPXO parameters so may not work correctly.
self.grid = Dataset(settings['tide_fes']+'grid_fes.nc')
height_z = np.array(np.rot90(self.grid.variables['hz']))
# set fill values to zero
height_z[height_z <0.00] = 0.00
else:
print('Don''t know that tide model')
......@@ -319,12 +321,10 @@ def bilinear_interpolation(lon, lat, data, lon_new, lat_new):
data_copy = data.copy()
data_copy[1:-1, 1:-1] = np.divide(height_temp, mask_temp)
nonmask_index = np.where(mask == 1)
data_copy[nonmask_index] = data[nonmask_index]
data_copy[data_copy == 0] = np.NaN
lonlat = np.concatenate((lon_new_copy, lat_new))
lonlat = np.reshape(lonlat, (lon_new_copy.size, 2), order='F')
result = interpolate.interpn((lon, lat), data_copy, lonlat)
return result
......
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