VOLTAR

In [ ]:
# Ao fim deste tutorial você saberá:
#     Como ler dados online com o OpenDap (PyDAP)
#     Como fazer um subset de dados
#     Como tirar flag values de um dado
#     Como calcular e referenciar geopotencial (Método Dinâmico Clássico)
#     Como calular função de corrente (Psi) através do geopotencial
#     Como calcular o campo de velocidade referente a um de Psi
#     Como plotar Psi e campo de velocidade sobre um mapa com batimetria[1]
    
# [1]: Para aprender como fazer o mapa, faça o tutorial "01 - Como fazer um mapa"
        
In [1]:
# De início precisamos importart a função para usar um servidor opendap
# e o pacote padrão de trabalho com matrizes no python
from pydap.client import open_url
import numpy as np

# Para não escrever tudo em uma linha só, quebramos o endereço em duas variáveis
base  = 'http://apdrc.soest.hawaii.edu'
adds = base+'/dods/public_data/WOA/WOA13/0.25_deg/annual/'
# Criamos uma função lamba que lê o dado a partir da entrada da propriedade
WOA13 = lambda prop:open_url(adds+prop)
# para obter temperatura faça WOA13('temp')
# para salinidade faça WOA13('salt')

# Definimos os limites do mapa que iremos plotar
# Ou seja, também são os limites do subset de dados
lati,latf = -20,-7
loni,lonf = -41,-30

# lemos os dados
TDAT,SDAT = WOA13('temp'),WOA13('salt')

# o que será que tem dentro dos dados?
TDAT.keys() # a função keys retorna o nome das variáveis que tem dentro de TDAT
Out[1]:
['tan', 'tmn', 'tdd', 'tsd', 'tse', 'toa', 'tgp', 'time', 'lev', 'lat', 'lon']
In [2]:
# Lemos então as latitudes e longitudes para fazer o subset dos dados
lon,lat = TDAT['lon'][:],TDAT['lat'][:] # [:] seleciona todos os valores

# Criamos as condições, ou seja, criamos variáveis
# booleanas que só retornem True quando estiverem nos
# limites que definimos
condlon = (lon>loni)&(lon<lonf)
condlat = (lat>lati)&(lat<latf)

# Para aplicar o subset precisamos aplicar uma condicional por vez
# Diferente seria se quiséssemos retornar um número definido de pontos
# Aí usaríamos índices ao invés de condidionais e faríamos ao mesmo tempo
TEMP = TDAT['tan'][:,:,condlat,:][:,:,:,condlon]
SALT = SDAT['san'][:,:,condlat,:][:,:,:,condlon]

# O que fizemos foi primeiro selecionar todos os tempos, todos os níveis,
# apenas as latitudes que queremos e todas as longitudes e depois restringimos
# a apenas as longitudes que queremos
In [3]:
# Vamos ver como é esse objeto TEMP
TEMP
Out[3]:
{'tan': <pydap.model.BaseType object at 0x7f8a867ddfd0>, 'time': <pydap.model.BaseType object at 0x7f8a867dde90>, 'lev': <pydap.model.BaseType object at 0x7f8a867dded0>, 'lat': <pydap.model.BaseType object at 0x7f8a86769050>, 'lon': <pydap.model.BaseType object at 0x7f8a86769410>}
In [ ]:
# Ele é um dicionário com vários objetos pydap.model
# Isso quer dizer que se quisermos o valor de cada uma dessas variáveis
# não será o suficiente apenas indexar o nome destas, mas também assinalar que
# queremos todos os seus elementos, ou seja, TEMP['key'][:]
In [4]:
# Vamos então definir todas as coordenadas dessa variável TEMP
lev,tlat,tlon = TEMP['lev'][:],TEMP['lat'][:],TEMP['lon'][:]

# E agora vamos definir a variável TEMPval que conterá os valores
# de temperatura em um formato que possamos plotar ou fazer cálculos
TEMPval = TEMP['tan'][:] #an significa que é a média climatológica

