{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Course website**: http://www.leouieda.com/geofisica1\n",
    "\n",
    "**Note**: This notebook is part of the course \"Geofísica 1\" of Geology program of the \n",
    "[Universidade do Estado do Rio de Janeiro](http://www.uerj.br/). \n",
    "All content can be freely used and adapted under the terms of the \n",
    "[Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).\n",
    "\n",
    "![Creative Commons License](https://i.creativecommons.org/l/by/4.0/88x31.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prática 8 - Magnetometria - Transformações da anomalia\n",
    "\n",
    "## Objetivos\n",
    "\n",
    "* Aplicar a Transformada de Fourier para calcular transformações da anomalia magnética de campo total: amplitude da derivada total, continuação para cima e redução ao polo.\n",
    "* Verificar como o comportamento dessas transformações é influenciado pela inclinação e declinação magnética.\n",
    "* Observar os efeitos do erro aleatório nas derivadas calculadas.\n",
    "* Aprender as utilidades e limites das transformações."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instruções\n",
    "\n",
    "O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.\n",
    "\n",
    "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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import seaborn\n",
    "import fatiando\n",
    "from fatiando import mesher, utils, gridder\n",
    "from fatiando.gravmag import prism, transform\n",
    "from fatiando.vis import myv, mpl\n",
    "seaborn.set_context('talk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Amplitude da Derivada Total (Total Gradient Amplitude - TGA)\n",
    "\n",
    "A Transformada de Fourier nos permite calcular derivadas dos nossos dados. Uma transformação muito utilizada na magnetometria é a Amplitude da Derivada Total, em inglês Total Gradient Amplitude (TGA):\n",
    "\n",
    "$$\n",
    "TGA = \\sqrt{\\left(\\frac{\\partial \\Delta T}{\\partial x}\\right)^2 + \\left(\\frac{\\partial\\Delta  T}{\\partial y}\\right)^2 + \\left(\\frac{\\partial\\Delta  T}{\\partial z}\\right)^2}\n",
    "$$\n",
    "\n",
    "**O TGA é famoso por concentrar a anomalia magnética de campo total sobre o corpo causador da anomalia.** \n",
    "\n",
    "As derivadas parciais na equação acima podem ser calculadas da anomalia de campo total ($\\Delta T$) usando a Transformada de Fourier, como vimos na aula passada. Para lembrar, podemos calcular a derivada de uma função $h(t)$ através de sua transformada $H(f)$\n",
    "\n",
    "$$\n",
    "W(f) = H(f) i 2 \\pi f\n",
    "$$\n",
    "\n",
    "em que $W(f)$ é a Transformada de Fourier de $\\partial h/ \\partial t$. Uma conhecendo a transformada, podemos calcular a derivada através da transformada inversa de Fourier\n",
    "\n",
    "$$\n",
    "\\frac{\\partial h}{\\partial t} = \\int\\limits_{-\\infty}^{\\infty} W(f) e^{i 2 \\pi f t} df\n",
    "$$\n",
    "\n",
    "### Mitos \n",
    "\n",
    "A Amplitude da Derivada Total é erroneamente chamada de Sinal Analítico 3D (3D Analytic Signal Amplitude - ASA). **Esse nome é errado** (Reid, 2012). Na verdade o que temos é o tamanho (amplitude) do gradiente da anomalia (vetor das derivadas em cada direção).\n",
    "\n",
    "Um mito persistente é o TGA não depende da direção do campo da Terra ou da magnetização do corpo (inclinação e declinação). Diversos autores já mostraram que esse **não é o caso** (ver Li, 2006). \n",
    "\n",
    "\n",
    "#### Referencias\n",
    "\n",
    "[Li, X. (2006), Understanding 3D analytic signal amplitude, GEOPHYSICS, 71(2), L13–L16, doi:10.1190/1.2184367.](http://library.seg.org/doi/abs/10.1190/1.2184367)\n",
    "\n",
    "[Reid, A. (2012), Forgotten truths, myths and sacred cows of Potential Fields Geophysics - II, in SEG Technical Program Expanded Abstracts 2012, pp. 1–3, Society of Exploration Geophysicists.](http://library.seg.org/doi/abs/10.1190/segam2012-0178.1)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulação\n",
    "\n",
    "Vamos explorar a cara do TGA causado por um modelo simples de 1 prisma alongado. Vocês poderão controlar a inclinação e declinação do campo magnético da Terra e o erro aleatório aplicado ao dado. Assumiremos somente magnetização induzida.\n",
    "\n",
    "A célula abaixo cria nosso modelo de 1 prisma e gera uma figura 3D. Feche a figura antes de continuar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shape = (100, 100)\n",
    "x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)\n",
    "dx, dy = 500, 5000\n",
    "cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)\n",
    "cubo_area = cubo.get_bounds()[:4][::-1]\n",
    "bounds = (-5000, 5000, -5000, 5000, 0, 5000)\n",
    "\n",
    "scene = myv.figure()\n",
    "myv.prisms([cubo])\n",
    "myv.outline(bounds)\n",
    "myv.wall_bottom(bounds)\n",
    "myv.wall_west(bounds)\n",
    "oa = myv.mlab.orientation_axes()\n",
    "oa.axes.x_axis_label_text = 'N'\n",
    "oa.axes.y_axis_label_text = 'E'\n",
    "oa.text_property.color = (0.0, 0.0, 0.0)\n",
    "myv.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def derivada_total(inc, dec, erro):\n",
    "    tf = prism.tf(x, y, z, [cubo], inc, dec, pmag=utils.ang2vec(1, inc, dec))\n",
    "    if erro > 0:\n",
    "        tf = utils.contaminate(tf, erro, seed=0)\n",
    "    total = transform.tga(x, y, tf, shape)\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(14, 9))\n",
    "    for ax, data, title in zip(axes.ravel(), [tf, total], ['Anomalia', 'Derivada Total']):\n",
    "        ax.set_aspect('equal')\n",
    "        plt.sca(ax)\n",
    "        plt.title(title)\n",
    "        mpl.square(cubo_area, style='--w')\n",
    "        scale = np.abs([data.min(), data.max()]).max()\n",
    "        plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape), \n",
    "                     30, cmap='RdBu_r', vmin=-scale, vmax=scale)\n",
    "        cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9) \n",
    "        cb.set_label('nT/m' if title != \"Anomalia\" else \"nT\")  \n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "    plt.tight_layout()    \n",
    "widgets.interactive(derivada_total, \n",
    "                    inc=widgets.IntSlider(min=-90, max=90, step=5, value=45),\n",
    "                    dec=widgets.IntSlider(min=-90, max=90, step=5, value=0),\n",
    "                    erro=widgets.FloatSlider(min=0, max=20, step=1, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Redução ao pólo\n",
    "\n",
    "A anomalia magnética é muito complicada de interpretar. Por isso existem tantos filtros e transformações. A melhor situação é quando a inclinação do campo magnético é $\\pm90^\\circ$ (nos polos). Nessas ocasiões, a anomalia se concentra sobre o corpo causador.\n",
    "\n",
    "Felizmente, há um jeito de calcular como a anomalia que medimos ficaria se estivesse nos polos. A técnica chama-se (numa explosão de criatividade) \"Redução ao Polo\". Um dos jeitos de se calcular a redução ao polo é usando a transformada de Fourier.\n",
    "\n",
    "**É necessário conhecer a direção de magnetização do corpo para aplicar a redução ao polo**. Isso é fácil se o corpo tiver somente magnetização induzida pelo campo geomagnético. A magnetização é paralela ao campo da Terra. A situação complica quando há magnetização remanente, aquela que os minerais ferromagnéticos guardam quando se resfriam abaixo da [temperatura de Curie](http://en.wikipedia.org/wiki/Curie_temperature). Por isso a redução ao polo não é tão facilmente utilizada quanto o TGA.\n",
    "\n",
    "Vamos utilizar a mesma simulação que fizemos acima para testar a redução ao polo. Você controlará a inclinação e declinação do campo da Terra (`inc` e `dec`). Vamos assumir **somente magnetização induzida**, logo a inclinação e declinação que será utilizada na redução ao polo é a **mesma que a da Terra**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "shape = (100, 100)\n",
    "x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)\n",
    "dx, dy = 500, 5000\n",
    "cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)\n",
    "cubo_area = cubo.get_bounds()[:4][::-1]\n",
    "def reducao_polo(erro, inc, dec):\n",
    "    tf = prism.tf(x, y, z, [cubo], inc, dec, pmag=utils.ang2vec(5, inc, dec))\n",
    "    if erro > 0:\n",
    "        tf = utils.contaminate(tf, erro, seed=0)\n",
    "    rtp = transform.reduce_to_pole(x, y, tf, shape, inc, dec, inc, dec)\n",
    "    fig, axes = mpl.subplots(1, 2, figsize=(14, 9))\n",
    "    for ax, data, title in zip(axes.ravel(), [tf, rtp], ['Anomalia', 'Reduzida ao polo']):\n",
    "        ax.set_aspect('equal')\n",
    "        plt.sca(ax)\n",
    "        plt.title(title)\n",
    "        mpl.square(cubo_area, style='--w')\n",
    "        scale = np.abs([data.min(), data.max()]).max()\n",
    "        plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape), \n",
    "                     30, cmap='RdBu_r', vmin=-scale, vmax=scale)\n",
    "        cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9) \n",
    "        cb.set_label('nT')  \n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "    plt.tight_layout()    \n",
    "widgets.interactive(reducao_polo,\n",
    "                    erro=widgets.FloatSlider(min=0, max=20, step=1, value=0), \n",
    "                    inc=widgets.IntSlider(min=-90, max=90, step=5, value=45),\n",
    "                    dec=widgets.IntSlider(min=-90, max=90, step=5, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuação para cima\n",
    "\n",
    "Em gravimetria e magnetometria, uma transformação que pode ser feita com os dados é a \"continuação para cima\". Essa transformação nos permite calcular como a anomalia medida seria se estivesse a uma altitude maior. Um dos jeitos de calcular a continuação é através da transformada de Fourier:\n",
    "\n",
    "$$\n",
    "H_{up} = H e^{-\\Delta z |k|}\n",
    "$$\n",
    "\n",
    "em que $H$ é a transformada dos dados a uma altitude $z$, $H_{up}$ é a transformada dos dados a uma altitude $z + \\Delta z$ e $|k|$ são as frequências."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos testar a continuação para cima em uma anomalia causada por dois corpos: um grande e profundo e outro razo e pequeno."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shape = (100, 100)\n",
    "x, y, z = gridder.regular((-6000, 6000, -6000, 6000), shape, z=0)\n",
    "modelo = [mesher.Prism(-200, 200, -200, 200, 100, 500),\n",
    "          mesher.Prism(-4000, 4000, -4000, 4000, 1000, 5000)]\n",
    "inc, dec = 45, -10\n",
    "mag = utils.ang2vec(1, inc, dec)\n",
    "modelo[0].props['magnetization'] = 1*mag\n",
    "modelo[1].props['magnetization'] = 1*mag\n",
    "bounds = (-6000, 6000, -6000, 6000, 0, 6000)\n",
    "\n",
    "myv.figure()\n",
    "myv.prisms(modelo)\n",
    "myv.outline(bounds)\n",
    "myv.wall_bottom(bounds)\n",
    "myv.wall_west(bounds)\n",
    "oa = myv.mlab.orientation_axes()\n",
    "oa.axes.x_axis_label_text = 'N'\n",
    "oa.axes.y_axis_label_text = 'E'\n",
    "oa.text_property.color = (0.0, 0.0, 0.0)\n",
    "myv.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def continuacao(altitude, erro):\n",
    "    cubo_area = modelo[0].get_bounds()[:4][::-1]\n",
    "    tf = prism.tf(x, y, z, modelo, inc, dec)\n",
    "    if erro > 0:\n",
    "        tf = utils.contaminate(tf, erro, seed=0)\n",
    "    cont = transform.upcontinue(x, y, tf, shape, altitude + 1e-5)\n",
    "    fig, axes = mpl.subplots(1, 2, figsize=(14, 9))\n",
    "    titles = ['Anomalia em 0 m', 'Continuada para {:.0f} m'.format(altitude)]\n",
    "    for ax, data, title in zip(axes.ravel(), [tf, cont], titles):\n",
    "        ax.set_aspect('equal')\n",
    "        plt.sca(ax)\n",
    "        plt.title(title)\n",
    "        for prisma in modelo:\n",
    "            mpl.square(prisma.get_bounds()[:4][::-1], style='--w')\n",
    "        scale = np.abs([data.min(), data.max()]).max()\n",
    "        plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape), \n",
    "                     30, cmap='RdBu_r', vmin=-scale, vmax=scale)\n",
    "        cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9) \n",
    "        cb.set_label('nT')  \n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "    plt.tight_layout()    \n",
    "widgets.interactive(continuacao, \n",
    "                    altitude=widgets.FloatSlider(min=0, max=4000, step=50, value=0),\n",
    "                    erro=widgets.FloatSlider(min=0, max=20, step=1, value=0))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}