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. All content can be freely used and adapted under the terms of the Creative Commons Attribution 4.0 International License.
Esse documento que você está usando é um Jupyter notebook. É 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).
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. 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.
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.
%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
print("Usando a versão do Fatiando: {}".format(fatiando.__version__))
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)
Nosso primeiro alvo é o Havaí:
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
.
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.
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 e escalas de cor apropriadas.
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')
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.
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.
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.
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.
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)
Vamos repetir o mesmo que fizemos antes mas para os dados do Japão:
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.
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.
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')
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.
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.
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)
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]
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)
Vamos carregar os dados e fazer a análise para a margem leste brasileira (um lugar que esperamos estar em equilíbrio isostático).
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>')
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)
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')
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)
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)
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)
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]
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.
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)
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.
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)
Vamos olhar como ficam as anomalias na região da Baía de Hudson, no Norte do Canadá:
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>')
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 para fazer os mapas. A projeção Mercator distorce demais as regiões de altas latitudes. Nesses casos, as projeções cônicas são melhores.
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')
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)
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.
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)
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]
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)