Using extract_netcdf_constants to read TPXO9-atlas-v5 files
Closed this issue · 12 comments
Hello,
I’m trying to use extract_netcdf_constants to calculate constituents from TPXO9-atlas-v5 files. It is working great when nearest neighbours interpolation is used, but not when trying bilinear interpolation. It still gives the correct phase, but the amplitudes are significantly off.
Example code:
from pyTMD.read_netcdf_model import extract_netcdf_constants
import glob
# Coordinates
lat = 50.103
lon = -5.54283
# Listing all constituent files
files = glob.glob(tpx_path + "/h/*.nc")
# Assigning grid file
grid = glob.glob(tpx_path + "/grid/*.nc")[0]
amp, ph, D, c = extract_netcdf_constants(lon, lat, grid, files, TYPE="z", GZIP=False, METHOD="bilinear",
EXTRAPOLATE=True, SCALE=1.0/1000)
Sorry to bother, I’m new to this field so it is likely something I am doing wrong. I must say I have found pyTMD extremely useful so far. thank you!
Hi @davidegregg,
Great question. I looked into this. I was finding that as long as extrapolate is set (since this coordinate is masked in FES), that the method of interpolation was largely negligible. Let me know if you were finding things different.
Hi @tsutterley, thanks for the response.
All works well for FES, but it's when I try to use the function for OSU TPXO9-atlas-v5 files that I start having problems (still with extrapolate set).
Right, sorry @davidegregg my brain has been in a bit of a fog over the past week.
There was a problem in the bilinear interpolator for interpolating the mask within TPXO9. I think everything should be fixed and ready to go.
Brilliant, thank you!
Hi again @tsutterley,
Sorry to open again, but even after this fix, I am getting much higher residuals for TPXO9-atlas-v5 when using "bilinear" compared to "nearest" which seems unusual.
I'll look into it today. Sorry for all the trouble
No problem at all, this package is extremely helpful for me! For reference, below I've put phasor maps I've made of the residuals in predictions of M2 compared to tide gauge observations to show the difference between the two interpolation methods. The first uses nearest neighbours, the second uses bilinear (a couple of other model predictions using bilinear are there too for comparison).
Hey @davidegregg,
I'm having a hard time recreating this difference. I took the location of the tide gauge that seemed to be one of the most prominent (Newport at 53.63018°N, -0.18742°E) and ran the tide prediction notebook.
Here's bilinear with TPXO9-v5-alpha:
And here's spline with TPXO9-v5-alpha:
And here's the tide prediction from the National Tide and Sea Level Facility:
I'm not really finding a difference between the interpolation methods here, and the model seems to be not too bad compared with the gauge. Can you help me figure out what's happening on your end?
Hey @tsutterley,
Sorry I might have explained the problem slightly wrong. For those coordinates (those are actually for the Immingham site I believe), the model performs reasonably similar for bilinear interpolation and nearest neighbours, this is just a site where TPXO9-atlas performs somewhat poorly.
But, if I take the North Shields tide gauge as an example (55.00744, -1.43978), I'm getting drastically different results between the two methods. Here's the code I'm using:
from pyTMD.read_netcdf_model import extract_netcdf_constants
import glob
# Coordinates for North Shields tide gauge
lat = 55.00744
lon = -1.43978
# Listing all constituent files
files = glob.glob(tpx_path + "/h/*.nc")
# Assigning grid file
grid = glob.glob(tpx_path + "/grid/*.nc")[0]
# Bilinear method
amp, ph, D, c = extract_netcdf_constants(lon, lat, grid, files, TYPE="z", GZIP=False, METHOD="bilinear",
EXTRAPOLATE=True, SCALE=1.0 / 1000)
# Nearest neighbours method
amp2, ph2, D2, c2 = extract_netcdf_constants(lon, lat, grid, files, TYPE="z", GZIP=False, METHOD="nearest",
EXTRAPOLATE=True, SCALE=1.0 / 1000)
I'm getting the following amplitudes for M2:
nearest: 1.603 m (near perfect prediction for the observed tide gauge data)
bilinear: 0.288 m (significantly off)
I hope that clears things up and thanks for all the help so far!
Thanks! This helped tremendously. I think I tracked down the root of the problem and will have it fixed in PR #105. Here's the updated results for North Shields now that I fixed the masking problem.
Bilinear:
Spline:
(and you were right that was the lat/lon of the Immingham site. I ran Immingham after Newport but the interpolation results were similar so I just included the plots for one site).
That's all working now, thanks!
Hi @tsutterley and @davidegregg,
I learned from you here to extract amplitudes for M2 at spacific location.
I am wondering is there any example I can learn to make a M2 map for a loacl area using pyTMD, like the following figure?
Looking forward to your help. Thanks.
Jiangchao