OPeNDAP Data URLs for TPX08 Data on Geoport THREDDS Server.



In [ ]:

    
cons = ['m2', 'n2', 's2', 'k2', 'k1', 'o1', 'p1', 'q1', 'm4']

url = dict()
for con in cons:
    url[con] = 'http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/tpx08/nc4/%s_nc4.nc' % con

Open the OPeNDAP Dataset using the NetCDF4-Python library.



In [ ]:

    
from netCDF4 import Dataset

nc = Dataset(url['m2']).variables
nc.keys()
print(nc['hRe'])

Specify region to extract [lon_min, lon_max, lat_min, lat_max].



In [ ]:

    
box = [-48.0 + 360, -44.0 + 360, -26.0, -22.0]  # Bacia de Santos.
#box = [-46.5 + 360, -46.2 + 360, -24.2, -23.7]  # Cidade de Santos.
#box = [-54.0 + 360, -30.0 + 360, -35.0,   5.0]  # Brasil



In [ ]:

    
lon, lat = nc['lon_z'][:], nc['lat_z'][:]
bi = (lon >= box[0]) & (lon <= box[1])
bj = (lat >= box[2]) & (lat <= box[3])

Extract the complex amplitude components for this subregion. The index dimensions in NetCDF file are reversed from usual [j,i] ordering, so transpose



In [ ]:

    
zi = nc['hIm'][bi, bj].T
zr = nc['hRe'][bi, bj].T
lon = nc['lon_z'][bi]
lat = nc['lat_z'][bj]

Compute amplitude and phase from real and imaginary parts.



In [ ]:

    
import numpy as np
import numpy.ma as ma

amp = np.abs(zr + 1j * zi) / 1000.  # Convert millibars to meters.
phase = np.rad2deg(np.arctan2(-zi, zr))
amp = ma.masked_equal(amp, 0.)
phase = ma.masked_equal(phase, 0.)



In [ ]:

    
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
from cartopy.feature import LAND
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

cmap = plt.cm.rainbow

cbarkw = dict(shrink=0.75, extend='both')
figkw = dict(ncols=2, figsize=(12, 4), subplot_kw=dict(projection=ccrs.PlateCarree()))
    
fig, (ax0, ax1) = plt.subplots(**figkw)
cs = ax0.pcolormesh(lon, lat, amp, cmap=cmap)
ax0.set_title('TPX08: M2 amplitude (m)')
ax0.coastlines('50m', color='k')
gl = ax0.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                   linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
fig.colorbar(cs, ax=ax0, **cbarkw)

levels = np.arange(-np.pi, np.pi + np.pi/2, np.pi/2)
cs = ax1.pcolormesh(lon, lat, phase, cmap=cmap)
cbar = fig.colorbar(cs, ax=ax1, **cbarkw)
cs = ax1.contour(lon, lat, np.deg2rad(phase), colors='k', levels=levels)
ax1.clabel(cs, fmt='%2.2f')
for t in cs.labelTexts:
    if float(t.get_text()) == -3.14:
        t.set_text(r'$-\pi$')
    elif float(t.get_text()) == -1.57:
        t.set_text(r'$-\pi/2$')
    elif float(t.get_text()) == 0:
        t.set_text(r'$2\pi$')
    elif float(t.get_text()) == 1.57:
        t.set_text(r'$\pi/2$')
    elif float(t.get_text()) == 3.14:
        t.set_text(r'$\pi$')
ax1.set_title('TPX08: M2 phase (degrees)')
ax1.coastlines('50m', color='k')
gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                   linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER



In [ ]:

    
import iris
import iris.plot as iplt
from brewer2mpl import brewer2mpl

constraint = iris.Constraint(coord_values={'longitude': lambda cell: box[0]-360 <= cell <= box[1]-360,
                                           'latitude': lambda cell: box[2] <= cell <= box[3]})

bathy = iris.load_cube('/home/filipe/00-NOBKP/OcFisData/ETOPO1_Bed_g_gmt4.grd', constraint)

land = brewer2mpl.get_map('Greens', 'sequential', 9)
ocean = brewer2mpl.get_map('Blues', 'sequential', 7, reverse=True)

colors = np.array(ocean.mpl_colors + land.mpl_colors[2:])

levels = [-4000, -2500, -1000, -700, -400, -145, -10, 0, 10, 145, 400, 800, 1200, 1600]

fig, ax = plt.subplots(figsize=(8, 8), subplot_kw=dict(projection=ccrs.PlateCarree()))
cs = iplt.contourf(bathy, levels, colors=colors, extend='both')
fig.colorbar(cs, **cbarkw)
ax.coastlines('50m', color='k')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER