#!/usr/bin/env python
# coding: utf-8

# **Course website**: http://www.leouieda.com/geofisica1
# 
# **Note**: This notebook is part of the course "Geofísica 1" of Geology program of the 
# [Universidade do Estado do Rio de Janeiro](http://www.uerj.br/). 
# All content can be freely used and adapted under the terms of the 
# [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).
# 
# ![Creative Commons License](https://i.creativecommons.org/l/by/4.0/88x31.png)

# Esse documento que você está usando é um [Jupyter notebook](http://jupyter.org/). É um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (números, texto, figuras, videos, etc).

# # Prática 4 - Gravimetria - Estudos de caso da anomalia da gravidade

# ## Objetivos
# 
# * Observar a topografia, distúrbio e anomalia Bouguer de: Havaí, Baía de Hudson, Japão e costa brasileira.
# * Fixar os conceitos de distúrbio da gravidade, anomalia Bouguer e isostasia.
# * Fazer gráficos de perfis cortando as principais feições de cada alvo.
# * Analisar a relação entre a topografia, isostasia, anomalias e tectônica.

# ## Instruções
# 
# O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.
# 
# As células com números ao lado, como `In [1]:`, são código [Python](http://python.org/). Algumas dessas células não produzem resultado e servem de preparação para os exemplos interativos. Outras, produzem gráficos interativos. **Você deve executar todas as células, uma de cada vez**, mesmo as que não produzem gráficos.
# 
# Para executar uma célula, clique em cima dela e aperte `Shift + Enter`. O foco (contorno verde ou cinza em torno da célula) deverá passar para a célula abaixo. Para rodá-la, aperte `Shift + Enter` novamente e assim por diante. Você pode executar células de texto que não acontecerá nada.

# ## Preparação
# 
# Exectute as células abaixo para carregar as componentes necessárias para nossa prática. Vamos utilizar várias *bibliotecas*, inclusive uma de geofísica chamada [Fatiando a Terra](http://www.fatiando.org).

# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import ipywidgets as widgets
from IPython.display import display, HTML
import seaborn
import fatiando
from fatiando.gravmag import normal_gravity
from fatiando import gridder
from icgem import load_icgem_gdf


# In[ ]:


print("Usando a versão do Fatiando: {}".format(fatiando.__version__))


# In[ ]:


def data_minmax(data):
    """Retorna vmin e vmax para que o valor 0 esteja no centro da escala de cor."""
    ranges = np.abs([data.min(), data.max()]).max()
    return dict(vmin=-ranges, vmax=ranges)


# # Havaí
# 
# Nosso primeiro alvo é o Havaí:
# 

# In[ ]:


HTML('<iframe src="https://www.google.com/maps/embed?pb=!1m10!1m8!1m3!1d5916036.911229154!2d-157.7151201!3d20.5932929!3m2!1i1024!2i768!4f13.1!5e1!3m2!1sen!2sbr!4v1413401199263" width="600" height="450" frameborder="0" style="border:0"></iframe>')


# Primeiro precisamos carregar os dados de topografia e gravidade do Havaí que estão nos arquivos `data/eigen-6c3stat-havai.gdf` e `data/etopo1-havai.gdf`.

# In[ ]:


havai = load_icgem_gdf('../data/eigen-6c3stat-havai.gdf')
havai['topo']  = load_icgem_gdf('../data/etopo1-havai.gdf', usecols=[-1])['topography_shm']


# Agora posso calcular o distúrbio e a anomalia Bouguer.

# In[ ]:


havai['disturbio'] = havai['gravity_earth'] - normal_gravity.gamma_closed_form(havai['latitude'], havai['h_over_geoid'])
havai['bouguer'] = havai['disturbio'] - normal_gravity.bouguer_plate(havai['topo'], density_rock=2670, density_water=1040)


# Vamos os mapas das anomalias e da topografia  usando a [projeção Mercator](http://en.wikipedia.org/wiki/Mercator_projection) e escalas de cor apropriadas.

# In[ ]:


s, n, w, e = havai['area']
havai_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, 
                   lat_ts=0.5*(s + n), resolution='h')