# Vamos ver como estão esses dados no tempo 1 e no nível 1
print TEMPval[0,0,:,:] # lembre que em python a indexação começa em 0
[[  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.57248898e+01
    2.56894913e+01   2.56483097e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.57675896e+01
    2.57245998e+01   2.56776905e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.58025913e+01
    2.57551003e+01   2.57032089e+01]
 ..., 
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.72427902e+01
    2.72305908e+01   2.72166100e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.72485104e+01
    2.72342911e+01   2.72192097e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   2.72522907e+01
    2.72373905e+01   2.72172089e+01]]
In [ ]:
# Parece que o dado está cheio de flags
# Flags são valores absurdos que são colocados para
# indicar alguma coisa:
#     Problema no sensor
#     Local onde não tem como ter dado (temperatura da água em um ponto de continente)
In [5]:
# Para retirar os flags, vamos substituí-los por NaN (Not a Number)
TEMPval[TEMPval==9.96920997e+36] = np.nan

# Vamos ver como ficou
print TEMPval[0,0,:,:]
[[         nan          nan          nan ...,  25.72488976  25.68949127
   25.64830971]
 [         nan          nan          nan ...,  25.76758957  25.72459984
   25.67769051]
 [         nan          nan          nan ...,  25.80259132  25.75510025
   25.70320892]
 ..., 
 [         nan          nan          nan ...,  27.24279022  27.23059082
   27.21660995]
 [         nan          nan          nan ...,  27.24851036  27.23429108
   27.21920967]
 [         nan          nan          nan ...,  27.25229073  27.23739052
   27.21720886]]
In [6]:
# Estamos trabalhando com uma média climatológica,
# então não faz sentido termos a variável tempo
TEMPval.shape
Out[6]:
(1, 102, 52, 44)
In [7]:
# Vamos reduzir a dimensão que só tem 1 unidade dos dados
# Para isso, usaremos a função squeeze
TEMPval = np.squeeze(TEMPval)

# Vamos ver como ficou
TEMPval.shape
Out[7]:
(102, 52, 44)
In [8]:
# Faremos o mesmo para salinidade
SALTval = SALT['san'][:]

# Será que a flag será a mesma
print SALTval[0,0,:,:]
[[  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.71901131e+01
    3.71890106e+01   3.71850128e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.72137108e+01
    3.72108116e+01   3.72051888e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.72350883e+01
    3.72308121e+01   3.72235870e+01]
 ..., 
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.62462883e+01
    3.62411003e+01   3.62360878e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.62205887e+01
    3.62149124e+01   3.62097130e+01]
 [  9.96920997e+36   9.96920997e+36   9.96920997e+36 ...,   3.61966133e+01
    3.61903114e+01   3.61847115e+01]]
In [9]:
# Substituímos
SALTval[SALTval==9.96920997e+36] = np.nan

# Como ficou?
print SALTval[0,0,:,:]
[[         nan          nan          nan ...,  37.19011307  37.18901062
   37.18501282]
 [         nan          nan          nan ...,  37.21371078  37.21081161
   37.20518875]
 [         nan          nan          nan ...,  37.23508835  37.23081207
   37.22358704]
 ..., 
 [         nan          nan          nan ...,  36.2462883   36.24110031
   36.2360878 ]
 [         nan          nan          nan ...,  36.22058868  36.21491241
   36.20971298]
 [         nan          nan          nan ...,  36.19661331  36.19031143
   36.18471146]]
In [10]:
# Da mesma forma vamos reduzir a dimensão tempo
SALTval.shape
Out[10]:
(1, 102, 52, 44)
In [11]:
# Reduzindo
SALTval = np.squeeze(SALTval)

# Vamos ver como ficou
SALTval.shape
# Espera-se que a salinidade e temperatura tenham a mesma dimensão
Out[11]:
(102, 52, 44)
In [12]:
#ignore este primeiro comando
#serve para plotar as imagens neste editor
%matplotlib inline

# Vamos ver como estão esses dados?
# Para isso, vamos ter que importar um módulo para plots
import matplotlib.pyplot as plt

