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.
In [ ]:
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
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
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
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:
dx, dy extraído;
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))