Extreme large u,v current near coast

[cywhale]

I use TPXO9-atlas-v5 model, and the u,v current get from transport divided by depth (U/D, V/D) in pyTMD code seems get extreme large u, v near coastal region (shallow water). For example: (Please correct me if I had any mistakes in code)

from pyTMD.io import model, ATLAS
from pyTMD import predict
from pyTMD.time import convert_calendar_dates, datetime_to_list
from datetime import datetime, timedelta

def get_current_model(model_name, model_directory, model_format, compressed=False):
    return model(model_directory, format=model_format, compressed=compressed).current(model_name)

tpxo_model_directory = './'
tpxo_model_format = 'netcdf'
tpxo_compressed = False
tpxo_model_name = 'TPXO9-atlas-v5'

tpxo_curr_model = get_current_model(tpxo_model_name, tpxo_model_directory, tpxo_model_format, tpxo_compressed)

x = [118.033332   118.19999866 118.36666533 118.53333199 118.69999866
 118.86666532 119.03333199 119.19999865]
y= [22.83333424 22.83333424 22.83333424 22.83333424 22.83333424 22.83333424
 22.83333424 22.83333424]
 
ampu, phu, D, c = ATLAS.extract_constants(
        x, y, tpxo_curr_model.grid_file,
        tpxo_curr_model.model_file['u'], type='u', method='spline',
        scale=tpxo_curr_model.scale, compressed=tpxo_curr_model.compressed
    )
print(ampu)
# [[29.881492163964104 58.07633346386415 259.1972825702797
#  368.21261192563594 633.5693187790492 367.0614087475642
#  9.981625336558242 334.24654747274866 377.5816565704873
# 461.38743839136697 157.4428182506624 24.65747654388396
#  107.1217650051865 0.05104605427265156 0.02552302713632578]....

ampv, phv, D, c = ATLAS.extract_constants(
        x, y, tpxo_curr_model.grid_file,
        tpxo_curr_model.model_file['v'], type='v', method='spline',
        scale=tpxo_curr_model.scale, compressed=tpxo_curr_model.compressed
    )

start_date = datetime(2023, 7, 25)
tide_time = []
dtlist = datetime_to_list(start_date)
tide_time.append(convert_calendar_dates(dtlist[0], dtlist[1], dtlist[2], dtlist[3]))
DELTAT = np.zeros_like(tide_time[0:1])

ny = 8
nx = 8
timelen = len(tide_time[0:1])
tide = {}

for TYPE in ['u', 'v']:
        if TYPE == 'u':
            amp = ampu
            ph = phu
        else:
            amp = ampv
            ph = phv

        shpx = amp.shape
        ampx = amp.reshape((shpx[0] * 1, shpx[1]))
        phx = ph.reshape((shpx[0] * 1, shpx[1]))
        # calculate complex phase in radians for Euler's
        cph = -1j * phx * np.pi / 180.0
        # calculate constituent oscillation
        hc = ampx * np.exp(cph)
        # Create a mask where values are NA or 0
        mask = np.isnan(hc) | (hc == 0)
        # Convert hc to a masked array
        hc = ma.array(hc, mask=mask) 

        TIDE = predict.map(tide_time[0], hc, c,
                           deltat=DELTAT[0], corrections=format)
        MINOR = predict.infer_minor(
            tide_time[0], hc, c, deltat=DELTAT[0], corrections=format)
        tx = TIDE+MINOR
        tx.data[tx.mask] = np.nan
        tide[TYPE] = tx.data

print(tide)
#{'u': array([  28.08796695,  289.27854635, 1033.66448834, 1633.89806405, 676.43379783,  859.09533155, 1763.9204305 , 1323.54183509]), 
#'v': array([-5833.34293851, -7885.86603018, -8757.43001407, -8943.87019027, -9169.67130404, -7707.99313235, -6800.35841252, -5431.25757682])}

print(np.sqrt((tide["u"] * 0.01)**2 + (tide["v"] * 0.01)**2))
[58.33411291 78.91169486 88.1822182  90.91889072 93.2165703  77.55720609 70.25402503 55.90198183]

As the result, it get extreme large u, v and thus magnitude > 50 m/s. Seems not reasonable values for tidal current, and that's not special cases, which occurred in most of the coastal regions. I'm confused about these results, and think if it's because I misused the TPXO9 models in code, and applied it in wrong use cases. Any comments about my code and results, or guides to the usage would be appreciated. Thanks a lot.

[tsutterley]

Hi @cywhale, try with the extrapolate option when extracting constants. That might get you more realistic values.

[cywhale]

Hi, thanks. I had tried

amp, ph, D, c = ATLAS.extract_constants(
        x, y, tpxo_curr_model.grid_file,
        tpxo_curr_model.model_file['u'], type='u', method='spline',
        scale=tpxo_curr_model.scale, compressed=tpxo_curr_model.compressed,
        extrapolate=True
    )

Do you mean this option?
But it produces almost the same results. I read the source code, the exptrapolate method in extract_constants() just apply to the points of invalid data points with mask is True. So I think it cannot affect the results when these data points are originally valid?
https://github.com/tsutterley/pyTMD/blob/main/pyTMD/io/ATLAS.py#L323

[tsutterley]

Was worth a shot. I hope that all the masked cases work, but I'm not 100% confident.
u and v are velocities in cm/s. We typically try to match the standards from the TMD tide program. I should update the documentation to note the typical units for outputs.

[cywhale]

Thanks for the reply. I check the mask is mostly on land, but it seems that some mased cases are for invalid data points area, right? I plot the masked area in this example:

I also check a u,v map, and the extreme large tidal current are shown. I want to make sure that the values are already multiplied by 0.01 to convert from cm/s to m/s: np.sqrt((tide["u"] * 0.01)**2 + (tide["v"] * 0.01)**2)) (so the unit of legend on this plot is m/s, around 3/4 values are under normally 5m/s, but other 1/4 values are unusual large, even > 50 m/s). I had tried to calculate/plot many other regions worldwide with TPXO9-atlas-v5 model by pyTMD, and it seems that extreme u,v value around coastal regions are general encountered.

[tsutterley]

You're right. I think the built-in scale value is wrong for TPXO9-atlas currents. I think it should be 1e-4 instead of 1e-2. I'll put a fix in today.

[tsutterley]

in the mean time, you can just set tpxo_curr_model.scale to 1e-4 before extracting the constants

[tsutterley]

this should be fixed in main after PR #243. Let me know if you need anything.

[cywhale]

Thanks for the quick fix. It's a great help.