# In[ ]:


x, y = havai_bm(havai['longitude'], havai['latitude'])
fig, axes = plt.subplots(1, 3, figsize=(14, 9))
ax = axes[0]
tmp = havai_bm.contourf(x, y, havai['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                        **data_minmax(havai['topo']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')
havai_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
ax = axes[1]
tmp = havai_bm.contourf(x, y, havai['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(havai['disturbio']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
havai_bm.drawcoastlines(ax=ax)
ax.set_title('Disturbio')
ax = axes[2]
tmp = havai_bm.contourf(x, y, havai['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(havai['bouguer']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
havai_bm.drawcoastlines(ax=ax)
ax.set_title('Anomalia Bouguer')
plt.tight_layout(pad=0, w_pad=0)


# Agora que temos os mapas, vamos ver como é o gráfico da **anomalia Bouguer pela topografia equivalente**.
# 
# Lembrando que a topografia equivalente é:
# 
# $$
# h_{eq} = 
# \begin{cases}
#     h & \quad \text{if } h  \ge 0 \\
#     h \frac{\rho_c - \rho_a}{\rho_c} & \quad \text{if } h < 0 \\
# \end{cases}
# $$
# 
# em que $h$ é a topografia, $\rho_c$ é a densidade da crosta e $\rho_a$ é a densidade da água do mar.

# In[ ]:


heq = havai['topo'].copy()
heq[heq < 0] *= (2670 - 1040)/2670
plt.figure(figsize=(7, 7))
plt.title('Bouguer x Topografia - Havai')
plt.plot(heq, havai['bouguer'], '.k')
plt.xlabel('Topografia equivalente (m)')
plt.ylabel('Anomalia Bouguer (mGal)')
plt.tight_layout(pad=0)


# É muito útil olhar os dados alongo de um **perfil** que corte as principais feições dos mapas acima. Vamos estrair o perfil entre os pontos `p1` e `p2` definidos abaixo.

# In[ ]:


p1, p2 = havai_bm(197.5, 15), havai_bm(207, 27.5)

plt.figure(figsize=(8, 7))
ax = plt.subplot(111)
tmp = havai_bm.contourf(x, y, havai['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                        **data_minmax(havai['topo']))
plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')
havai_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
xp, yp = np.transpose([p1, p2])
havai_bm.plot(xp, yp, 'o-k', linewidth=2)
plt.tight_layout(pad=0)


# Rode a célula abaixo para extrair o perfil dos dados.

# In[ ]:


tmp = gridder.profile(x, y, havai['topo'], p1, p2, 100)
perfil_havai = dict(distancia=tmp[2], topo=tmp[3])
perfil_havai['bouguer'] = gridder.profile(x, y, havai['bouguer'], p1, p2, 100)[-1]
perfil_havai['disturbio'] = gridder.profile(x, y, havai['disturbio'], p1, p2, 100)[-1]


# Agora podemos fazer um gráfico com os 3 dados (topografia e anomalias) ao longo do perfil.

# In[ ]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))
ax1, ax2 = axes
d = perfil_havai['distancia']*0.001
ax2.fill_between([d.min(), d.max()], [0, 0], -6000, color='#2780E3')
ax2.fill_between(d, perfil_havai['topo'], -6000, color='#333333')
ax2.set_ylabel('Topografia (m)')
ax2.set_xlabel('Distancia (km)')

ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')
ax1.plot(d, perfil_havai['bouguer'], '-', label='Bouguer')
ax1.plot(d, perfil_havai['disturbio'], '-', label='Disturbio')
ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')
ax1.legend(loc='center right')
ax1.set_xlim(d.min(), d.max())
plt.tight_layout(h_pad=0, pad=0)


# ### Para pensar
# 
# * As ilhas estão em equilíbrio isostático?
# * Por que o distúrbio é positivo em cima das ilhas?
# * Por que o distúrbio tem seu mínimo logo ao lado das ilhas?
# * Olhando para o perfil de topografia, por que após o baixo ao lado da ilha há um alto na topografia antes de estabilizar entre -5000 m e -6000 m?
# * Como e por que as feições do distúrbio se correlacionam com a topografia?
# * Por que a anomalia Bouguer é sempre positiva?
# * Por que a anomalia Bouguer é menos positiva sobre as ilhas?
# 

# # Japão

# Vamos repetir o mesmo que fizemos antes mas para  os dados do Japão:

# In[ ]:


HTML('<iframe src="https://www.google.com/maps/embed?pb=!1m10!1m8!1m3!1d10208277.426130258!2d134.3822639!3d36.1346696!3m2!1i1024!2i768!4f13.1!5e1!3m2!1sen!2sbr!4v1413403088567" width="600" height="450" frameborder="0" style="border:0"></iframe>')


# Primeiro, carregar os dados e calcular o distúrbio e anomalia Bouguer.

# In[ ]:


japao = load_icgem_gdf('../data/eigen-6c3stat-japao.gdf')
japao['topo']  = load_icgem_gdf('../data/etopo1-japao.gdf', usecols=[-1])['topography_shm']
japao['disturbio'] = japao['gravity_earth'] - normal_gravity.gamma_closed_form(japao['latitude'], japao['h_over_geoid'])
japao['bouguer'] = japao['disturbio'] - normal_gravity.bouguer_plate(japao['topo'], density_rock=2670, density_water=1040)


# Vamos ver como ficam as anomalias e a topografia para essa região.

# In[ ]:


s, n, w, e = japao['area']
japao_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, 
                   lat_ts=0.5*(s + n), resolution='l')


# In[ ]:


x, y = japao_bm(japao['longitude'], japao['latitude'])
fig, axes = plt.subplots(1, 3, figsize=(15, 10))
ax = axes[0]
tmp = japao_bm.contourf(x, y, japao['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                        **data_minmax(japao['topo']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')
japao_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
ax = axes[1]
tmp = japao_bm.contourf(x, y, japao['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(japao['disturbio']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
japao_bm.drawcoastlines(ax=ax)
ax.set_title('Disturbio')
ax = axes[2]
tmp = japao_bm.contourf(x, y, japao['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(japao['bouguer']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
japao_bm.drawcoastlines(ax=ax)
ax.set_title('Anomalia Bouguer')
plt.tight_layout(pad=0, w_pad=0)


# Agora podemos fazer o gráfico de anomalia Bouguer por topografia equivalente.

# In[ ]:


heq = japao['topo'].copy()
heq[heq < 0] *= (2670 - 1040)/2670
plt.figure(figsize=(7, 7))
plt.title('Bouguer x Topografia - Japao')
plt.plot(heq, japao['bouguer'], '.k')
plt.xlabel('Topografia equivalente (m)')
plt.ylabel('Anomalia Bouguer (mGal)')
plt.tight_layout(pad=0)


# Vamos extrair um perfil cortando a zona de subducção, ilhas do Japão, Mar do Japão e chegando no continente.

# In[ ]:


p1, p2 = japao_bm(130, 43), japao_bm(148, 38)

plt.figure(figsize=(9, 5))
ax = plt.subplot(111)
tmp = japao_bm.contourf(x, y, japao['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                        **data_minmax(japao['topo']))
plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')
japao_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
xp, yp = np.transpose([p1, p2])
japao_bm.plot(xp, yp, 'o-k', linewidth=2)
plt.tight_layout(pad=0)


# In[ ]:


tmp = gridder.profile(x, y, japao['topo'], p1, p2, 100)
perfil_japao = dict(distancia=tmp[2], topo=tmp[3])
perfil_japao['bouguer'] = gridder.profile(x, y, japao['bouguer'], p1, p2, 100)[-1]
perfil_japao['disturbio'] = gridder.profile(x, y, japao['disturbio'], p1, p2, 100)[-1]


# In[ ]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))
ax1, ax2 = axes
d = perfil_japao['distancia']*0.001
ax2.fill_between([d.min(), d.max()], [0, 0], -10000, color='#2780E3')
ax2.fill_between(d, perfil_japao['topo'], -10000, color='#333333')
ax2.set_ylabel('Topografia (m)')
ax2.set_xlabel('Distancia (km)')

ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')
ax1.plot(d, perfil_japao['bouguer'], '-', label='Bouguer')
ax1.plot(d, perfil_japao['disturbio'], '-', label='Disturbio')
ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')
ax1.legend(loc='center right')
ax1.set_xlim(d.min(), d.max())
plt.tight_layout(h_pad=0, pad=0)


# ### Para pensar
# 
# * Por que o distúrbio é praticamente zero na placa do Pacífico e no Mar do Japão?
# * Por que o distúrbio apresenta um par negativo-positivo quando passamos da crosta oceânica para as ilhas do Japão?
# * Qual é a relação entre o distúrbio, a topografia e a subducção?
# 

# # Costa brasileira
# 
# Vamos carregar os dados e fazer a análise para a margem leste brasileira (um lugar que esperamos estar em equilíbrio isostático).

# In[ ]:


HTML('<iframe src="https://www.google.com/maps/embed?pb=!1m14!1m12!1m3!1d9620784.768866993!2d-50.015484842383465!3d-26.80125719375305!2m3!1f0!2f0!3f0!3m2!1i1024!2i768!4f13.1!5e1!3m2!1sen!2sus!4v1443710325602" width="600" height="450" frameborder="0" style="border:0" allowfullscreen></iframe>')


# In[ ]:


brasil = load_icgem_gdf('../data/eigen-6c3stat-brasil.gdf')
brasil['topo']  = load_icgem_gdf('../data/etopo1-brasil.gdf', usecols=[-1])['topography_shm']
brasil['disturbio'] = brasil['gravity_earth'] - normal_gravity.gamma_closed_form(brasil['latitude'], brasil['h_over_geoid'])
brasil['bouguer'] = brasil['disturbio'] - normal_gravity.bouguer_plate(brasil['topo'], density_rock=2670, density_water=1040)


# In[ ]:


s, n, w, e = brasil['area']
brasil_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, 
                    lat_ts=0.5*(s + n), resolution='l')


# In[ ]:


x, y = brasil_bm(brasil['longitude'], brasil['latitude'])
fig, axes = plt.subplots(1, 3, figsize=(15, 10))
ax = axes[0]
tmp = brasil_bm.contourf(x, y, brasil['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                        **data_minmax(brasil['topo']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Topografia')
ax = axes[1]
tmp = brasil_bm.contourf(x, y, brasil['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(brasil['disturbio']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Disturbio')
ax = axes[2]
tmp = brasil_bm.contourf(x, y, brasil['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(brasil['bouguer']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Anomalia Bouguer')
plt.tight_layout(pad=0, w_pad=0)


# In[ ]:


heq = brasil['topo'].copy()
heq[heq < 0] *= (2670 - 1040)/2670
plt.figure(figsize=(7, 7))
plt.title('Bouguer x Topografia - Brasil')
plt.plot(heq, brasil['bouguer'], '.k')
plt.xlabel('Topografia equivalente (m)')
plt.ylabel('Anomalia Bouguer (mGal)')
plt.tight_layout(pad=0)


# In[ ]:


p1, p2 = brasil_bm(360 - 58.5, -20.5), brasil_bm(360 - 33.5, -31)

plt.figure(figsize=(9, 5))
ax = plt.subplot(111)
tmp = brasil_bm.contourf(x, y, brasil['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                         **data_minmax(brasil['topo']))
plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')
brasil_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
xp, yp = np.transpose([p1, p2])
brasil_bm.plot(xp, yp, 'o-k', linewidth=2)
plt.tight_layout(pad=0)


# In[ ]:


tmp = gridder.profile(x, y, brasil['topo'], p1, p2, 200)
perfil_brasil = dict(distancia=tmp[2], topo=tmp[3])
perfil_brasil['bouguer'] = gridder.profile(x, y, brasil['bouguer'], p1, p2, 200)[-1]
perfil_brasil['disturbio'] = gridder.profile(x, y, brasil['disturbio'], p1, p2, 200)[-1]


# In[ ]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))
ax1, ax2 = axes
d = perfil_brasil['distancia']*0.001
ax2.fill_between([d.min(), d.max()], [0, 0], -6000, color='#2780E3')
ax2.fill_between(d, perfil_brasil['topo'], -6000, color='#333333')
ax2.set_ylabel('Topografia (m)')
ax2.set_xlabel('Distancia (km)')

ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')
ax1.plot(d, perfil_brasil['bouguer'], '-', label='Bouguer')
ax1.plot(d, perfil_brasil['disturbio'], '-', label='Disturbio')
ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')
ax1.legend(loc='center right')
ax1.set_xlim(d.min(), d.max())
plt.tight_layout(h_pad=0, pad=0)


# Vamos tentar remover o efeito da Moho isostática de nossos dados. Para fazer isso, vamos aproveitar da relação linear entre a anomalia Bouguer e a topografia equivalente. Se conseguirmos determinar a reta que melhor ajusta nossos dados, podemos utilizar essa reta como sendo o efeito gravitacional da Moho. Ao subtrairmos os dados que caem na reta dos dados obserados, estamos retirando o efeito da Moho isostática.

# In[ ]:


heq = brasil['topo'].copy()
heq[heq < 0] *= (2670 - 1040)/2670
A = np.ones((heq.size, 2))
A[:, 0] = heq
p = np.linalg.solve(A.T.dot(A), A.T.dot(brasil['bouguer']))
efeito_moho = A.dot(p)


# In[ ]:


plt.figure(figsize=(7, 7))
plt.title('Bouguer x Topografia - Brasil')
plt.plot(heq, brasil['bouguer'], '.k', label='Dados observados')
plt.plot(heq, efeito_moho, '-', label='y = {:.3f} heq + {:.3f}'.format(*p))
plt.legend(loc='upper right')
plt.xlabel('Topografia equivalente (m)')
plt.ylabel('Anomalia Bouguer (mGal)')
plt.tight_layout(pad=0)


# Agora que determinamos o efeito da Moho isostática, podemos removê-lo da nossa anomalia Bouguer.

# In[ ]:


x, y = brasil_bm(brasil['longitude'], brasil['latitude'])
fig, axes = plt.subplots(1, 3, figsize=(15, 11))
ax = axes[0]
tmp = brasil_bm.contourf(x, y, brasil['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(brasil['disturbio']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Disturbio')
ax = axes[1]
tmp = brasil_bm.contourf(x, y, brasil['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                        **data_minmax(brasil['bouguer']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Anomalia Bouguer')
ax = axes[2]
tmp = brasil_bm.contourf(x, y, brasil['bouguer'] - efeito_moho, 80, tri=True, ax=ax, cmap='RdBu_r', 
                        vmin=-110, vmax=110)
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
brasil_bm.drawcoastlines(ax=ax)
brasil_bm.drawcountries(ax=ax)
ax.set_title('Anomalia Bouguer - Moho')
plt.tight_layout(pad=0, w_pad=0)


# ### Para pensar
# 
# * A margem brasileira está em equilíbrio isostático?
# * O que estamos presumindo quando removemos o efeito da Moho utilizando o método acima?
# * O que sobra na anomalia Bouguer após removermos o efeito da Moho isostática?
# * Qual é origem da anomalia residual negativa que vemos na porção continental?
# * Qual é origem das anomalias residuais negativas que vemos na porção oceânica?

# # Hudson Bay

# Vamos olhar como ficam as anomalias na região da Baía de Hudson, no Norte do Canadá:

# In[ ]:


HTML('<iframe src="https://www.google.com/maps/embed?pb=!1m10!1m8!1m3!1d12110458.266253604!2d-85.99188428241526!3d61.3758016298575!3m2!1i1024!2i768!4f13.1!5e1!3m2!1sen!2sbr!4v1413404706398" width="600" height="450" frameborder="0" style="border:0"></iframe>')


# In[ ]:


hudson = load_icgem_gdf('../data/eigen-6c3stat-hudson.gdf')
hudson['topo']  = load_icgem_gdf('../data/etopo1-hudson.gdf', usecols=[-1])['topography_shm']
hudson['disturbio'] = hudson['gravity_earth'] - normal_gravity.gamma_closed_form(hudson['latitude'], hudson['h_over_geoid'])
hudson['bouguer'] = hudson['disturbio'] - normal_gravity.bouguer_plate(hudson['topo'], density_rock=2670, density_water=1040)


# Dessa vez, vamos usar uma [projeção polar esterográfica](http://en.wikipedia.org/wiki/Stereographic_projection) para fazer os mapas. 
# A projeção Mercator distorce demais as regiões de altas latitudes.
# Nesses casos, as [projeções cônicas](http://en.wikipedia.org/wiki/Map_projection#Conic) são melhores.

# In[ ]:


s, n, w, e = hudson['area']
hudson_bm = Basemap(projection='laea', width=4000000, height=3500000,
                   lat_ts=0.5*(s + n), lat_0=0.5*(s + n), lon_0=0.5*(e  + w), resolution='l')


# In[ ]:


x, y = hudson_bm(hudson['longitude'], hudson['latitude'])
fig, axes = plt.subplots(1, 3, figsize=(15, 12))
ax = axes[0]
tmp = hudson_bm.contourf(x, y, hudson['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                         **data_minmax(hudson['topo']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')
hudson_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
ax = axes[1]
tmp = hudson_bm.contourf(x, y, hudson['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                         **data_minmax(hudson['disturbio']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
hudson_bm.drawcoastlines(ax=ax)
ax.set_title('Disturbio')
ax = axes[2]
tmp = hudson_bm.contourf(x, y, hudson['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', 
                         **data_minmax(hudson['bouguer']))
plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')
hudson_bm.drawcoastlines(ax=ax)
ax.set_title('Anomalia Bouguer')
plt.tight_layout(pad=0, w_pad=0)


# In[ ]:


heq = hudson['topo'].copy()
heq[heq < 0] *= (2670 - 1040)/2670
plt.figure(figsize=(7, 7))
plt.title('Bouguer x Topografia - Bahia de Hudson')
plt.plot(heq, hudson['bouguer'], '.k')
plt.xlabel('Topografia equivalente (m)')
plt.ylabel('Anomalia Bouguer (mGal)')
plt.tight_layout(pad=0)


# O perfil que vamos extrair vai do continente ao mar do Norte, cortando a baía.

# In[ ]:


p1, p2 = hudson_bm(260.5, 45.5), hudson_bm(299, 70)

plt.figure(figsize=(9, 5))
ax = plt.subplot(111)
tmp = hudson_bm.contourf(x, y, hudson['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', 
                         **data_minmax(hudson['topo']))
plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')
hudson_bm.drawcoastlines(ax=ax)
ax.set_title('Topografia')
xp, yp = np.transpose([p1, p2])
hudson_bm.plot(xp, yp, 'o-k', linewidth=2)
plt.tight_layout(pad=0)


# In[ ]:


tmp = gridder.profile(x, y, hudson['topo'], p1, p2, 200)
perfil_hudson = dict(distancia=tmp[2], topo=tmp[3])
perfil_hudson['bouguer'] = gridder.profile(x, y, hudson['bouguer'], p1, p2, 200)[-1]
perfil_hudson['disturbio'] = gridder.profile(x, y, hudson['disturbio'], p1, p2, 200)[-1]


# In[ ]:


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))
ax1, ax2 = axes
d = perfil_hudson['distancia']*0.001
ax2.fill_between([d.min(), d.max()], [0, 0], -2000, color='#2780E3')
ax2.fill_between(d, perfil_hudson['topo'], -2000, color='#333333')
ax2.set_ylabel('Topografia (m)')
ax2.set_xlabel('Distancia (km)')

ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')
ax1.plot(d, perfil_hudson['bouguer'], '-', label='Bouguer')
ax1.plot(d, perfil_hudson['disturbio'], '-', label='Disturbio')
ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')
ax1.legend(loc='center right')
ax1.set_xlim(d.min(), d.max())
plt.tight_layout(h_pad=0, pad=0)


# ### Para pensar
# 
# * Esta região está em equilíbrio isostático?
# * Por que o distúrbio fica cada vez mais negativo quando mais progredimos para o Norte do perfil?
# 
#