In [ ]:

    

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from oceans.datasets import etopo_subset


from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from mpl_toolkits.axes_grid1.inset_locator import inset_axes


    

In [ ]:

    

def map_limits(m):
    llcrnrlon = min(m.boundarylons)
    urcrnrlon = max(m.boundarylons)
    llcrnrlat = min(m.boundarylats)
    urcrnrlat = max(m.boundarylats)
    return llcrnrlon, urcrnrlon, llcrnrlat, urcrnrlat


def make_map(llcrnrlon=-42, urcrnrlon=-33, llcrnrlat=-21.5,
             urcrnrlat=-15.5, projection='merc', resolution='i',
             figsize=(6, 6), inset=True):
    m = Basemap(llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon,
                llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat,
                projection=projection, resolution=resolution)
    fig, ax = plt.subplots(figsize=figsize)
    m.drawstates()
    m.drawcoastlines()
    m.fillcontinents(color='0.85')
    meridians = np.arange(llcrnrlon, urcrnrlon + 2, 2)
    parallels = np.arange(llcrnrlat, urcrnrlat + 1, 1)
    m.drawparallels(parallels, linewidth=0, labels=[1, 0, 0, 0])
    m.drawmeridians(meridians, linewidth=0, labels=[0, 0, 0, 1])
    m.llcrnrlon = llcrnrlon
    m.urcrnrlon = urcrnrlon
    m.llcrnrlat = llcrnrlat
    m.urcrnrlat = urcrnrlat
    m.ax = ax

    if inset:
        axin = inset_axes(m.ax, width="40%", height="40%", loc=5)
        # Global inset map.
        inmap = Basemap(projection='ortho', lon_0=-44.5, lat_0=-25.5,
                        ax=axin, anchor='NE')
        inmap.drawcountries(color='white')
        inmap.fillcontinents(color='gray')
        bx, by = inmap(m.boundarylons, m.boundarylats)
        xy = list(zip(bx, by))
        mapboundary = Polygon(xy, edgecolor='k', linewidth=1, fill=False)
        inmap.ax.add_patch(mapboundary)
    return fig, m


    

In [ ]:

    

pos = '../../data/abrolhos/posicoes_abrolhos.dat'
latd, latm, lond, lonm, depth, station = np.loadtxt(pos, unpack=True)

lat = -(latd + latm/60.) 
lon = -(lond + lonm/60.)


    

In [ ]:

    

fig, m = make_map(llcrnrlon=-42, urcrnrlon=-33,
                  llcrnrlat=-21.5, urcrnrlat=-15.5, figsize=(7, 5))

# Topography.
topo = '../../data/ETOPO1_Bed_g_gmt4.grd'
x, y, topo = etopo_subset(*map_limits(m), smoo=True, tfile=topo)
topo = np.where(topo > -1., 1.e10, topo)
topo = ma.masked_values(topo, 1.e10)

# Coutour topography.
levels = [100, 1500, 2000, 3000, 5000]
cs = m.contour(x, y, -topo, levels, colors='k',
               latlon=True, alpha=0.5)

m.ax.clabel(cs, fmt='%1.0f m', fontsize=8, inline=1)

# Plot stations.
kw = dict(marker='.', markerfacecolor='r', markeredgecolor='k',
          linestyle='none', latlon=True)
_ = m.plot(lon, lat, **kw)


    

In [ ]:

    

import gsw

dist = gsw.distance(lon, lat)
dist = np.r_[0, dist.cumsum()] / 1e3


    

In [ ]:

    

fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(dist, depth, 'k.-')
ax.plot(dist[2:9], depth[2:9], 'r.')
ax.plot(dist[10:20], depth[10:20], 'r.')
ax.set_ylim(-200, max(depth))
ax.fill_betweenx(depth, dist, color='lightblue', alpha=0.5)
ax.invert_yaxis()
ax.set_ylabel('Depth [m]')
ax.set_xlabel('Distance [km]')


    

In [ ]:

    

import os
from glob import glob
from collections import OrderedDict

from pandas import Panel
from ctd import DataFrame, Series 

def proc_ctd(fname):
    """Quick-n-dirty CTD processing."""
    cast = DataFrame.from_fsi(fname, below_water=True,
                              compression='bz2').split()[0]
    # Removed unwanted columns.
    keep = set(['POT*', 'SAL*'])
    drop = keep.symmetric_difference(cast.columns)
    cast.drop(drop, axis=1, inplace=True)
    cast = cast.apply(Series.bindata, **dict(delta=1.))
    cast = cast.apply(Series.interpolate)
    cast.name = name
    return cast

section = OrderedDict()
for fname in sorted(glob('../../data/abrolhos/*.txt.bz2')):
    name = os.path.basename(fname).split('.')[0]
    cast = proc_ctd(fname)
    section.update({name: cast})
    
section = Panel.fromDict(section)
SAL = section.minor_xs('SAL*')
TEMP = section.minor_xs('POT*')

SAL.lon, SAL.lat = lon, lat
TEMP.lon, TEMP.lat = lon, lat


    

In [ ]:

    

fig, ax, cb = TEMP.plot_section(figsize=(8, 4))
fig, ax, cb = SAL.plot_section(figsize=(8, 4))


    

In [ ]:

    

p = SAL.index.values[..., None]
N2, p_mid = gsw.Nsquared(SAL.values, TEMP.values, p)
N = np.sqrt(N2)
N = ma.masked_invalid(N)


    

In [ ]:

    

g = 9.8

def Rd(lat0, h):
    """Raio de deformação de Rossby."""
    return np.sqrt(g*h) / abs(gsw.f(lat0))

def Rdi(lat0, N):
    """Compute the first baroclinic Rossby
    radius after Chelton et al. 1998 J.Phys.Oc."""
    # Compute phase speed, equation 2.2
    C = N / np.pi

    # Compute radius.
    return C / np.abs(gsw.f(lat0))


    

In [ ]:

    

#N0 = N.sum()
N0 = N[600:-1].sum()
lat0 = lat.mean()

h = max(depth[2:9])

R = dist[8] - dist[2]
Rtrop = Rd(lat0, h) / 1e3
Rclim = Rdi(lat0, N0) / 1e3

print("R:    %i\nRdi: %i\nRd: %i" % (R, Rclim, Rtrop))


    

In [ ]:

    

h = max(depth[10:20])

R = dist[19] - dist[10]
Rtrop = Rd(lat0, h) / 1e3
Rclim = Rdi(lat0, N0) / 1e3

print("R:   %i\nRdi: %i\nRd: %i" % (R, Rclim, Rtrop))


    

In [ ]:

    

fig, ax = plt.subplots(figsize=(3, 4))
ax.plot(N.sum(axis=1), p_mid[..., 0])
ax.invert_yaxis()


    

In [ ]: