Notebook para calcular a corrente geostrófica entre dois pontos usando o os dados de Mean Dynamic Height do arquivo mdt_cnes_cls2009_global_v1.1.nc.

Esse notebook visa explorar dados de Altura Dinâmica Média calculando a corrente geostrófica gerado pela elevação no nível do mar.

O arquivo mdt_cnes_cls2009_global_v1.1.nc* consiste de uma grade global de Altura Dinâmica Média. Esse conjunto de dados foi criado pelo do grupo AVISO.

* Apesar de ser gratuíto é necessário fazer um pedido através do formulário online explicando como os dados serão utilizados.

Nas células abaixo vamos importar os módulos/funções que utilizaremos. Começando com NumPy e Matplolib que são os nosso pacotes default para trabalhar com m arrays e plotagem.


In [ ]:
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt

Agora vamos import os pacotes iris e cartopy, ambos desenvolvidos pelo UK Met Office e extremamente úteis para carregar, plotar e manipular dados oceanográficos em diversos formatos.


In [ ]:
import iris
import iris.plot as iplt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

Agora importamos as funções de wrap angles para converter os ângulos de longitude do formato -180--180 para 0--360 e vice-e-versa, e para alisar os dados.


In [ ]:
def wrap_lon180(lon):
    lon = np.atleast_1d(lon).copy()
    angles = np.logical_or((lon < -180), (180 < lon))
    lon[angles] = wrap_lon360(lon[angles] + 180) - 180
    return lon


def wrap_lon360(lon):
    lon = np.atleast_1d(lon).copy()
    positive = lon > 0
    lon = lon % 360
    lon[np.logical_and(lon == 0, positive)] = 360
    return lon


def smoo1(datain, window_len=11, window='hanning'):
    r"""Smooth the data using a window with requested size.

    Parameters
    ----------
    datain : array_like
             input series
    window_len : int
                 size of the smoothing window; should be an odd integer
    window : str
             window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'.
             flat window will produce a moving average smoothing.

    Returns
    -------
    data_out : array_like
            smoothed signal
    """

    datain = np.asarray(datain)

    if datain.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")

    if datain.size < window_len:
        raise ValueError("Input vector needs to be bigger than window size.")

    if window_len < 3:
        return datain

    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError("""Window is on of 'flat', 'hanning', 'hamming',
                         'bartlett', 'blackman'""")

    s = np.r_[2 * datain[0] - datain[window_len:1:-1], datain, 2 *
              datain[-1] - datain[-1:-window_len:-1]]

    if window == 'flat':  # Moving average.
        w = np.ones(window_len, 'd')
    else:
        w = eval('np.' + window + '(window_len)')

    data_out = np.convolve(w / w.sum(), s, mode='same')
    return data_out[window_len - 1:-window_len + 1]

Vamos precisar também do módulo seawater EOS-80, mas utilizaremos apenas os cálculos do parâmetro de Coriolis (sw.f), gravidade (sw.g), e distância (sw.dist).


In [ ]:
import seawater as sw

Por último, mas não menos importante, vamos importar o sub-módulo KDTree do módulo scipy. Esse algorítimo nos permite encontrar rapidamento pontos próximos uns dos outros. Utilizaremos ele para achar os dados mais próximos de onde clicarmos com o mouse.


In [ ]:
from scipy.spatial import KDTree

Definição das funções.

Primeiro vamos definir uma função para calcular a velocidade geostrófica em função da inclinação do nível do mar de acordo com a equação:

$$v = i_x \frac{g}{f}$$

In [ ]:
def geostrophic_current(ix, lat):
    g = sw.g(lat.mean())
    f = sw.f(lat.mean())
    v = ix * g / f
    return v

Agora vamos definir algumas funções que vão ajudar a plotar os dados e extrair a informações que precisamos para aplicar a equação acima.


In [ ]:
def fix_axis(lims, p=0.1):
    """Ajusta eixos + ou - p dos dados par exibir melhor os limites."""
    min = lims.min() * (1 - p) if lims.min() > 0 else lims.min() * (1 + p)
    max = lims.max() * (1 + p) if lims.max() > 0 else lims.max() * (1 - p)
    return min, max


def plot_mdt(mdt, projection=ccrs.PlateCarree(), figsize=(12, 10)):
    """Plota 'Mean Dynamic Topography' no mapa global."""
    fig = plt.figure(figsize=figsize)
    ax = plt.axes(projection=projection)
    ax.add_feature(cfeature.LAND, facecolor='0.75')
    cs = iplt.pcolormesh(mdt, cmap=plt.cm.RdYlBu_r)
    ax.coastlines()
    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
    cbar = fig.colorbar(cs, extend='both', orientation='vertical', shrink=0.6)
    cbar.ax.set_title('[m]')
    return fig, ax


def get_position(fig, ax):
    """Escolhe dois pontos para fazer o cálculo."""
    points = np.array(fig.ginput(2))
    lon, lat = points[:, 0], points[:, 1]
    kw = dict(marker='o', markerfacecolor='k', markeredgecolor='w',
              linestyle='none', alpha=0.65, markersize=5)
    ax.plot(lon, lat, transform=ccrs.Geodetic(), **kw)
    ax.set_title('')
    plt.draw()
    return lon, lat