# criamos uma figura
plt.figure()
# plotamos longitudesXlatitudesXtemperatura no nível 0
plt.contourf(tlon,tlat,TEMPval[0,:,:])
# barra de cores
plt.colorbar()

# criamos uma figura
plt.figure()
# plotamos longitudesXlatitudesXsalinidade no nível 0
plt.contourf(tlon,tlat,SALTval[0,:,:])
# barra de cores
plt.colorbar()

# comando para mostrar as figuras
plt.show()

# OBS: Usamos as mesmas latitudes e logitudes da 
#     temperatura para plotar salinidade, isso tem que
#     ser feito com cuidado. Fizemos porque por serem
#     dados climatológicos, devem estar sobre o mesmo grid
In [13]:
# Para plotar função de corrente, vamos precisar calcular 
# o Geopotencial primeiro

# Para isso, importamos o módulo seawater
import seawater as sw

# Vamos ver como funciona essa função
help(sw.gpan)
Help on function gpan in module seawater.geostrophic:

gpan(s, t, p)
    Geopotential Anomaly calculated as the integral of svan from the
    the sea surface to the bottom. THUS RELATIVE TO SEA SURFACE.
    
    Adapted method from Pond and Pickard (p76) to calculate gpan relative to
    sea surface whereas Pond and Pickard calculated relative to the deepest
    common depth.  Note that older literature may use units of "dynamic
    decimeter" for above.
    
    
    Parameters
    ----------
    s(p) : array_like
           salinity [psu (PSS-78)]
    t(p) : array_like
           temperature [℃ (ITS-90)]
    p : array_like
        pressure [db].
    
    Returns
    -------
    gpan : array_like
           geopotential anomaly
           [m :sup:`3` kg :sup:`-1`
           Pa = m :sup:`2` s :sup:`-2` = J kg :sup:`-1`]
    
    Examples
    --------
    >>> # Data from Unesco Tech. Paper in Marine Sci. No. 44, p22.
    >>> import seawater as sw
    >>> s = [[0, 1, 2], [15, 16, 17], [30, 31, 32], [35,35,35]]
    >>> t = [[15]*3]*4
    >>> p = [[0], [250], [500], [1000]]
    >>> sw.gpan(s, t, p)
    array([[   0.        ,    0.        ,    0.        ],
           [  56.35465209,   54.45399428,   52.55961152],
           [  84.67266947,   80.92724333,   77.19028933],
           [ 104.95799186,   99.38799979,   93.82834339]])
    
    References
    ----------
    .. [1] S. Pond & G.Pickard 2nd Edition 1986 Introductory Dynamical
       Oceanography Pergamon Press Sydney. ISBN 0-08-028728-X

In [ ]:
# Parece que precisamos mudar o formato dos nossos dados
# Fazendo testes descobri que essa função não calcula com dados
# em 3 dimensões.No máximo 2 (profXdistância)
# Então temos que juntar as dimensões latitude e logitude em uma só
# antes de calular o geopotencial
In [14]:
# obtendo as dimensões de TEMPval
L,M,N = TEMPval.shape
# Precisamos ter uma propriedade de profundidade
# no mesmo formato dos dados com as dimensões horizontais juntas

# O que eu faço aqui é utilizar de uma propriedade das listas []
# Quando se multiplica uma lista por um valor, a gente junta ela com ela mesma
# o número de valor que você multiplicou. Assim [1,2,3]*2 = [1,2,3,1,2,3]
# Como o objeto lev é da classe numpy.array, quando colocarmos ele dentro de uma lista
# e multiplicarmos por um valor juntará ele com ele mesmo como se fosse em uma nova linha
# [np.array([1,2,3])]*2 = [array([1, 2, 3]), array([1, 2, 3])]
# mas para isso precisamos converter para array essa lista de arrays
# np.array([np.array([1,2,3])]*2) = array([[1, 2, 3],
#                                          [1, 2, 3]])

