(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.
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
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.
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.
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).
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.
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}$$
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.
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.
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."
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.
%pylab --no-import-all tk
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.
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.
%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,dyextraído; - Calcular a corrente geostrófica.
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.
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))
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
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/.