def mid_point(arr):
    return (arr[1:] + arr[:-1]) / 2.

Por último vamos fazer a função que acha os dados na reta definida pelos pontos que escolhemos.


In [ ]:
def get_nearest(xi, yi, cube):
    """Encontra os dados mais próximos aos pontos escolhidos."""
    x = cube.coord('longitude').points
    y = cube.coord('latitude').points
    X, Y = np.meshgrid(x, y)

    tree = KDTree(zip(X.ravel(), Y.ravel()))
    dist, indices = tree.query(np.array([xi, yi]).T)
    indices = np.unravel_index(indices, X.shape)
    lon, lat = X[indices], Y[indices]

    maskx = np.logical_and(x >= min(lon), x <= max(lon))
    masky = np.logical_and(y >= min(lat), y <= max(lat))
    maxnp = len(np.nonzero(maskx)[0]) + len(np.nonzero(masky)[0])

    lons = np.linspace(lon[0], lon[1], maxnp)
    lats = np.linspace(lat[0], lat[1], maxnp)

    # Find all x, y, data in that line using the same KDTree obj.
    dist, indices = tree.query(np.array([lons, lats]).T)
    indices = np.unique(indices)
    indices = np.unravel_index(indices, X.shape)

    X, Y = X[indices], Y[indices]
    elvs = cube.data.T[indices]
    
    # Sort Y with X.
    Y = np.array([y for (x, y) in sorted(zip(X, Y))])

    return X, Y, elvs


def compute_distance(lats, lons):
    dist, angle = sw.dist(lats, lons, 'km')
    dist = np.r_[0, dist.cumsum()]
    return dist, angle

Carregar os dados e "limpá-los."


In [ ]:
cube = iris.load_cube('../../../data/mdt_cnes_cls2009_global_v1.1.nc',
                      iris.Constraint('Mean Dynamic Topography'))
print(cube)

# Coloca uma máscara sobre os dados inválidos.
cube.data = ma.masked_equal(cube.data, 9999.0)

Precisamos plotar os dados em uma figura pop-up para escolher os dois pontos onde vamos calcular a corrente geostrófica.


In [ ]:
%pylab --no-import-all wx

In [ ]:
fig, ax = plot_mdt(cube, projection=ccrs.PlateCarree(), figsize=(12, 10))
_ = ax.set_title('Escolha dois pontos.')

Rode novamente a célula abaixo até escolher seus pontos e quando estiver pronto feche a figura.


In [ ]:
lon, lat = get_position(fig, ax)
print('Longitude: %s\nLatitude: %s' % (lon, lat))
lon = wrap_lon360(lon)

Vamos voltar ao modo "não pop-up" e re-plotar a nossa figura com os pontos que escolhemos.


In [ ]:
%pylab --no-import-all inline
fig, ax = plot_mdt(cube, projection=ccrs.PlateCarree(), figsize=(10, 9))
kw = dict(marker='o', markerfacecolor='k', markeredgecolor='w',
          linestyle='none', alpha=0.65, markersize=5)
_ = ax.plot(lon, lat, transform=ccrs.PlateCarree(), **kw)

Os passos a seguir são:

  1. Extrair os dados em uma linha definida por esses dois pontos;
  2. Suavizamos as elevações com uma janela do tipo hanning de 5 pontos;
  3. Calcular a inclinação para cada dx, dy extraído;
  4. Calcular a corrente geostrófica.

In [ ]:
lons, lats, elvs = get_nearest(lon, lat, cube)

dist, angle = compute_distance(lons, lats)
elvs = smoo1(elvs, window_len=11, window='hanning')

ix = np.diff(elvs) / np.diff(dist * 1e3)
v = geostrophic_current(ix, lats.mean())

Vamos plotar o perfil da inclinação marcando o local de maior gradiente/corrente.


In [ ]:
arrowprops = dict(connectionstyle="angle3,angleA=0,angleB=-90",
                  arrowstyle="->", alpha=0.65)

cdist = mid_point(dist)
idx = np.abs(v).argmax()
maximum = ix == ix.max()
vmax = v[idx]
symbol = r'$\bigotimes$' if vmax > 0 else r'$\bigodot$'


fig, (ax0, ax1) = plt.subplots(nrows=2, figsize=(10, 4), sharex=True)
ax0.plot(dist, elvs)
ax0.axis('tight')
ax0.set_ylabel('Height [m]')

ax0.text(dist[maximum], elvs[maximum], symbol, va='center', ha='center')
ax0.annotate(r'%2.2f m s$^{-1}$' % vmax, xy=(dist[maximum], elvs[maximum]),
             xycoords='data', xytext=(-50, 30), textcoords='offset points',
             arrowprops=arrowprops)
ax0.set_ylim(fix_axis(elvs))
ax0.set_xlim(fix_axis(dist))


ax1.axis('tight')
ax1.set_ylabel(r'Velocity [m s$^{-1}$]')
ax1.set_xlabel('Distance [km]')

kw = dict(scale_units='xy', angles='xy', scale=1)
qk = ax1.quiver(cdist, [0]*len(cdist), [0]*len(v), v, **kw)
ax1.set_ylim(fix_axis(v))
_ = ax1.set_xlim(fix_axis(cdist))