# Faremos então isso com lev
LEV = np.array([lev]*M*N)
# Vamos ver como está o formato (precisam estar iguais)
print SALTval.reshape(L,M*N).shape
print np.array(LEV).shape
(102, 2288)
(2288, 102)
In [15]:
# teremos que transpor LEV
LEV = LEV.T
print SALTval.reshape(L,M*N).shape
print np.array(LEV).shape
(102, 2288)
(102, 2288)
In [16]:
# Agora podemos calcular o geopotencial
GPAN = sw.gpan(SALTval.reshape(L,M*N),TEMPval.reshape(L,M*N),LEV)
# Não esqueça de mudar o formato do Geopotencial para o formato de TEMPval e SALTval
GPAN = GPAN.reshape(L,M,N)

# Vamos ver como ficou
print GPAN[0,:,:]
[[ nan  nan  nan ...,   0.   0.   0.]
 [ nan  nan  nan ...,   0.   0.   0.]
 [ nan  nan  nan ...,   0.   0.   0.]
 ..., 
 [ nan  nan  nan ...,   0.   0.   0.]
 [ nan  nan  nan ...,   0.   0.   0.]
 [ nan  nan  nan ...,   0.   0.   0.]]
In [17]:
# Só ver os números não é o suficiente, mas
# de acordo com o help da função, espera-se que
# o padrão seja referenciar na superfície

# cria figura
plt.figure()
# plota nível 1 de geopotencial
plt.contourf(tlon,tlat,GPAN[0,:,:])
# barra de cores
plt.colorbar()
plt.show()
In [18]:
# Bom, parece que todos os valores do nível 1 são zero
# Precisamos referenciar em algum nível de não movimento
# coerente com a região

# Eu estudo essa região e apartir de dados de velocidade
# observados, é coerente referenciar o nível de não movimento
# e 1000 metros para representar as correntes acima deste nível


# Vamos ver se existe esse nível nos dados

# cria figura
plt.figure()
# plota profundidades com marcador de bolinha
plt.plot(lev,marker='o')
# modifica os limites em y para algo perto de 1000 metros
plt.ylim([990,1010])
Out[18]:
(990, 1010)
In [19]:
# Parece que tem o nível de 1000 metros

# Vamos ver se é o valor exato
ind = np.where(lev==1000)
print ind
(array([46]),)
In [20]:
# Temos então que o nível de referência será o índice 46

GPANref = GPAN-GPAN[46,:,:] #referenciando

# Não esqueça que a integral do cálculo do geopotencial antes
# começava de cima para baixo e ao referenciar você inverteu o
# sentido. Assim, vamos mudar o sinal de GPANref
GPANref = -GPANref
In [21]:
# Vamos ver agora como ficou o geopotencial para superfície

# cria figura
plt.figure()
# plota nível 1
plt.contourf(tlon,tlat,GPANref[0,:,:])
# barra de cores
plt.colorbar()
# mostra figura
plt.show()
In [22]:
# Através da teoria, sabe-se que a função de corrente
# geostrófica é a razão entre o geopotencial e o parâmetro
# de coriolis

# Assim, vamos calcular a função de corrente geostrófica
# a partir do f médio da região
PSI = GPANref/sw.f(tlat.mean())

# Vamos selecionar um nível para plotar (superfície)
PSIplot = PSI[0,:,:]

# Plotando
# cria figura
plt.figure()
# plota contornos preenchidos de PSIplot
plt.contourf(tlon,tlat,PSIplot)
# barra de cores
plt.colorbar()

# mostra figura
plt.show()
In [23]:
# Os valores de função de corrente são arbitrários
# o que importa são as suas derivadas
# pois definem  o campo de velocidade

# Assim, para facilitar, vamos deixar PSI começar do 0
PSIplot = PSIplot-np.nanmin(PSIplot)

# Para calcular o campo de velocidade através do PSI
# Teremos que importar uma função do OceanLab
from OceanLab.DYN import psi2uv

# A função só funciona se longitude e latitudes estiverem
# em formato de grid, então:
LON,LAT = np.meshgrid(tlon,tlat)
# Calculando U e V
U,V = psi2uv(LON,LAT,PSIplot)

