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

Abaixo importamos o módulo brewer2mpl que nos permite criar e utilizar colormaps usando a excelente ferramenta colorbrewer2.



In [ ]:

    
from brewer2mpl import brewer2mpl
cmap = brewer2mpl.get_map('RdYlGn', 'diverging', 11, reverse=True).mpl_colormap

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 [ ]:

    
from oceans.ff_tools import wrap_lon360, wrap_lon180, smoo1

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=cmap)
    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 tk



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:

Extrair os dados em uma linha definida por esses dois pontos;
Suavizamos as elevações com uma janela do tipo hanning de 5 pontos;
Calcular a inclinação para cada dx, dy extraído;
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))

Podemos carregar apenas os dados nas proximidades da reta para agilizar a nossa figura.



In [ ]:

    
dx = dy = 10
lon = wrap_lon360(lon)
xmin, xmax = map(int, [lon[0]-dx, lon[1]+dx])
ymin, ymax = map(int, (lat[0]-dy, lat[1]+dy))
coord_values=dict(latitude=lambda cell: ymin <= cell <= ymax,
                  longitude=lambda cell: xmin <= cell <= xmax)
cube = iris.load_cube('../../data/mdt_cnes_cls2009_global_v1.1.nc',
                      iris.Constraint(name='Mean Dynamic Topography', coord_values=coord_values))

E finalmente, uma figura com a corrente sobreposta a MDT.



In [ ]:

    
rad = np.deg2rad(angle)
mean_rad = np.arctan2(np.sin(rad).mean(), np.cos(rad).mean())

if mean_rad < 0:
    mean_rad += np.pi
else:
    mean_rad = np.pi - mean_rad

# Gira as correntes antes de plotar.
ui, vi = v * np.cos(mean_rad), v * np.sin(mean_rad)

fig, ax = plot_mdt(cube, projection=ccrs.PlateCarree(), figsize=(10, 10))

kw = dict(marker='o', markeredgecolor='w', linestyle='none', alpha=0.65, markersize=5)
ax.plot(lons, lats, transform=ccrs.PlateCarree(), markerfacecolor='r', **kw)
x, y = map(mid_point, (lons, lats))
kw = dict(color='k', units='inches', alpha=0.65)
Q = ax.quiver(x, y, ui, vi, transform=ccrs.PlateCarree(), **kw)
ax.axis([wrap_lon180(xmin), wrap_lon180(xmax), ymin, ymax])
qk = ax.quiverkey(Q, 0.5, 0.05, 0.5, r'0.5 m s${-1}$', fontproperties={'weight': 'bold'})



In [ ]:

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.

Definição das funções.

Carregar os dados e "limpá-los."

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`.