notebook.community

Edit and run 





In [ ]:



    


from IPython.core.display import HTML

with open('../../common/creativecommons.html', 'r') as f:
    html = f.read()
    
with open('../../common/custom.css', 'r') as f:
    styles = f.read()
    
HTML(styles)

text = 'Check this post at'
uri = 'http://nbviewer.ipython.org/urls/raw.github.com/ocefpaf/python4oceanographers/master/content/downloads/notebooks'
name = get_notebook_name()
link = """<p>%s <a href="%s/%s"><em>nbviewer</em>.</a></p>""" % (text, uri, name)
html += str(link)

Explorar dados de altura dinâmica média calculando a corrente geostrófica usando a inclinação da superfície livre.

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.

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.

Imports.

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 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 versa. Importamos também as funções para decompor velocidade e direção em suas componentes u e v e uma para alisar os dados.





In [ ]:



    


from oceans.ff_tools import wrap_lon360, wrap_lon180
from oceans.ff_tools.ocfis import uv2spdir, spdir2uv, 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, y = cube.dim_coords
    X, Y = np.meshgrid(x.points, y.points)
    xi = wrap_lon360(xi)

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

    maskx = np.logical_and(x.points >= min(lon), x.points <= max(lon))
    masky = np.logical_and(y.points >= min(lat), y.points <= 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)

    lons, lats = X[indices], Y[indices]
    elvs = cube.data.T[indices]

    dist, angle = sw.dist(lats, lons, 'km')
    dist *= 1e3
    dist = np.r_[0, dist.cumsum()]
    return (lons, lats), (elvs, dist, angle)

Carregar os dados e "limpá-los."





In [ ]:



    


cube = iris.load_cube('mdt_cnes_cls2009_global_v1.1.nc',
                      iris.Constraint('Mean Dynamic Topography'))
print(cube)



    











In [ ]:



    


data = cube.data.filled(fill_value=np.NaN).copy()
data[data == 9999.0] = np.NaN
data = np.ma.masked_invalid(data)
cube.data = data

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 qt4



    











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))

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=(8, 6))
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, dist, angle) = get_nearest(lon, lat, cube)
elvs = smoo1(elvs, window_len=5, window='hanning')
ix = np.diff(elvs) / np.diff(dist)
v = geostrophic_current(ix, lats.mean())
maximum = ix == ix.max()
dist *= 1e-3  # Converte para km.

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





In [ ]:



    


fig, ax = plt.subplots(figsize=(10, 2))
fig.subplots_adjust(bottom=0.25)
ax.plot(dist, elvs)
ax.axis('tight')
ax.set_ylabel('Height [m]')
ax.set_xlabel('Distance [km]')
ax.set_title('Sea Surface Slope')
vmax = v.max() if v.max() > np.abs(v.min()) else v.min()
symbol = r'$\bigotimes$' if vmax > 0 else r'$\bigodot$'
_ = ax.text(dist[maximum], elvs[maximum], symbol, va='center', ha='center')

_ = ax.annotate(r'%2.2f m s$^{-1}$' % vmax, xy=(dist[maximum], elvs[maximum]),
                xycoords='data', xytext=(-50, 30), textcoords='offset points',
                arrowprops=dict(arrowstyle="->", alpha=0.65,
                                connectionstyle="angle3,angleA=0,angleB=-90"))

Podemos também plota um stick-plot com o perfil do jato de corrente na secção que escolhemos.





In [ ]:



    


fig, ax = plt.subplots(figsize=(6, 6))
ax.set_title('Jet profile')
ax.set_xlabel('Distance [m]')
ax.set_ylabel(r'Velocity [m $^{-1}$]')
xm = mid_point(dist)
kw = dict(scale_units='xy', angles='xy', scale=1)
qk = ax.quiver(xm, [0]*len(xm), [0]*len(v), v, **kw)
_ = ax.set_ylim(fix_axis(v))
_ = ax.set_xlim(fix_axis(xm))

Vamos rodar as corrente calculadas para plotar elas no mapa. (# FIXME!)





In [ ]:



    


rot = angle.mean()
ang, spd = uv2spdir(0, v, rot=rot)
ui, vi = spdir2uv(spd, ang, deg=False)

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={'latitude':lambda cell: ymin <= cell <= ymax,
              'longitude': lambda cell: xmin <= cell <= xmax}
cube = iris.load_cube('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 [ ]:



    


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 = quiverkey(Q, 0.5, 0.05, 0.5, r'0.5 m s${-1}$', fontproperties={'weight': 'bold'})



    











In [ ]:



    


HTML(html)

Content source: ocefpaf/descriptive_oceanography

Similar notebooks:

notebook.community | gallery | about