fatiando-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4f46b488-41bd-489a-ba32-8d47cdc22cb7",
   "metadata": {},
   "source": [
    "# CoCo: Learning how to use Fatiando a Terra tools\n",
    "\n",
    "_by [Santiago Soler](https://github.com/santisoler/), [Agustina Pesce](https://github.com/aguspesce/) and [Andrea Balza](https://github.com/andieie)_\n",
    "\n",
    "In this notebook we will learn how we can use the Fatiando a Terra tools to\n",
    "process gravity data (or any spatial data).\n",
    "\n",
    "Processing gravity data might seem easy, although there are a lot of details\n",
    "that we must to take into account.\n",
    "Nevertheless, the idea of this notebook is not to become experts on gravity\n",
    "modelling, but to use it as an excuse for learning more Python and how should\n",
    "we approach learning a new library.\n",
    "\n",
    "A very nice place to start learning a new tool is its website and its\n",
    "documentation pages.\n",
    "\n",
    "The Fatiando a Terra project has its own website: https://www.fatiando.org\n",
    "\n",
    "The project is composed by different open-source Python libraries, each one\n",
    "with a specific goal and set of problems that they tackle.\n",
    "Each one of these libraries has its own documentation pages:\n",
    "- Pooch: https://www.fatiando.org/pooch\n",
    "- Boule: https://www.fatiando.org/boule\n",
    "- Verde: https://www.fatiando.org/verde\n",
    "- Harmonica: https://www.fatiando.org/harmonica\n",
    "\n",
    "## Goals\n",
    "\n",
    "1. Fetch and cache the data using [Pooch](https://github.com/fatiando/pooch)\n",
    "1. Calculate the gravity disturbance using [Boule](https://github.com/fatiando/boule)\n",
    "1. Interpolate data and apply map projections to grids using [Verde](https://github.com/fatiando/verde)\n",
    "1. Perform topographic correction and equivalent-source interpolation using [Harmonica](https://github.com/fatiando/harmonica)\n",
    "---\n",
    "En este notebook vamos a aprender cómo podemos utilizar las herramientas de Fatiando a Terra para procesar datos de gravedad (o cualquier tipo de dato espacial).\n",
    "\n",
    "El procesamiento de datos gravimétricos puede parecer sencillo, sin embargo hay muchos detalles que se deben tener en cuenta. Sin embargo, la idea de este notebook no es convertirnos en expertxs en el modelado gravimétrico, sino utilizarlo como excusa para aprender más sobre Python y cómo debemos proceder a la hora de aprender una nueva librería.\n",
    "\n",
    "Un muy buen lugar para comenzar a aprender una nueva librería es su sitio web y sus páginas de documentación.\n",
    "\n",
    "El proyecto Fatiando a Terra posee su propia página web: https://www.fatiando.org\n",
    "\n",
    "El proyecto está compuesto por diferentes librerías open-source escritas en Python, cada una con objetivos específicos y problemas que permiten resolver. Cada una de estas librerías posee sus propias páginas de documentación:\n",
    "\n",
    "- Pooch: https://www.fatiando.org/pooch\n",
    "- Boule: https://www.fatiando.org/boule\n",
    "- Verde: https://www.fatiando.org/verde\n",
    "- Harmonica: https://www.fatiando.org/harmonica\n",
    "    \n",
    "## Objetivos\n",
    "\n",
    "1. Importar librerías nuevas\n",
    "1. Leer documentación (online y docstrings)\n",
    "1. Descargar datos con Pooch\n",
    "1. Manipular datos espaciales con Verde\n",
    "1. Calcular gravedad normal con Boule\n",
    "1. Procesar datos gravimétricos con Harmonica\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f649168-92ef-421a-b8c1-27f2908fc929",
   "metadata": {},
   "source": [
    "## Study area \n",
    "\n",
    "The **Bushveld Igneous Complex** is located in South Africa and is the largest known layered igneous intrusion. It has been tilted and eroded forming the outcrops around what appears to be the edge of a great geological basin: the Transvaal Basin. It is approximately 2 billion years old and is divided into four different limbs: the northern, southern, eastern, and western limbs. The Bushveld Complex comprises the Rustenburg Layered suite, the Lebowa Granites and the Rooiberg Felsics, that are overlain by the Karoo sediments. See [Webb et al. (2004)](https://doi.org/10.2113/107.1-2.207) for an overview and previous interpretations of the Bushveld in depth.\n",
    "\n",
    "<figure>\n",
    "<img src=\"images/bushveld_igneous_complex_geology.jpg\" alt=\"Geologic map of the Bushveld Complex\">\n",
    "<figcaption><em>\n",
    "    Geology and mineral exploration sites of the Bushveld Igneous Complex \n",
    "    (<a href=\"https://en.wikipedia.org/wiki/File:Bushveld_Igneous_Complex.png\">public domain</a>).\n",
    "</em></figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c5b1a97-5d7a-4c51-9a45-0b42002606db",
   "metadata": {},
   "source": [
    "## Import packages | Importar paquetes\n",
    "\n",
    "Import Numpy, Pandas and Matplotlib for managing and plotting our data.\n",
    "\n",
    "Importemos Numpy, Pandas y Matplotlib para gestionar nuestros datos y graficarlos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b627a1a-d116-499b-8aab-97a2687d9ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0e29b43-0edd-48ec-a685-157d9b09ed74",
   "metadata": {},
   "source": [
    "Lets import **Xarray**, a library similar to Pandas, but it allows us to work with multidimensional data besides tabular data.\n",
    "\n",
    "Importemos **Xarray**, una librería similar a Pandas, pero que nos permite trabajar con datos multidimensionales en vez de datos tabulados."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e657b39-4462-4997-8be7-e462676a0a65",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469cbb09-bee8-4aab-831a-fd0085dc0028",
   "metadata": {},
   "source": [
    "Lets import **pyproj**: it will help us to convert geographic coordinates into Cartesian coordinates through map projections.\n",
    "\n",
    "Importemos **pyproj** para poder convertir coordeanadas geográficas en coordenadas Cartesianas mediante proyecciones.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae146cba-7be6-45f2-b1ee-1ce376bb0fc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf6f26dc-b99b-4659-b502-aea7ada1d3a2",
   "metadata": {},
   "source": [
    "Finally, lets import some Fatiando libraries.\n",
    "\n",
    "Finalmente, importemos algunas librerías de Fatiando."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29141fa2-1461-4f1a-a76b-d87dfc699294",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pooch\n",
    "import verde as vd\n",
    "import boule as bl\n",
    "import harmonica as hm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1016f8b1-c789-406f-92a1-e0a1bba3a1e5",
   "metadata": {},
   "source": [
    "## Fetch data | Descarguemos datos\n",
    "\n",
    "In this notebook we will use three different types of data:\n",
    "\n",
    "- **gravity data**: we will use a set of gravity data over South Africa. These data consist of measurements of the gravity acceleration produced by Earth at different points of observation located on South Africas surface topography. Gravity data are in mGal, a unit which is equivalent to 1e5 m/s2. \n",
    "\n",
    "- **topography of South Africa** : the last step of this notebook aims to remove the gravity effect of the topography from our data. In order to do this, we need the DEM (digital elevation model) or a topography grid.\n",
    "\n",
    "These data are available with a Creative Common Attribution license in https://doi.org/10.5281/zenodo.5167357.\n",
    "Let's see how we can use *Pooch* to download this data to our computer to then load it to our notebook. \n",
    "\n",
    "Lets define the Digital Object Identified (DOI) of the Zenodo file where our data are located.\n",
    "\n",
    "---\n",
    "En este notebook vamos a utilizar tres tipos de datos diferentes:\n",
    "\n",
    "- **datos gravimétricos**: vamos a utilizar un conjunto de datos gravimétricos sobre Sudáfrica.Estos datos consisten en mediciones de la aceleración gravitatoria producida por el planeta Tierra sobre diferentes puntos de observación ubicados sobre la superficie de Sudáfrica. Los datos gravimétricos vienen dados en mGal, una unidad que equivale a 1e5 m/s^2.\n",
    "- **topografía de Sudáfrica**: el último paso de este notebook consistirá en remover el efecto gravitatorio de la topografía de Sudáfrica de nuestros datos. Para ello necesitamos un DEM (modelo de elevación digital) o una grilla de topografía.\n",
    "\n",
    "Todos estos datos se encuentran disponibles bajo licencia Creative Commons Attribution en https://doi.org/10.5281/zenodo.5167357. Veamos cómo podemos utilizar Pooch para descargar estos datos a nuestra computadora para luego poder cargarlos en nuestros notebooks.\n",
    "\n",
    "Definimos el Digital Object Identified (DOI) del archivo de Zenodo en el cual se encuentran nuestros datos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd8d3452-b166-4ba8-8e2a-ac23f55b1fac",
   "metadata": {},
   "outputs": [],
   "source": [
    "doi = \"10.5281/zenodo.5167357\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e03c0b-5126-4ea2-8a45-05cdf3bc7a2f",
   "metadata": {},
   "source": [
    "Let's use the function `pooch.retrieve` to download the data file `southern-africa-gravity.csv.xz` and `earth-topography-10arcmin.nc` . Also, we will use the checksum so Pooch could verify if these files are corrupted or not. \n",
    "\n",
    "Utilicemos la función `pooch.retrieve` para descargar los archivos de datos `southern-africa-gravity.csv.xz` y `earth-topography-10arcmin.nc`. Además, vamos a utilizar la suma de comprobación de estos archivos, de forma tal que Pooch pueda verificar que los archivos que se descarguen no estén corruptos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44670ff7-6e40-4df0-8e72-7efc6df2a040",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity = f\"doi:{doi}/southern-africa-gravity.csv.xz\"\n",
    "gravity_hash = \"md5:1dee324a14e647855366d6eb01a1ef35\"\n",
    "\n",
    "path_gravity = pooch.retrieve(gravity, known_hash=gravity_hash, progressbar=True)\n",
    "path_gravity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee72ba7-219e-4c29-9da7-2852a52e1e19",
   "metadata": {},
   "outputs": [],
   "source": [
    "topography = f\"doi:{doi}/earth-topography-10arcmin.nc\"\n",
    "topography_hash = \"md5:c43b61322e03669c4313ba3d9a58028d\"\n",
    "\n",
    "path_topography = pooch.retrieve(\n",
    "    topography, known_hash=topography_hash, progressbar=True\n",
    ")\n",
    "path_topography"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e9477f8-533c-4274-9106-7529a6268fb7",
   "metadata": {},
   "source": [
    "## Load the data | Cargar los datos\n",
    "\n",
    "Let's start with the **gravity data**. Because the gravity data that we downloaded is a CSV file (zipped in an .xz), we will use Pandas to read them.\n",
    "\n",
    "Comencemos por los **datos gravimétricos**. Dado que los datos de gravedad que descargamos consiste en un archivo CSV (comprimido en un .xz), vamos a utilizar Pandas para leerlo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fd52f6d-5ae2-4cfb-8e70-789e19d24b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(path_gravity)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392e07ce-bf2a-4ece-90f9-9056b6706db5",
   "metadata": {},
   "source": [
    "The topography data are in regular grids as netCDF (.nc) files. To load these data we will use Xarray.\n",
    "\n",
    "Los datos de **topografía** vienen dados en grillas regulares en archivos netCDF (.nc). Para cargar estos datos vamos a hacer uso de Xarray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c03a942c-dbd5-430d-97ad-9ff48079d307",
   "metadata": {},
   "outputs": [],
   "source": [
    "topography = xr.load_dataarray(path_topography)\n",
    "topography"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e638ab90-0550-461d-8c07-f0ff9fb797e6",
   "metadata": {},
   "source": [
    "## Plot the data | Grafiquemos los datos\n",
    "\n",
    "All of the data that we downloaded are georeferenced in eliptical coordinates (lat/long). We will use matplotlib to visualize the data that we just downloaded.\n",
    "\n",
    "> **Note**:\n",
    "> To plot georeferenced data, we recommend libraries that are more suitable, like **pyGMT**\n",
    "> In order to simpligy this notebook, we will use matplotlib\n",
    "\n",
    "---\n",
    "Todos los datos que descargamos se encuentran georeferenciados en coordenadas elipsoidales (latitud, longitud). \n",
    "Utilicemos matplotlib para poder visualizar los datos que descargamos.\n",
    "\n",
    "> **Nota**\n",
    "> Para graficar datos georeferenciados es recomendable utilizar librerías más adecuadas, como **pyGMT**.\n",
    "> Con el objetivo de mantener este notebook más simple, utilizaremos matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "280105cb-a324-4668-ab66-d2c570372c67",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 12))\n",
    "plt.scatter(data.longitude, data.latitude, c=data.gravity_mgal, s=1)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00594538-49f6-4baa-8a36-d21f72307727",
   "metadata": {},
   "source": [
    "## Crop the data to the Bushveld | Recortemos los datos sobre Bushveld\n",
    "\n",
    "As we can see, the gravity data are of all of South Africa. However, our region of interest is the Bushveld complex. We will zoom in to our region and add more external data, so we also avoid boundary effects of previous processes. \n",
    "\n",
    "Lets define the Bushveld region in degrees ( west, east, south, north):\n",
    "\n",
    "Como podemos ver, los datos de gravedad pertenecen a toda Sudáfrica, sin embargo nuestra zona de interés es la región de Bushveld. Las grillas de geoide y topografía vamos a recortarlas a regiones un poco más extensas, con el objetivo de evitar efectos de borde en los posteriores procesos.\n",
    "\n",
    "Definamos la región de Bushveld en grados (oeste, este, sur, norte):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e42ce0e7-1702-4d92-b174-33d44591e92a",
   "metadata": {},
   "outputs": [],
   "source": [
    "region = (25, 32, -27, -23)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15ebff27-1e7b-4861-bf7b-592258461186",
   "metadata": {},
   "source": [
    "We define the region and then add padding cells to the grid of interest. We will use `vd.pad_region()` function from Verde:\n",
    "\n",
    "Agrandemos la region para utilizar con las grillas. Para ello haremos uso de la función `vd.pad_region()` de Verde:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5f71fb6-9f49-469a-8459-f5a63932127f",
   "metadata": {},
   "outputs": [],
   "source": [
    "region_pad = vd.pad_region(region, pad=5)\n",
    "region_pad"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "329a2171-5508-4bce-bd69-e60e715ccfde",
   "metadata": {},
   "source": [
    "We can use `vd.inside` function in Verde to obtain the data that are inside of our region and then cut them using Pandas.\n",
    "\n",
    "Podemos utilizar la función `vd.inside` de Verde para obtener qué datos se encuentran dentro de nuestra región y luego recortarlos con Pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dab30e19-7525-4459-9f83-33613f9bd0a9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "inside_points = vd.inside((data.longitude, data.latitude), region)\n",
    "\n",
    "data = data[inside_points]\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f44146be-2c64-45d8-ac2f-c6a0fed110bd",
   "metadata": {},
   "source": [
    "To cut the grid directly we could use the `sel()` method in the `xarray.DataArray`. \n",
    "\n",
    "Para recortar grillas podemos utilizar directamente el método `sel()` de los `xarray.DataArray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56565e88-0c82-4b19-8aec-611cddde9270",
   "metadata": {},
   "outputs": [],
   "source": [
    "topography = topography.sel(\n",
    "    longitude=slice(region_pad[0], region_pad[1]),\n",
    "    latitude=slice(region_pad[2], region_pad[3])\n",
    ")\n",
    "topography"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29c45e84-57e9-4c2f-8632-f995d9e88288",
   "metadata": {},
   "source": [
    "## Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "606381d4-572d-42a1-89fe-adb94fbcf9a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 12))\n",
    "plt.scatter(data.longitude, data.latitude, c=data.gravity_mgal, s=10)\n",
    "plt.colorbar(label=\"mGal\", shrink=0.6)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b703a92b-a91e-4d57-9252-f5c91f0d7ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 12))\n",
    "topography.plot()\n",
    "plt.scatter(data.longitude, data.latitude, c=\"black\", s=0.2)\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8952de6-1417-4bbf-b4e9-e3a8826301e6",
   "metadata": {},
   "source": [
    "## Project the data and the grid | Proyectando los datos y la grilla\n",
    "\n",
    "\n",
    "Because our region of interest is small enough, we can use a projection to transform our data from geographic coordinates to cartesian coordinates. This simplifies our tasks like interpolating and topographic corrections, that tend to be alot faster in castesian coordinates than espheric of geographic coordinates.\n",
    "\n",
    "We will make use of `pyproj` to create a Mercator projection type with latitude in \"real scale\" centered on the data.\n",
    "\n",
    "Dado que nuestra región es los suficientemente pequeña podemos utilizar una proyección para transformar nuestros datos en coordenadas geográficas (longitud, latitud) a coordenadas Cartesianas. Esto simplifica mucho tareas como  la interpolación y las correcciones topográficas, que suelen ser mucho más rápidas en coordenadas Cartesianas que en coordenadas esféricas o geográficas.\n",
    "\n",
    "Haremos uso de `pyproj` para crear una proyección de tipo Mercator con la latitud de la \"escala real\" centrada en los datos.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b0ed76d-082e-4db1-bebb-06d88207536b",
   "metadata": {},
   "outputs": [],
   "source": [
    "projection = pyproj.Proj(proj=\"merc\", lat_ts=data.latitude.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "285e770b-8afd-4ca2-8ab5-123ece556deb",
   "metadata": {},
   "source": [
    "Let's use this `projection` to transform the gravity data coordinates and then include them in the dame `pandas.DataFrame`\n",
    "Utilicemos esta `projection` para transformar las coordenadas de nuestros datos de gravedad y luego incluirlas en el mismo `pandas.DataFrame`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abb95ef3-a760-4a22-918f-8adadba00c32",
   "metadata": {},
   "outputs": [],
   "source": [
    "easting, northing = projection(data.longitude.values, data.latitude.values)\n",
    "data = data.assign(easting_m=easting, northing_m=northing)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5e4ee6d-1965-4e6e-808d-5395980a1ca6",
   "metadata": {},
   "source": [
    "Projecting a grid is more complicated, it involves interpolating with the goal of having equidistance nodes. Verde has a function called `project_grid` which calculates this automatically. We will use a nearest neighbor interpolation because it is one the faster methods and the projections that we choose will not deform our data significantly. \n",
    "\n",
    "Projectar grillas es más complicado, ya que involucra realizar algunas interpolaciones a fines de mantener los nodos equiespaciados. Verde posee la función `project_grid` que nos permite calcular todo esto de forma automática. Utilizaremos una interpolación a primeros vecinos dado que es una de las más rápidas y la proyección que elegimos no deforma tanto nuestros datos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ac9c9d-9f35-497e-98e7-41f69843649a",
   "metadata": {},
   "outputs": [],
   "source": [
    "topography_proj = vd.project_grid(topography, projection, method=\"nearest\")\n",
    "topography_proj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "252e21d3-fbce-44a3-a02a-2c856a7ddf88",
   "metadata": {},
   "source": [
    "## Gravity disturbance | Disturbio de gravedad\n",
    "\n",
    "The gravity disturbance is defined as:\n",
    "\n",
    "$$ \\delta g(\\lambda, \\phi, h) = g(\\lambda, \\phi, h) - \\gamma(\\lambda, \\phi, h) $$\n",
    "\n",
    "where $(\\lambda, \\phi, h)$ are the longitudinal, latitudinal and altitude coordinates of the observation point, and $\\gamma$ is the **normal gravity** generated by the reference ellipsoid. Note that these three quantities are located at the same observation point.\n",
    "\n",
    "> **Note**\n",
    "> The altitudes of the data we download are referred to sea level, that is, to the geoid. However, it is recommended to work with altitudes above the reference ellipsoid.\n",
    "> In order to keep this notebook short, we are going to skip the step of transforming altitudes on geoid to altitudes on ellipsoid.\n",
    "> Since the differences between the two are small, we won't see any noticeable effects of skipping that step.\n",
    "> For more information on how to perform this step, refer to: https://github.com/fatiando/tutorials\n",
    "\n",
    "The values ​​of $g$ on the observation points are given in the data that we downloaded, therefore we have to calculate the normal gravity in the same points.\n",
    "Boule allows us to define a reference ellipsoid (we will use WGS84) and calculate its gravitational effect at any external point.\n",
    "\n",
    "---\n",
    "\n",
    "El disturbio de gravedad se define como:\n",
    "\n",
    "$$ \\delta g(\\lambda, \\phi, h) = g(\\lambda, \\phi, h) - \\gamma(\\lambda, \\phi, h) $$\n",
    "\n",
    "donde $(\\lambda, \\phi, h)$ son las coordenadas longitudinales, latitudinales y la altitud del punto de observación, y $\\gamma$ es la **gravedad normal** generada por el elipsoide de referencia. Nótese que estas tres cantidades se ubican en el mismo punto de observación.\n",
    "\n",
    "> **Nota**\n",
    "> Las altitudes de los datos que descargamos están referidas al nivel del mar, es decir, al geoide. Sin embargo, es recomendable trabajar con altitudes sobre el elipsoide de referencia.\n",
    "> Con el objetivo de mantener este notebook breve, vamos a saltarnos el paso de transformar las altitudes sobre geoide a altitudes sobre elipsoide.\n",
    "> Debido a que las diferencias entre ambas son pequeñas, no veremos efectos apreciables al omitir ese paso.\n",
    "> Para más información sobre cómo se puede realizar este paso referirse a: https://github.com/fatiando/tutorials\n",
    "\n",
    "Los valores de $g$ sobre los puntos de observación vienen dados en los datos que descargamos, por lo tanto nos resta calcular la gravedad normal en los mismos puntos.\n",
    "Boule nos permite definir un elipsoide de referencia (usaremos el WGS84) y calcular su efecto gravitatorio en cualquier punto exterior."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95770f9f-4d6e-4243-aa9c-2741684d419c",
   "metadata": {},
   "source": [
    "Let's define the reference elipsoid\n",
    "\n",
    "Definimos el elipsoide de referencia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d60a20d5-81c6-4e3c-a2ff-c7e85bf456ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipsoid = bl.WGS84\n",
    "ellipsoid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46145f60-a179-49cd-b7eb-fa9ee2d92faf",
   "metadata": {},
   "source": [
    "We will now calculate normal gravity over our observation points\n",
    "\n",
    "Calculamos gravedad normal sobre los puntos de observación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1bcfccd-c8b1-4c91-932a-dbb5877aaad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "normal_gravity = ellipsoid.normal_gravity(data.latitude, data.height_sea_level_m)\n",
    "normal_gravity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01caed11-ebd0-48ea-885a-3a049713ab2b",
   "metadata": {},
   "source": [
    "Calculamos el disturbio de gravedad y lo agregamos a nuestro `DataFrame`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "890b2229-5493-4be7-af7f-f138577ea82a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity_disturbance = data.gravity_mgal - normal_gravity\n",
    "\n",
    "data = data.assign(gravity_disturbance_mgal=gravity_disturbance)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34cd15a1-4cce-42ec-857d-35b12e08a2bb",
   "metadata": {},
   "source": [
    "Lets plot the gravity disturbance. Since we will make many plots throughout this notebook, it is better to create a function for this purpose:\n",
    "\n",
    "Grafiquemos el disturbio de gravedad. Como vamos a realizar múltiples gráficas a lo largo de este notebook, es mejor crear una funcion para ello:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac914bce-5acc-43ab-a9c5-a1cf4e7896b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter_plot_gravity(easting, northing, gravity_data, title=None):\n",
    "    \"\"\"\n",
    "    Generate scatter plots of gravity data\n",
    "    \n",
    "    The colorbar will be centered around cero and using the seismic colormap.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    easting: array\n",
    "        Easting coordinates of the gravity data in meters.\n",
    "    northing : array\n",
    "        Northing coordinates of the gravity data in meters.\n",
    "    gravity_data : array\n",
    "        Values of the gravity data in mgal.\n",
    "    title : str (optional)\n",
    "        Title of the plot.\n",
    "    \"\"\"\n",
    "    maxabs = vd.maxabs(gravity_data)\n",
    "    plt.figure(figsize=(12, 12))\n",
    "    plt.scatter(\n",
    "        easting,\n",
    "        northing,\n",
    "        c=gravity_data,\n",
    "        s=10,\n",
    "        cmap=\"seismic\",\n",
    "        vmin=-maxabs,\n",
    "        vmax=maxabs,\n",
    "    )\n",
    "    plt.xlabel(\"easting [m]\")\n",
    "    plt.ylabel(\"northing [m]\")\n",
    "    if title is not None:\n",
    "        plt.title(title)\n",
    "    plt.colorbar(label=\"mGal\", shrink=0.6)\n",
    "    plt.gca().set_aspect(\"equal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df491559-2489-47aa-97c1-a6c21815f925",
   "metadata": {},
   "source": [
    "Now we can plot the gravity disturbance using the function we wrote.\n",
    "\n",
    "Ahora sí, grafiquemos el disturbio de gravedad usando la función que escribimos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "328c08ea-4ef1-462f-9ae2-ba946bf9e7d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter_plot_gravity(\n",
    "    data.easting_m,\n",
    "    data.northing_m,\n",
    "    data.gravity_disturbance_mgal,\n",
    "    title=\"Gravity disturbance\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99f641d2-d493-4e8a-97f2-46d0f2bab089",
   "metadata": {},
   "source": [
    "## Topographic correction | Corrección topográfica\n",
    "\n",
    "The gravity disturbance contains all the information about the gravitational acceleration generated by anomalous masses. One of the things that contributes the most to the gravity disturbance is the topography: those bodies that are above the reference ellipsoid.\n",
    "The effect of the topography hides the effect of buried bodies that are of greater interest to us (in this case the igneous rocks of Bushveld).\n",
    "Let us then remove the effect of the topography of the gravity disturbance.\n",
    "\n",
    "Harmonica allows us to calculate the gravitational effects of topography by modeling it with simple geometries such as rectangular prisms.\n",
    "\n",
    "The `prism_layer` function creates a prism model from a topographical surface that undulates around a reference level. In our case the reference level will be the ellipsoid (zero altitude) and the surface will be the topography.\n",
    "\n",
    "The hardest part is assigning the density to each of these prisms:\n",
    "\n",
    "* $\\rho = 2670\\ kg/m^3$ for continental regions with positive topography (above the ellipsoid)\n",
    "* $\\rho = 1040 - 2670\\ kg/m^3$ in the oceans (1040 is the approximate density of sea water)\n",
    "\n",
    "\n",
    "<figure>\n",
    "<img src=\"images/topographic-correction.svg\" alt=\"Sketch of the surfaces and masses involved in topographic correction of gravity disturbances.\">\n",
    "<figcaption><em>\n",
    "    Sketch of the surfaces and masses involved in topographic correction of gravity disturbances.\n",
    "    (CC BY).\n",
    "</em></figcaption>\n",
    "</figure>\n",
    "\n",
    "> **Note**\n",
    "> There are some details to consider when making this type of topographic corrections: the sea level does not coincide with the surface of the ellipsoid, and therefore there are bodies of water above and below the ellipsoid.\n",
    "> In the interest of keeping this notebook brief, we have left out these details.\n",
    "> For more information on this, you can refer to: https://github.com/fatiando/tutorials\n",
    "\n",
    "--- \n",
    "\n",
    "El disturbio de gravedad contiene toda la información acerca de la aceleración gravitatoria generada por masas anómalas. Uno de los conjuntos que mayor contribuye al disturbio de gravedad es la topografía: aquellos cuerpos que se encuentran por encima del elipsoide de referencia.\n",
    "El efecto de la topografía nos esconde el efecto de cuerpos enterrados que nos son de mayor de interés (en este caso las rocas ígneas de Bushveld).\n",
    "Quitemos entonces el efecto de la topografía del disturbio de gravedad.\n",
    "\n",
    "Harmonica nos permite calcular los efectos gravitatorios de la topografía modelándola con geometrías sencillas como son los prismas rectangulares.\n",
    "\n",
    "La función `prism_layer` crea un modelo de prismas a partir de una superficie topográfica que ondula alrededor de un nivel de referencia. En nuestro caso el nivel de referencia será el elipsoide (altitud cero) y la superficie será la topografía.\n",
    "\n",
    "La parte más difícil es asignar la densidad a cada uno de estos prismas:\n",
    "\n",
    "* $\\rho = 2670\\ kg/m^3$ para regiones continentales con topografía positiva (por encima del elipsoide)\n",
    "* $\\rho = 1040 - 2670\\ kg/m^3$ en los oceanos (1040 es la densidad aproximada del agua del mar)\n",
    "\n",
    "> **Nota**\n",
    "> Existen algunos detalles a considerar al realizar este tipo de correcciones topográficas: el nivel del mar no coincide con la superficie del elipsoide, y por ende hay masas de agua por encima y por debajo del elipsoide.\n",
    "> Con el objetivo de mantener este notebook breve, hemos obviado estos detalles.\n",
    "> Para más información al respecto puede referirse a: https://github.com/fatiando/tutorials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "268ba0ee-186f-4fb0-9b0f-62b9acee7b1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "density = np.where(topography_proj > 0, 2670, 1040 - 2670)\n",
    "\n",
    "topography_model = hm.prism_layer(\n",
    "    coordinates=(topography_proj.easting, topography_proj.northing),\n",
    "    surface=topography_proj,\n",
    "    reference=0,\n",
    "    properties={\"density\": density},\n",
    ")\n",
    "topography_model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1f47cb4-de82-4ed9-9853-04f914f4ab9b",
   "metadata": {},
   "source": [
    "We now calculate the gravity effect of these prisms over our observation points\n",
    "\n",
    "Calculemos ahora el efecto gravitatorio de estos prismas sobre los puntos de observación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0e204c7-6ffd-402e-a127-61c5d0337fab",
   "metadata": {},
   "outputs": [],
   "source": [
    "coordinates = (data.easting_m, data.northing_m, data.height_sea_level_m)\n",
    "\n",
    "terrain_effect = topography_model.prism_layer.gravity(coordinates, field=\"g_z\")\n",
    "terrain_effect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf76692b-12e7-4494-aa3e-b3f70bd5b364",
   "metadata": {},
   "source": [
    "We store these values in our `DataFrame` \n",
    "\n",
    "Guardemos estos valores en nuestro `DataFrame`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4982e659-7139-4de8-bdd8-5d08dc0abadb",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.assign(terrain_effect_mgal=terrain_effect)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbdc2feb-f04b-43ab-b815-886c2f56c616",
   "metadata": {},
   "source": [
    "We calculate our **Bouguer gravity disturbance** that is defined as the gravity disturbance without the effect of topography. \n",
    "\n",
    "Y calculemos el **disturbio de Bouguer** que se define como el disturbio de gravedad menos el efecto de la topografía"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e7af63d-6519-447e-ac85-daf3e3b5e2a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity_bouguer = data.gravity_disturbance_mgal - data.terrain_effect_mgal\n",
    "\n",
    "data = data.assign(gravity_bouguer_mgal=gravity_bouguer)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69cfae00-215f-45a3-bb52-e02ab409e354",
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter_plot_gravity(\n",
    "    data.easting_m,\n",
    "    data.northing_m,\n",
    "    data.gravity_bouguer_mgal,\n",
    "    title=\"Bouguer gravity disturbance\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "515e14ae-2b9e-4b91-85e1-6ea43e942a63",
   "metadata": {},
   "source": [
    "## Regional-residual separation\n",
    "\n",
    "The Bouguer gravity disturbance is mainly dominated by the effects of very deep masses located at the interface between the crust and the mantle (Mohorovicic discontinuity).\n",
    "If we want to observe the signal from the Bushveld igneous bodies (called **residuals**) we need to remove the effect of these deep masses (called **regional**).\n",
    "\n",
    "For this we are going to use the technique called **equivalent source**:\n",
    "- we will define a set of deep point sources\n",
    "- we will adjust their masses in such a way that they generate a field approximately equal to that of the Bouguer disturbance\n",
    "- we will calculate the **regional** effect as the effect generated by these sources in the observation points\n",
    "\n",
    "<figure>\n",
    "<img src=\"images/equivalent-sources.svg\" alt=\"Sketch equivalent sources method to perform interpolation or predict residuals.\">\n",
    "<figcaption><em>\n",
    "    Sketch equivalent sources method to perform interpolation or predict residuals.\n",
    "    (CC BY).\n",
    "</em></figcaption>\n",
    "</figure>\n",
    "\n",
    "For this we will make use of Harmonica's `EquivalentSources`, which allow us to carry out all these steps through a simple interface to use.\n",
    "\n",
    "--- \n",
    "El disturbio de gravedad de Bouguer está dominado principalmente por los effectos de masas muy profundas ubicadas en la interfaz entre la corteza y el manto (discontinuidad de Mohorovicic).\n",
    "Si queremos observar la señal de los cuerpos ígneos de Bushveld (llamados **residuales**) es necesario que quitemos el efecto de estas masas profundas (llamdas **regionales**).\n",
    "\n",
    "Para ello vamos a utilizar la técnica de **fuentes equivalentes**:\n",
    "- definiremos un conjunto de fuentes puntuales profundas\n",
    "- ajustaremos sus masas de forma tal que generen un campo aproximadamente igual al del disturbio de Bouguer\n",
    "- calcularemos el efecto **regional** como el efecto que generan estas fuentes en los puntos de observación\n",
    "\n",
    "Para ello haremos uso de las `EquivalentSources` de Harmonica, que nos permiten llevar a cabo todos estos pasos mediante una interfaz simple de utilizar."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8302300f-b4a8-42ca-8df3-04835f1aaca6",
   "metadata": {},
   "source": [
    "We define the equivalent sources located at 500 km depth ( the damping parameter lets you smooth the masses of the sources. A larger value generates smoother solutions)\n",
    "\n",
    "Definimos las fuentes equivalentes ubicadas a 500km de profundidad (el damping es un parámetro de amortiguamiento que permite \"suavizar\" las masas de las fuentes, un valor alto genera soluciones más suaves)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a79359e2-584f-4855-b204-abe30279f8f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "deep_sources = hm.EquivalentSources(damping=1000, depth=500e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3f84b79-50e3-45a4-af16-89356ff55804",
   "metadata": {},
   "source": [
    "We fit the deep sources using the Bouguer gravity disturbance values. \n",
    "\n",
    "Ajustamos las masas de las fuentes con los datos del disturbio de Bouguer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e52d4e5f-e2dd-4210-870b-453979374be1",
   "metadata": {},
   "outputs": [],
   "source": [
    "coordinates = (data.easting_m, data.northing_m, data.height_sea_level_m)\n",
    "\n",
    "deep_sources.fit(coordinates, data.gravity_bouguer_mgal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1ca810f-55e8-4599-a6b0-0a7163334d4e",
   "metadata": {},
   "source": [
    "Now we calculate the **regional** effect of the gravity of these deep sources\n",
    "\n",
    "Calculamos el efecto **regional** como el efecto gravitatorio de estas fuentes profundas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01933390-339c-4c47-88eb-2e900e47492e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity_regional = deep_sources.predict(coordinates)\n",
    "\n",
    "data = data.assign(gravity_regional_mgal=gravity_regional)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8716b71d-68b4-4686-9888-db99aabcf92e",
   "metadata": {},
   "source": [
    "We then calculate the **residual** effect by subtracting the regional effect from the Bouguer gravity disturbance\n",
    "\n",
    "Calculamos el efecto **residual** como el disturbio de Bouguer menos el efecto regional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a91962b-bb1d-4253-97ac-40d29473b99f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity_residual = data.gravity_bouguer_mgal - gravity_regional\n",
    "\n",
    "data = data.assign(gravity_residual_mgal=gravity_residual)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388a5dfe-dc06-4e8e-b00d-8b74c2051635",
   "metadata": {},
   "source": [
    "Graficamos un mapa del campo residual. Este mapa debe tener semejanzas con lo que conocemos acerca de la geología de Bushveld.\n",
    "\n",
    "We plot a map of the residual field. This map should bear some resemblance to what we know about the geology of the Bushveld."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e136f41c-cf44-4363-8046-7716bfc5dc49",
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter_plot_gravity(\n",
    "    data.easting_m,\n",
    "    data.northing_m,\n",
    "    data.gravity_residual_mgal,\n",
    "    title=\"Residual gravity\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b18284d-bc7f-497f-8736-75cfe053791c",
   "metadata": {},
   "source": [
    "## Grid the residuals | Grillando los campos residuales\n",
    "\n",
    "One way to improve the display of our results is to interpolate them into a regular grid at constant altitude.\n",
    "\n",
    "The best way to do this is by **equivalent sources**. We are going to use the same technique as in the previous step, but now the sources will not be located at great depths, but rather shallower.\n",
    "Once we fit the coefficients of the sources, we are going to predict their gravitational effects on a regular grid at constant altitude.\n",
    "\n",
    "Una forma de mejorar la visualización de nuestros resultados es interpolarlos en una grilla regular a altitud constante.\n",
    "\n",
    "La mejor forma de hacer esto es mediante **fuentes equivalentes**. Vamos a utilizar la misma técnica que en el paso anterior, pero ahora las fuentes no estarán ubicadas a grandes profundidades, sino más someras.\n",
    "Una vez que ajustemos los coeficientes de las fuentes, vamos a predecir sus efectos gravitatorios en una grilla regular a altitud constante."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "098180bf-efb6-4728-8674-c65b71f8b50c",
   "metadata": {},
   "source": [
    "We define the equivalent sources located at a depth of 10km and with a damping of 10 (we do not want very smooth solutions) and adjust them with the data\n",
    "\n",
    "Definimos las fuentes equivalentes ubicadas a 10km de profundidad y con un damping de 10 (no queremos soluciones muy suaves) y las ajustamos con los datos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90482559-31fb-4af2-8d8c-21c993dadd31",
   "metadata": {},
   "outputs": [],
   "source": [
    "eqs = hm.EquivalentSources(damping=10, depth=10e3)\n",
    "\n",
    "coordinates = (data.easting_m, data.northing_m, data.height_sea_level_m)\n",
    "eqs.fit(coordinates, data.gravity_residual_mgal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88145550-7f1c-46e4-ac1a-58ceb4bea85c",
   "metadata": {},
   "source": [
    "We calculate the effect of the sources on a regular grid at 2200 meters altitude, with a resolution of 2 arc minutes\n",
    "\n",
    "Calculamos el efecto de las fuentes sobre una grilla regular a 2200 metros de altitud, con una resolución de 2 minutos de arco."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc7e7e8f-f4eb-4e91-a5b6-a6b5b71f521a",
   "metadata": {},
   "outputs": [],
   "source": [
    "residual_grid = eqs.grid(\n",
    "    upward=2200,  # height in meters\n",
    "    spacing=3000,  # spacing in meters\n",
    "    data_names=[\"gravity_residual\"],\n",
    ")\n",
    "residual_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f14f07fa-7788-4a63-a6d0-21a174991b45",
   "metadata": {},
   "source": [
    "Graficamos la grilla\n",
    "\n",
    "Plot the grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cabb96b-5fd4-40c8-b4c4-21216ff16748",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 12))\n",
    "residual_grid.gravity_residual.plot()\n",
    "plt.title(\"Residual gravity grid\")\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ca2ed8-c6db-48cf-bdb2-5104a1c0a1e4",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Homework | Tarea\n",
    "\n",
    "**EN**:\n",
    "\n",
    "Now it's your turn to process gravity data!\n",
    "\n",
    "To do that, we are going to use a set of gravity data over Alice Springs, Australia.\n",
    "You can find the data file in https://raw.githubusercontent.com/santisoler/coco-fatiando-tutorial/main/homework_data/australia_gravity.csv.\n",
    "The checksum hash of this file is: `sha256:b0f44e5b523a3ca034d07820022c23e5b65875956a88e6eda4cd69f9b66da9fc`\n",
    "\n",
    "The goal is to write a new notebook that carry out the following tasks:\n",
    "\n",
    "1. Import the packages needed to download, process and plot the data.\n",
    "1. Download the gravity data of Alice Springs using Pooch.\n",
    "1. Download the Earth topography grid in the `earth-topography-10arcmin.nc` (doi: 10.5281/zenodo.5167357) file using Pooch.\n",
    "1. Load the gravity data using Pandas and the topography grid using Xarray.\n",
    "1. Crop the topography grid to the data region plus a 5 degrees padding.\n",
    "    > Help: we can use the `verde.get_region()` function to obtain the data region and then use `verde.pad_region()` to pad it.\n",
    "1. Project the gravity data and the topography grid using a Mercator projection. Use the mean latitude of the data as the central latitude.\n",
    "1. Compute the **gravity disturbance** by removing the normal gravity from the data. To do so, compute the **normal gravity** generated by the WGS84 ellipsoid through Boule.\n",
    "1. Compute the **terrain effect** through Harmonica. You can use the same density values we used in this notebook for the prisms. Then, compute the **Bouguer gravity** as the difference between the **gravity disturbance** and the **terrain effect**.\n",
    "1. Split the Bouguer gravity in the **regional field** and the **residual field**. You can compute the regional field through deep equivalent sources (`harmonica.EquivalentSources`).\n",
    "1. **Grid** the residual field onto a regular grid at a 1300m height with a spacing of 9km.\n",
    "\n",
    "Finally, **create a new** `coco-fatiando-2022` **repo** in your GitHub account and upload the notebook you wrote.\n",
    "\n",
    "You can plot the data while you perform you computations. Besides, we recommend to add every new result to the DataFrame and wrote Markdown cells to briefly explain what is the code doing.\n",
    "\n",
    "**ES**:\n",
    "\n",
    "¡Ahora es tu turno de procesar datos gravimétricos!\n",
    "\n",
    "Para ello vamos a usar un conjunto de datos gravimétricos sobre Alice Springs, Australia.\n",
    "Puedes encontrar este archivo en https://raw.githubusercontent.com/santisoler/coco-fatiando-tutorial/main/homework_data/australia_gravity.csv.\n",
    "La suma de comprobación (hash) de ese archivo es: `sha256:b0f44e5b523a3ca034d07820022c23e5b65875956a88e6eda4cd69f9b66da9fc`\n",
    "\n",
    "El objetivo es escribir un nuevo notebook en el cual se realicen las siguientes tareas:\n",
    "\n",
    "1. Importar las librerías necesarias para descargar, procesar y graficar los datos.\n",
    "1. Descargar los datos gravimétricos de Alice Springs utilizando Pooch.\n",
    "1. Descargar la grilla de topografía presente en el archivo `earth-topography-10arcmin.nc` (doi: 10.5281/zenodo.5167357) utilizando Pooch.\n",
    "1. Cargar los datos gravimétricos con Pandas y la grilla de topografía con Xarray.\n",
    "1. Recortar la grilla de topografía a la región de los datos gravimétricos, incluyendo 5 grados de padding a cada lado.\n",
    "    > Ayuda: podemos utilizar la función `verde.get_region()` para obtener la región de un conjunto de coordenadas, y luego `verde.pad_region()` para agregar padding a la región.\n",
    "1. Proyectar los datos gravimétricos y la grilla de topografía utilizando una proyección de tipo Mercator, con la latitud central ubicada en la latitud media de los datos.\n",
    "1. Calcular el **disturbio de gravedad** quitando la gravedad normal de los datos observados. Para ello, calcula la gravedad normal generada por el elipsoide WGS84 mediante Boule.\n",
    "1. Calcula el **efecto gravitatorio de la topografía** mediante Harmonica. Utiliza la mismas densidades que elegimos en este notebook. Calcula luego la **gravedad de Bouguer** como la diferencia entre el **disturbio de gravedad** y el **efecto topográfico**.\n",
    "1. Separa la gravedad de Bouguer en el **campo regional** y el **campo residual**. Puedes calcular el campo regional mediante fuentes equivalentes (`harmonica.EquivalentSources`) profundas.\n",
    "1. **Grilla** el campo residual en una grilla regular a una altitud de 1300m y con una resolución de 9km.\n",
    "\n",
    "Finalmente, **crea un repositorio** `coco-fatiando-2022` en tu cuenta de GitHub y sube el notebook que escribiste.\n",
    "\n",
    "Puedes graficar los datos a medida que vas realizando las tareas. Además, recomendamos agregar cada nuevo cálculo al DataFrame de datos y escribir celdas Markdown en las que detalles brevemente qué realizas en las celdas de código."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cd8ccf-4404-4384-a4c2-fd91d02535a4",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## License\n",
    "\n",
    "All Python source code is made available under the BSD 3-clause license. You\n",
    "can freely use and modify the code, without warranty, so long as you provide\n",
    "attribution to the authors.\n",
    "\n",
    "Unless otherwise specified, all figures and Jupyter notebooks are available\n",
    "under the Creative Commons Attribution 4.0 License (CC-BY).\n",
    "\n",
    "The full text of these licenses is provided in the [`LICENSE.txt`](LICENSE.txt)\n",
    "file."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}