{ "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 4 - Gravimetria - Estudos de caso da anomalia da gravidade" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Objetivos\n", "\n", "* Observar a topografia, distúrbio e anomalia Bouguer de: Havaí, Baía de Hudson, Japão e costa brasileira.\n", "* Fixar os conceitos de distúrbio da gravidade, anomalia Bouguer e isostasia.\n", "* Fazer gráficos de perfis cortando as principais feições de cada alvo.\n", "* Analisar a relação entre a topografia, isostasia, anomalias e tectônica." ] }, { "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": "markdown", "metadata": {}, "source": [ "## Preparação\n", "\n", "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)." ] }, { "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", "from mpl_toolkits.basemap import Basemap\n", "import ipywidgets as widgets\n", "from IPython.display import display, HTML\n", "import seaborn\n", "import fatiando\n", "from fatiando.gravmag import normal_gravity\n", "from fatiando import gridder\n", "from icgem import load_icgem_gdf" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"Usando a versão do Fatiando: {}\".format(fatiando.__version__))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def data_minmax(data):\n", " \"\"\"Retorna vmin e vmax para que o valor 0 esteja no centro da escala de cor.\"\"\"\n", " ranges = np.abs([data.min(), data.max()]).max()\n", " return dict(vmin=-ranges, vmax=ranges)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Havaí\n", "\n", "Nosso primeiro alvo é o Havaí:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "HTML('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "havai = load_icgem_gdf('../data/eigen-6c3stat-havai.gdf')\n", "havai['topo'] = load_icgem_gdf('../data/etopo1-havai.gdf', usecols=[-1])['topography_shm']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora posso calcular o distúrbio e a anomalia Bouguer." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "havai['disturbio'] = havai['gravity_earth'] - normal_gravity.gamma_closed_form(havai['latitude'], havai['h_over_geoid'])\n", "havai['bouguer'] = havai['disturbio'] - normal_gravity.bouguer_plate(havai['topo'], density_rock=2670, density_water=1040)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "s, n, w, e = havai['area']\n", "havai_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, \n", " lat_ts=0.5*(s + n), resolution='h')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x, y = havai_bm(havai['longitude'], havai['latitude'])\n", "fig, axes = plt.subplots(1, 3, figsize=(14, 9))\n", "ax = axes[0]\n", "tmp = havai_bm.contourf(x, y, havai['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(havai['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "havai_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "ax = axes[1]\n", "tmp = havai_bm.contourf(x, y, havai['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(havai['disturbio']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "havai_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Disturbio')\n", "ax = axes[2]\n", "tmp = havai_bm.contourf(x, y, havai['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(havai['bouguer']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "havai_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Anomalia Bouguer')\n", "plt.tight_layout(pad=0, w_pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora que temos os mapas, vamos ver como é o gráfico da **anomalia Bouguer pela topografia equivalente**.\n", "\n", "Lembrando que a topografia equivalente é:\n", "\n", "$$\n", "h_{eq} = \n", "\\begin{cases}\n", " h & \\quad \\text{if } h \\ge 0 \\\\\n", " h \\frac{\\rho_c - \\rho_a}{\\rho_c} & \\quad \\text{if } h < 0 \\\\\n", "\\end{cases}\n", "$$\n", "\n", "em que $h$ é a topografia, $\\rho_c$ é a densidade da crosta e $\\rho_a$ é a densidade da água do mar." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "heq = havai['topo'].copy()\n", "heq[heq < 0] *= (2670 - 1040)/2670\n", "plt.figure(figsize=(7, 7))\n", "plt.title('Bouguer x Topografia - Havai')\n", "plt.plot(heq, havai['bouguer'], '.k')\n", "plt.xlabel('Topografia equivalente (m)')\n", "plt.ylabel('Anomalia Bouguer (mGal)')\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "É 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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1, p2 = havai_bm(197.5, 15), havai_bm(207, 27.5)\n", "\n", "plt.figure(figsize=(8, 7))\n", "ax = plt.subplot(111)\n", "tmp = havai_bm.contourf(x, y, havai['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(havai['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "havai_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "xp, yp = np.transpose([p1, p2])\n", "havai_bm.plot(xp, yp, 'o-k', linewidth=2)\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para extrair o perfil dos dados." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmp = gridder.profile(x, y, havai['topo'], p1, p2, 100)\n", "perfil_havai = dict(distancia=tmp[2], topo=tmp[3])\n", "perfil_havai['bouguer'] = gridder.profile(x, y, havai['bouguer'], p1, p2, 100)[-1]\n", "perfil_havai['disturbio'] = gridder.profile(x, y, havai['disturbio'], p1, p2, 100)[-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora podemos fazer um gráfico com os 3 dados (topografia e anomalias) ao longo do perfil." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))\n", "ax1, ax2 = axes\n", "d = perfil_havai['distancia']*0.001\n", "ax2.fill_between([d.min(), d.max()], [0, 0], -6000, color='#2780E3')\n", "ax2.fill_between(d, perfil_havai['topo'], -6000, color='#333333')\n", "ax2.set_ylabel('Topografia (m)')\n", "ax2.set_xlabel('Distancia (km)')\n", "\n", "ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')\n", "ax1.plot(d, perfil_havai['bouguer'], '-', label='Bouguer')\n", "ax1.plot(d, perfil_havai['disturbio'], '-', label='Disturbio')\n", "ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')\n", "ax1.legend(loc='center right')\n", "ax1.set_xlim(d.min(), d.max())\n", "plt.tight_layout(h_pad=0, pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* As ilhas estão em equilíbrio isostático?\n", "* Por que o distúrbio é positivo em cima das ilhas?\n", "* Por que o distúrbio tem seu mínimo logo ao lado das ilhas?\n", "* 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?\n", "* Como e por que as feições do distúrbio se correlacionam com a topografia?\n", "* Por que a anomalia Bouguer é sempre positiva?\n", "* Por que a anomalia Bouguer é menos positiva sobre as ilhas?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Japão" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos repetir o mesmo que fizemos antes mas para os dados do Japão:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "HTML('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Primeiro, carregar os dados e calcular o distúrbio e anomalia Bouguer." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "japao = load_icgem_gdf('../data/eigen-6c3stat-japao.gdf')\n", "japao['topo'] = load_icgem_gdf('../data/etopo1-japao.gdf', usecols=[-1])['topography_shm']\n", "japao['disturbio'] = japao['gravity_earth'] - normal_gravity.gamma_closed_form(japao['latitude'], japao['h_over_geoid'])\n", "japao['bouguer'] = japao['disturbio'] - normal_gravity.bouguer_plate(japao['topo'], density_rock=2670, density_water=1040)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos ver como ficam as anomalias e a topografia para essa região." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s, n, w, e = japao['area']\n", "japao_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, \n", " lat_ts=0.5*(s + n), resolution='l')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x, y = japao_bm(japao['longitude'], japao['latitude'])\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 10))\n", "ax = axes[0]\n", "tmp = japao_bm.contourf(x, y, japao['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(japao['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "japao_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "ax = axes[1]\n", "tmp = japao_bm.contourf(x, y, japao['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(japao['disturbio']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "japao_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Disturbio')\n", "ax = axes[2]\n", "tmp = japao_bm.contourf(x, y, japao['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(japao['bouguer']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "japao_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Anomalia Bouguer')\n", "plt.tight_layout(pad=0, w_pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora podemos fazer o gráfico de anomalia Bouguer por topografia equivalente." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "heq = japao['topo'].copy()\n", "heq[heq < 0] *= (2670 - 1040)/2670\n", "plt.figure(figsize=(7, 7))\n", "plt.title('Bouguer x Topografia - Japao')\n", "plt.plot(heq, japao['bouguer'], '.k')\n", "plt.xlabel('Topografia equivalente (m)')\n", "plt.ylabel('Anomalia Bouguer (mGal)')\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos extrair um perfil cortando a zona de subducção, ilhas do Japão, Mar do Japão e chegando no continente." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1, p2 = japao_bm(130, 43), japao_bm(148, 38)\n", "\n", "plt.figure(figsize=(9, 5))\n", "ax = plt.subplot(111)\n", "tmp = japao_bm.contourf(x, y, japao['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(japao['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "japao_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "xp, yp = np.transpose([p1, p2])\n", "japao_bm.plot(xp, yp, 'o-k', linewidth=2)\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmp = gridder.profile(x, y, japao['topo'], p1, p2, 100)\n", "perfil_japao = dict(distancia=tmp[2], topo=tmp[3])\n", "perfil_japao['bouguer'] = gridder.profile(x, y, japao['bouguer'], p1, p2, 100)[-1]\n", "perfil_japao['disturbio'] = gridder.profile(x, y, japao['disturbio'], p1, p2, 100)[-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))\n", "ax1, ax2 = axes\n", "d = perfil_japao['distancia']*0.001\n", "ax2.fill_between([d.min(), d.max()], [0, 0], -10000, color='#2780E3')\n", "ax2.fill_between(d, perfil_japao['topo'], -10000, color='#333333')\n", "ax2.set_ylabel('Topografia (m)')\n", "ax2.set_xlabel('Distancia (km)')\n", "\n", "ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')\n", "ax1.plot(d, perfil_japao['bouguer'], '-', label='Bouguer')\n", "ax1.plot(d, perfil_japao['disturbio'], '-', label='Disturbio')\n", "ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')\n", "ax1.legend(loc='center right')\n", "ax1.set_xlim(d.min(), d.max())\n", "plt.tight_layout(h_pad=0, pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* Por que o distúrbio é praticamente zero na placa do Pacífico e no Mar do Japão?\n", "* Por que o distúrbio apresenta um par negativo-positivo quando passamos da crosta oceânica para as ilhas do Japão?\n", "* Qual é a relação entre o distúrbio, a topografia e a subducção?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Costa brasileira\n", "\n", "Vamos carregar os dados e fazer a análise para a margem leste brasileira (um lugar que esperamos estar em equilíbrio isostático)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "HTML('')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "brasil = load_icgem_gdf('../data/eigen-6c3stat-brasil.gdf')\n", "brasil['topo'] = load_icgem_gdf('../data/etopo1-brasil.gdf', usecols=[-1])['topography_shm']\n", "brasil['disturbio'] = brasil['gravity_earth'] - normal_gravity.gamma_closed_form(brasil['latitude'], brasil['h_over_geoid'])\n", "brasil['bouguer'] = brasil['disturbio'] - normal_gravity.bouguer_plate(brasil['topo'], density_rock=2670, density_water=1040)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s, n, w, e = brasil['area']\n", "brasil_bm = Basemap(projection='merc', llcrnrlat=s, urcrnrlat=n, llcrnrlon=w, urcrnrlon=e, \n", " lat_ts=0.5*(s + n), resolution='l')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x, y = brasil_bm(brasil['longitude'], brasil['latitude'])\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 10))\n", "ax = axes[0]\n", "tmp = brasil_bm.contourf(x, y, brasil['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(brasil['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Topografia')\n", "ax = axes[1]\n", "tmp = brasil_bm.contourf(x, y, brasil['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(brasil['disturbio']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Disturbio')\n", "ax = axes[2]\n", "tmp = brasil_bm.contourf(x, y, brasil['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(brasil['bouguer']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Anomalia Bouguer')\n", "plt.tight_layout(pad=0, w_pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "heq = brasil['topo'].copy()\n", "heq[heq < 0] *= (2670 - 1040)/2670\n", "plt.figure(figsize=(7, 7))\n", "plt.title('Bouguer x Topografia - Brasil')\n", "plt.plot(heq, brasil['bouguer'], '.k')\n", "plt.xlabel('Topografia equivalente (m)')\n", "plt.ylabel('Anomalia Bouguer (mGal)')\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1, p2 = brasil_bm(360 - 58.5, -20.5), brasil_bm(360 - 33.5, -31)\n", "\n", "plt.figure(figsize=(9, 5))\n", "ax = plt.subplot(111)\n", "tmp = brasil_bm.contourf(x, y, brasil['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(brasil['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "xp, yp = np.transpose([p1, p2])\n", "brasil_bm.plot(xp, yp, 'o-k', linewidth=2)\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmp = gridder.profile(x, y, brasil['topo'], p1, p2, 200)\n", "perfil_brasil = dict(distancia=tmp[2], topo=tmp[3])\n", "perfil_brasil['bouguer'] = gridder.profile(x, y, brasil['bouguer'], p1, p2, 200)[-1]\n", "perfil_brasil['disturbio'] = gridder.profile(x, y, brasil['disturbio'], p1, p2, 200)[-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))\n", "ax1, ax2 = axes\n", "d = perfil_brasil['distancia']*0.001\n", "ax2.fill_between([d.min(), d.max()], [0, 0], -6000, color='#2780E3')\n", "ax2.fill_between(d, perfil_brasil['topo'], -6000, color='#333333')\n", "ax2.set_ylabel('Topografia (m)')\n", "ax2.set_xlabel('Distancia (km)')\n", "\n", "ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')\n", "ax1.plot(d, perfil_brasil['bouguer'], '-', label='Bouguer')\n", "ax1.plot(d, perfil_brasil['disturbio'], '-', label='Disturbio')\n", "ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')\n", "ax1.legend(loc='center right')\n", "ax1.set_xlim(d.min(), d.max())\n", "plt.tight_layout(h_pad=0, pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "heq = brasil['topo'].copy()\n", "heq[heq < 0] *= (2670 - 1040)/2670\n", "A = np.ones((heq.size, 2))\n", "A[:, 0] = heq\n", "p = np.linalg.solve(A.T.dot(A), A.T.dot(brasil['bouguer']))\n", "efeito_moho = A.dot(p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(7, 7))\n", "plt.title('Bouguer x Topografia - Brasil')\n", "plt.plot(heq, brasil['bouguer'], '.k', label='Dados observados')\n", "plt.plot(heq, efeito_moho, '-', label='y = {:.3f} heq + {:.3f}'.format(*p))\n", "plt.legend(loc='upper right')\n", "plt.xlabel('Topografia equivalente (m)')\n", "plt.ylabel('Anomalia Bouguer (mGal)')\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora que determinamos o efeito da Moho isostática, podemos removê-lo da nossa anomalia Bouguer." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x, y = brasil_bm(brasil['longitude'], brasil['latitude'])\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 11))\n", "ax = axes[0]\n", "tmp = brasil_bm.contourf(x, y, brasil['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(brasil['disturbio']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Disturbio')\n", "ax = axes[1]\n", "tmp = brasil_bm.contourf(x, y, brasil['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(brasil['bouguer']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Anomalia Bouguer')\n", "ax = axes[2]\n", "tmp = brasil_bm.contourf(x, y, brasil['bouguer'] - efeito_moho, 80, tri=True, ax=ax, cmap='RdBu_r', \n", " vmin=-110, vmax=110)\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "brasil_bm.drawcoastlines(ax=ax)\n", "brasil_bm.drawcountries(ax=ax)\n", "ax.set_title('Anomalia Bouguer - Moho')\n", "plt.tight_layout(pad=0, w_pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* A margem brasileira está em equilíbrio isostático?\n", "* O que estamos presumindo quando removemos o efeito da Moho utilizando o método acima?\n", "* O que sobra na anomalia Bouguer após removermos o efeito da Moho isostática?\n", "* Qual é origem da anomalia residual negativa que vemos na porção continental?\n", "* Qual é origem das anomalias residuais negativas que vemos na porção oceânica?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Hudson Bay" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos olhar como ficam as anomalias na região da Baía de Hudson, no Norte do Canadá:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "HTML('')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hudson = load_icgem_gdf('../data/eigen-6c3stat-hudson.gdf')\n", "hudson['topo'] = load_icgem_gdf('../data/etopo1-hudson.gdf', usecols=[-1])['topography_shm']\n", "hudson['disturbio'] = hudson['gravity_earth'] - normal_gravity.gamma_closed_form(hudson['latitude'], hudson['h_over_geoid'])\n", "hudson['bouguer'] = hudson['disturbio'] - normal_gravity.bouguer_plate(hudson['topo'], density_rock=2670, density_water=1040)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dessa vez, vamos usar uma [projeção polar esterográfica](http://en.wikipedia.org/wiki/Stereographic_projection) para fazer os mapas. \n", "A projeção Mercator distorce demais as regiões de altas latitudes.\n", "Nesses casos, as [projeções cônicas](http://en.wikipedia.org/wiki/Map_projection#Conic) são melhores." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s, n, w, e = hudson['area']\n", "hudson_bm = Basemap(projection='laea', width=4000000, height=3500000,\n", " lat_ts=0.5*(s + n), lat_0=0.5*(s + n), lon_0=0.5*(e + w), resolution='l')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x, y = hudson_bm(hudson['longitude'], hudson['latitude'])\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 12))\n", "ax = axes[0]\n", "tmp = hudson_bm.contourf(x, y, hudson['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(hudson['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "hudson_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "ax = axes[1]\n", "tmp = hudson_bm.contourf(x, y, hudson['disturbio'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(hudson['disturbio']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "hudson_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Disturbio')\n", "ax = axes[2]\n", "tmp = hudson_bm.contourf(x, y, hudson['bouguer'], 80, tri=True, ax=ax, cmap='RdBu_r', \n", " **data_minmax(hudson['bouguer']))\n", "plt.colorbar(tmp, ax=ax, orientation='horizontal', pad=0.01, aspect=30, shrink=1).set_label('mGal')\n", "hudson_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Anomalia Bouguer')\n", "plt.tight_layout(pad=0, w_pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "heq = hudson['topo'].copy()\n", "heq[heq < 0] *= (2670 - 1040)/2670\n", "plt.figure(figsize=(7, 7))\n", "plt.title('Bouguer x Topografia - Bahia de Hudson')\n", "plt.plot(heq, hudson['bouguer'], '.k')\n", "plt.xlabel('Topografia equivalente (m)')\n", "plt.ylabel('Anomalia Bouguer (mGal)')\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "O perfil que vamos extrair vai do continente ao mar do Norte, cortando a baía." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p1, p2 = hudson_bm(260.5, 45.5), hudson_bm(299, 70)\n", "\n", "plt.figure(figsize=(9, 5))\n", "ax = plt.subplot(111)\n", "tmp = hudson_bm.contourf(x, y, hudson['topo'], 80, tri=True, ax=ax, cmap='BrBG_r', \n", " **data_minmax(hudson['topo']))\n", "plt.colorbar(tmp, ax=ax, orientation='vertical', pad=0.01, aspect=30, shrink=1).set_label('m')\n", "hudson_bm.drawcoastlines(ax=ax)\n", "ax.set_title('Topografia')\n", "xp, yp = np.transpose([p1, p2])\n", "hudson_bm.plot(xp, yp, 'o-k', linewidth=2)\n", "plt.tight_layout(pad=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmp = gridder.profile(x, y, hudson['topo'], p1, p2, 200)\n", "perfil_hudson = dict(distancia=tmp[2], topo=tmp[3])\n", "perfil_hudson['bouguer'] = gridder.profile(x, y, hudson['bouguer'], p1, p2, 200)[-1]\n", "perfil_hudson['disturbio'] = gridder.profile(x, y, hudson['disturbio'], p1, p2, 200)[-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 1, sharex=True, figsize=(14, 6))\n", "ax1, ax2 = axes\n", "d = perfil_hudson['distancia']*0.001\n", "ax2.fill_between([d.min(), d.max()], [0, 0], -2000, color='#2780E3')\n", "ax2.fill_between(d, perfil_hudson['topo'], -2000, color='#333333')\n", "ax2.set_ylabel('Topografia (m)')\n", "ax2.set_xlabel('Distancia (km)')\n", "\n", "ax1.set_ylabel('Disturbio e Anomalia Bouguer (mGal)')\n", "ax1.plot(d, perfil_hudson['bouguer'], '-', label='Bouguer')\n", "ax1.plot(d, perfil_hudson['disturbio'], '-', label='Disturbio')\n", "ax1.hlines(0, d.min(), d.max(), linestyles='--', color='#333333')\n", "ax1.legend(loc='center right')\n", "ax1.set_xlim(d.min(), d.max())\n", "plt.tight_layout(h_pad=0, pad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* Esta região está em equilíbrio isostático?\n", "* Por que o distúrbio fica cada vez mais negativo quando mais progredimos para o Norte do perfil?\n", "\n" ] } ], "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 }