python4oceanographers

Turning ripples into waves

Mean Dynamic Topography

(Updating my laboratory lecture on Geostrophic Currents. This post was written in Portuguese.)

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 [2]:
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 [3]:
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 [4]:
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 [5]:
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 [6]:
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 [7]:
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 [8]:
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 [9]:
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 [10]:
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 [11]:
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)
Mean Dynamic Topography / (m)       (longitude: 1440; latitude: 720)
     Dimension coordinates:
          longitude                           x               -
          latitude                            -               x
     Attributes:
          CreatedBy: rio@px-124.cls.fr
          CreatedOn: 11-MAR-2010 16:50:55:000000
          FileType: GRID_DOTS
          OriginalName: MDT_CNES-CLS09_v1.1.nc
          title: MDT CNES-CLS09 V1.1

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

In [12]:
%pylab --no-import-all tk
Populating the interactive namespace from numpy and matplotlib

In [13]:
fig, ax = plot_mdt(cube, projection=ccrs.PlateCarree(), figsize=(12, 10))
_ = ax.set_title('Escolha dois pontos.')
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/coords.py:770: UserWarning: Coordinate 'longitude' is not bounded, guessing contiguous bounds.
  'contiguous bounds.'.format(self.name()))
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/coords.py:770: UserWarning: Coordinate 'latitude' is not bounded, guessing contiguous bounds.
  'contiguous bounds.'.format(self.name()))

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

In [14]:
lon, lat = get_position(fig, ax)
print('Longitude: %s\nLatitude: %s' % (lon, lat))
lon = wrap_lon360(lon)
Longitude: [-57.41000392 -50.95416632]
Latitude: [-39.86813163 -44.26359553]

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/matplotlib/backend_bases.py:2407: MatplotlibDeprecationWarning: Using default event loop until function specific to this GUI is implemented
  warnings.warn(str, mplDeprecation)

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

In [15]:
%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)
Populating the interactive namespace from numpy and matplotlib

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 [16]:
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 [17]:
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))
In [18]:
HTML(html)
Out[18]:

This post was written as an IPython notebook. It is available for download or as a static html.

Creative Commons License
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

Comments