hydrogeoscience / pygtide

A Python module and wrapper for ETERNA PREDICT to compute gravitational tides on Earth
Mozilla Public License 2.0
50 stars 19 forks source link

Displacement artifact? #35

Open TraziBatine opened 1 year ago

TraziBatine commented 1 year ago

Hi. I was wondering if this shown here on the plot is an artefact or am I doing something wrong... I calculated displacement with:

pt.predict(lat, lon, startdate, duration, samprate, tidalpoten, tidalcompo=3, statazimut=0)

As you can see, there is a horizontal line at around -35, +35. test

Thanks for the info!

trichter commented 1 year ago

Hi. Interesting! Which observable did you plot? Can you share the code to generate this plot? A guess: Maybe +-35° is the latitude for which diurnal and semidiurnal constituents have equal strength? For the region including the equator semidiurnal tides dominate, for the regions including the poles diurnal tides dominate. There has to be a region in which the transition might produce such effects.

TraziBatine commented 1 year ago

Hi. Thanks for the answer. Here I plot N (statazimut=0) component displacement. On the E component (statazimut=90) such lines are not present, thats why I thought it to be an artefact.

Calculations are done on hourly basis.

Later I will test statazimut=1 to see if it is still present.

Code for plot:

`step = 2.0 grid_lat = np.arange(-90.0, 90.0 + step, step) grid_lon = np.arange(-180.0, 180.0 + step, step)

Construct a dense grid of lon/lat points.

inds_lon, inds_lat = np.meshgrid(grid_lon, grid_lat)

Create an empty array of correct size.

all_vals = np.empty((grid_lat.shape[0], grid_lon.shape[0])) all_vals.fill(np.nan)

iextent = (-180, 180, -90, 90)

inorm = Normalize(vmin=tides_min, vmax=tides_max)

land_color = '#f5f5f3' water_color = '#cdd2d4'

all_invals = list(zip(lons, lats, tides))

Fill empty array with values from input data.

for inval in all_invals: ilon = np.where(inds_lon == inval[0])[1][0] ilat = np.where(inds_lat == inval[1])[0][0] all_vals[ilat, ilon] = inval[2]

fig = plt.figure(figsize=(16,9))

ax = fig.add_subplot(1,1,1, projection=crs.PlateCarree())

ax.add_feature(cfeature.OCEAN, color=water_color, linewidth=1) ax.add_feature(cfeature.LAND, color=land_color, linewidth=1) ax.add_feature(cfeature.COASTLINE, color='black', linewidth=1) ax.gridlines() ax.set_global()

im = ax.imshow(all_vals, extent=iextent, origin='lower', alpha=0.75, interpolation='bilinear', aspect='equal', resample=False, cmap='jet')

fig.tight_layout() plt.show()`

Cheers

hydrogeoscience commented 7 months ago

This seems to be a bug within the Fortran code. I will need some time to dig into this issue ...

hydrogeoscience commented 7 months ago

Just saw that the Fortran code has a comment "Attention: this component has never been tested" - which speaks to this possibly being buggy.