# Vamos ver como ficou
# cria a figura
plt.figure()
# contornos de PSI
plt.contourf(tlon,tlat,PSIplot)
# barra de cores
plt.colorbar()

# Vamos plotar setas para identificar o campo de velocidade
# Para isso vamos usar a função quiver. Entramos com valores
# de longitude, latitude, velocidade zonal e velocidade meridional
Q = plt.quiver(LON,LAT,U,V,scale=1) #scale escalonará o tamanho das setas

# Agora vamos plotar uma identificação do valor das setas
# Isso se chama quiverkey
plt.quiverkey(Q,-39,-10,0.1,label=r'10 cm$\cdot$s$^{-2}$',coordinates='data')

# O primeiro arqumento é o objeto quiver
# O segundo e o terceiro argumentos são a posição
# O quarto argumento é o tamanho (0.1 m/s)
# label informará o que vai ficar escrito
# e coordinates dirá que coordenadas você usará para o 2o e 3o argumentos


plt.show()

# OBS:
#     Quando se coloca r antes de uma string você está informando
#     que o que está dentro deve ser interpretado pelo LaTeX
    
#     As coordenadas podem ser dos dados, ou da própria figura,
#     indo de 0 a 1, começando no canto esquerdo baixo da figura
In [24]:
# Vamos agora plotar isso sobre o estilo de mapa que criamos com o tutorial:
# 01 - Como fazer um mapa

# Modificamos só a resolução do mapa para intermediária ('i')
############################
# Resultado do tutorial:
############################

from mpl_toolkits.basemap import Basemap # módulo de mapas
import cmocean.cm as cm

# Cria uma figura
plt.figure(figsize=(10,10)) #aqui definimos um tamanho maior para a figura
# Cria o mapa
M = Basemap(llcrnrlon=loni, llcrnrlat=lati, 
            urcrnrlon=lonf, urcrnrlat=latf,
            resolution='i',projection='merc')

# Preenchimento dos continentes
M.fillcontinents(color='#BDA973');
# Linhas de costa
M.drawcoastlines(linewidth=1);

# Meridianos e Paralelos
lbs=[1, 0, 0, 1]
M.drawmeridians(range(-180,180,2),labels=lbs);
M.drawparallels(range(-90,90,2),labels=lbs);

# Escala
M.drawmapscale(-39,-9,-35,-14,300,barstyle='fancy',fontsize=12);
############################
############################

# Agora vamos plotar os contornos de PSI
C = M.contourf(LON,LAT,PSIplot*1e-4,50,vmin=0,vmax=3.5,cmap=cm.velocity,latlon=True)
plt.colorbar(C)
# dividimos os valores de psi por 1000 e deixamos o matplotlib escolher 50
# níveis para plotar. Escolhemos o valor mínimo de 0 e máximo de 3.5
# O pacote cmocean não tem uma colormap para função de corrente
# Como os valores de função de corrente são arbitrários, usamos um colormap
# que desse para mostrar a bifurcação que ocorre na região

# Toda vez que plotar algo em um objeto Basemap, não esqueça
# de informar que as coordenadas que vc está usando são
# em latitude e longitude, ou seja, latlon=True

# Aqui faremos o quiver
Q = M.quiver(LON,LAT,U,V,scale=1,latlon=True)

# Aqui vamos fazer o quiverkey, mas definiremos a fonte como negrito
# e em tamnho 15
fontquiver = dict(weight= 'bold',size=15)
plt.quiverkey(Q,0.15,0.7,0.1,label=r'10 cm$\cdot$s$^{-2}$',fontproperties=fontquiver)

# como não definimos as coordenas, por padrão são usadas
# as coordenadas da figura
Out[24]:
<matplotlib.quiver.QuiverKey at 0x7f8a517cff50>
In [ ]:
# Está faltando plotar a batimetria.
# Para baixar a batimetria teremos que importar a função
# para tal do OceanLab

# No tutorial 01 - Como fazer um mapa
# também explica como baixar batimetria
In [25]:
# importando função
from OceanLab.utils import download_bathy

