{
"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
}