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.
Objetivos:
Esse documento que você está usando é um IPython notebook. É um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (que pode ser números, texto, figuras, videos, etc). Esta prática usará a biblioteca Fatiando a Terra de modelagem geofísica e também o módulo numpy.
O notebook é divido em células (como esta). Para editar o conteúdo de uma célula, clique nela (clique nesta para editar esse texto). Para executar uma célula, aperte Shift + Enter
. Execute as duas células abaixo.
%matplotlib inline
from __future__ import division
from IPython.html import widgets
import numpy as np
from fatiando.gravmag import prism, fourier
from fatiando import utils, gridder, mesher
from fatiando.vis import mpl
import fatiando
mpl.rc('lines', linewidth=2)
mpl.rc('font', size=12)
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
areacubo = [-400, 400, 400, 1200]
cubo = mesher.Prism(areacubo[0], areacubo[1], -50000, 50000, areacubo[2], areacubo[3])
area = (-1200, 1200, -400, 2000)
x, y, z = gridder.regular(area, (13, 13), z=0)
def total_cubo(remanente, inc):
def mask(area, v):
x1, x2, y1, y2 = area
v[(x >= x1) & (x <= x2) & (y >= y1) & (y <= y2)] = 0
return v
terra = utils.ang2vec(5000, 45, 0)
remanente = utils.ang2vec(remanente, inc, 0)
induzido = utils.ang2vec(20, 45, 0)
total = induzido + remanente
cubo.addprop('magnetization', total)
bxi = mask(areacubo, prism.bx(x, z, y, [cubo], pmag=induzido))
bzi = mask(areacubo, prism.bz(x, z, y, [cubo], pmag=induzido))
bxr = mask(areacubo, prism.bx(x, z, y, [cubo], pmag=remanente))
bzr = mask(areacubo, prism.bz(x, z, y, [cubo], pmag=remanente))
xp = np.linspace(area[0], area[1], 100)
zp = np.zeros_like(xp)
yp = zp
tf = prism.tf(xp, yp, zp, [cubo], 45, 0)
tx = bxi + bxr + terra[0]
tz = bzi + bzr + terra[2]
fig, axes = mpl.subplots(2, 2, sharex='col', figsize=(10, 8))
ax1, ax2, ax3, ax4 = axes.ravel()
for ax in [ax1, ax2, ax3]:
ax.set_aspect('equal')
mpl.sca(ax)
mpl.set_area(area)
ax.invert_yaxis()
fig.text(0.05, 0.96, '(a)', fontsize=16)
fig.text(0.52, 0.96, '(b)', fontsize=16)
fig.text(0.05, 0.51, '(c)', fontsize=16)
fig.text(0.52, 0.51, '(d)', fontsize=16)
mpl.sca(ax1)
mpl.square(areacubo, fill='grey')
mpl.quiver(x, y, bxi, bzi, linewidth=0.5,
scale=50, pivot='middle', angles='xy', scale_units='xy')
mpl.hlines(0, area[0], area[1], linewidth=3)
mpl.arrow(0, 800, 5*induzido[0], 5*induzido[2], width=2, color='k')
mpl.ylabel('z (m)')
mpl.title('Campo induzido')
mpl.sca(ax2)
mpl.square(areacubo, fill='grey')
mpl.quiver(x, y, bxr, bzr, linewidth=0.5,
scale=50, pivot='middle', angles='xy', scale_units='xy',
color='r')
mpl.hlines(0, area[0], area[1], linewidth=3)
mpl.arrow(0, 800, 5*remanente[0], 5*remanente[2], width=2, color='r')
mpl.title('Campo remanente')
mpl.sca(ax3)
mpl.square(areacubo, fill='grey')
mpl.quiver(x, y, tx, tz, linewidth=0.5, color='b',
scale=50, pivot='middle', angles='xy', scale_units='xy')
mpl.hlines(0, area[0], area[1], linewidth=3)
mpl.arrow(0, 800, 5*total[0], 5*total[2], width=2, color='b')
mpl.xlabel('x (m)')
mpl.ylabel('z (m)')
mpl.title('Campo total = corpo + Terra')
mpl.sca(ax4)
mpl.plot(xp, tf, '-k')
mpl.title('Anomalia de campo total')
mpl.grid(True)
mpl.hlines(0, area[0], area[1], linewidth=1, linestyles='dashed')
mpl.xlabel('x (m)')
mpl.tight_layout()
widgets.interactive(total_cubo,
remanente=widgets.IntSliderWidget(min=0, max=50, step=5, value=0),
inc=widgets.IntSliderWidget(min=-180, max=180, step=5, value=45))
remanente
controla a intensidade da magnetização remanente.inc
controla a inclinação da magnetização remanente. Note que ela pode variar de -180 a 180.shape = (100, 100)
x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)
dx, dy = 500, 5000
cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)
cubo_area = cubo.get_bounds()[:4][::-1]
def mapa_cubo(remanente, inc, dec):
induzida = utils.ang2vec(10, 45, 0)
remanente = utils.ang2vec(remanente, inc, dec)
tf = prism.tf(x, y, z, [cubo], 45, 0, pmag=induzida + remanente)
dt = fourier.ansig(x, y, tf, shape)
fig, axes = mpl.subplots(1, 2, figsize=(14, 6))
titles = ['Induzida: inc={} dec={} \n Remanente: int={} inc={} dec={}'.format(45, 0, np.linalg.norm(remanente), inc, dec),
'Derivada total']
for data, ax, title in zip([tf, dt], axes, titles):
ax.set_aspect('equal')
mpl.sca(ax)
mpl.title(title)
#mpl.square(cubo_area)
scale = np.abs([data.min(), data.max()]).max()
mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)
mpl.colorbar(pad=0).set_label('nT')
mpl.m2km()
mpl.xlabel('y (km)')
mpl.ylabel('x (km)')
mpl.tight_layout(h_pad=0, w_pad=0)
widgets.interactive(mapa_cubo,
remanente=widgets.FloatSliderWidget(min=0, max=50, step=1, value=0),
inc=widgets.IntSliderWidget(min=-180, max=180, step=5, value=45),
dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=0))
remanente
controla a intensidade da magnetização remanente.inc
e dec
controlam, respectivamente, a inclinação e declinação da magnetização remanente.shape = (100, 100)
x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)
dx, dy = 500, 5000
cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)
cubo_area = cubo.get_bounds()[:4][::-1]
def reducao_polo(inc, dec, remanente):
induzida = utils.ang2vec(5, 45, 0)
remanente = utils.ang2vec(remanente, -45, 35)
total = induzida + remanente
tf = prism.tf(x, y, z, [cubo], 45, 0, pmag=total)
rtp = fourier.reduce_to_pole(x, y, tf, shape, 45, 0, inc, dec)
positivo = 100*np.sum(rtp >= 0)/rtp.size
titles = ['Induzida: inc={} dec={} \n Remanente: int={}'.format(45, 0, np.linalg.norm(remanente)),
u'Reduzida ao polo com inc={} dec={} \n Proporção positiva/negativa = {}%'.format(inc, dec, positivo)]
fig, axes = mpl.subplots(1, 2, figsize=(14, 6))
for ax, data, title in zip(axes.ravel(), [tf, rtp], titles):
ax.set_aspect('equal')
mpl.sca(ax)
mpl.title(title)
#mpl.square(cubo_area)
scale = np.abs([data.min(), data.max()]).max()
mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)
mpl.colorbar(pad=0).set_label('nT')
mpl.m2km()
mpl.tight_layout(h_pad=0, w_pad=0)
return utils.vec2ang(total)
w = widgets.interactive(
reducao_polo,
remanente=widgets.FloatSliderWidget(min=0, max=50, step=5, value=0),
inc=widgets.IntSliderWidget(min=-90, max=90, step=5, value=20),
dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=-30))
w
remanente
controla a intensidade da magnetização remanente.inc
e dec
controlam, respectivamente, a inclinação e declinação utilizadas na redução ao polo.Rode a célula abaixo para conferir se o valor de magnetização que você determinou está correto.
print(u'Magnetização total verdadeira: inc={:.2f} dec={:.2f}'.format(w.result[1], w.result[2]))