# Baixando batimetria para os limites do mapa
bLON,bLAT,BAT = download_bathy(lnd=loni,lnu=lonf,ltd=lati,ltu=latf)
# Isso costuma demorar um pouco então você só precisa fazer uma vez
# Pois pode salvar esses dados em um formato que você possa ler sempre que precisar
In [26]:
# Para ler e salvar dados no formato pickle do python
# importe estas duas funções
from OceanLab.utils import save_pickle,load_pickle

# Criando dicionário para salvar os dados
BDATA = dict(LON=bLON,LAT=bLAT,BAT=BAT)
# Salvando os dados em uma pasta escolhida com um determinado nome
save_pickle(BDATA,'/home/iury/Downloads/'+'minha_batimetria')
In [27]:
# Agora toda vez que quiser a batimetria deste lugar você usa
# a função load_pickle
BDATA = load_pickle('/home/iury/Downloads/'+'minha_batimetria')
bLON,bLAT,BAT = BDATA['LON'],BDATA['LAT'],BDATA['BAT']

# Vamos ver como estão esses dados
BAT
Out[27]:
masked_array(data =
 [[501.0 648.0 827.0 ..., -4946.0 -4938.0 -4940.0]
 [602.0 670.0 642.0 ..., -4944.0 -4939.0 -4945.0]
 [434.0 623.0 526.0 ..., -4936.0 -4936.0 -4945.0]
 ..., 
 [411.0 348.0 344.0 ..., -5400.0 -5408.0 -5423.0]
 [361.0 306.0 354.0 ..., -5409.0 -5416.0 -5430.0]
 [280.0 291.0 325.0 ..., -5416.0 -5422.0 -5432.0]],
             mask =
 [[False False False ..., False False False]
 [False False False ..., False False False]
 [False False False ..., False False False]
 ..., 
 [False False False ..., False False False]
 [False False False ..., False False False]
 [False False False ..., False False False]],
       fill_value = nan)
In [28]:
###################################################################################
# PARTE REPETIDA DO ÚLTIMO PLOT
###################################################################################
# Cria uma figura
plt.figure(figsize=(10,10)) #aqui definimos um tamanho maior para a figura
# Cria o mapa
M = Basemap(llcrnrlon=loni, llcrnrlat=lati, 
            urcrnrlon=lonf, urcrnrlat=latf,
            resolution='i',projection='merc')

# Preenchimento dos continentes
M.fillcontinents(color='#BDA973');
# Linhas de costa
M.drawcoastlines(linewidth=1);

# Meridianos e Paralelos
lbs=[1, 0, 0, 1]
M.drawmeridians(range(-180,180,2),labels=lbs);
M.drawparallels(range(-90,90,2),labels=lbs);

# Escala
M.drawmapscale(-39,-9,-35,-14,300,barstyle='fancy',fontsize=12);

# Contornos de PSI
C = M.contourf(LON,LAT,PSIplot*1e-4,50,vmin=0,vmax=3.5,cmap=cm.velocity,latlon=True)
# barra de cores
cbar = plt.colorbar(C)

# campo de velocidade
Q = M.quiver(LON,LAT,U,V,scale=1,latlon=True)
# Quiverkey
fontquiver = dict(weight= 'bold',size=15)
plt.quiverkey(Q,0.15,0.7,0.1,label=r'10 cm$\cdot$s$^{-2}$',fontproperties=fontquiver)
#######################################################################################

# Agora vamos plotar a batimetria
# Para isso faremos primeiro um contorno preenchido da profundidade escolhida
# até o zero com uma cor qualquer
# Eu escolhi uma cor na escala cinza. Você coloca entre aspas um valor
# de 0 (preto) até 1 (branco)
M.contourf(bLON,bLAT,BAT,[-1200,0],colors='0.6',latlon=True)
# Agora vamos plotar uma linha nesse nível de 1200 metros
M.contour(bLON,bLAT,BAT,[-1200],colors='k',latlon=True)
Out[28]:
<matplotlib.contour.QuadContourSet at 0x7f8a3ef77090>