{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "eddf2e87-d54b-4d44-903c-24124afbd262",
   "metadata": {},
   "source": [
    "# Live demo: Processing gravity data with Fatiando a Terra\n",
    "\n",
    "This notebook is based on the [Harmonica tutorial at Transform 2021](https://github.com/fatiando/transform21)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae028f3f-b962-49d0-b401-e93e2eefa897",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "\n",
    "Start by loading everything we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce93d2aa-cefa-4d3a-9518-ab7c806d5cfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The standard Python science stack\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "# For projections (wrapped for Proj)\n",
    "import pyproj\n",
    "\n",
    "# Plotting maps using GMT\n",
    "import pygmt\n",
    "\n",
    "# The Fatiando stack\n",
    "import pooch\n",
    "import verde as vd\n",
    "import boule as bl\n",
    "import harmonica as hm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "617b0940-7288-4460-904c-c65cf498434d",
   "metadata": {},
   "source": [
    "## Data from the Bushveld Igneous Complex (South Africa)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11d465c5-bac1-4113-b568-54684c004a01",
   "metadata": {},
   "source": [
    "We can use [Pooch](https://www.fatiando.org/pooch) to download data files from anywhere on the web. Let's download some public domain gravity data from the Bushveld Igneous Complex that we have on some GitHub repositories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b8c430-73e0-4f98-a0cb-6cfa9d65e30e",
   "metadata": {},
   "outputs": [],
   "source": [
    "url_grav = \"https://github.com/fatiando/2021-gsh/raw/main/data/bushveld_gravity.csv\"\n",
    "md5_grav = \"md5:45539f7945794911c6b5a2eb43391051\"\n",
    "\n",
    "url_topo = \"https://github.com/fatiando/transform21/raw/main/data/bushveld_topography.nc\"\n",
    "md5_topo = \"md5:62daf6a114dda89530e88942aa3b8c41\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1e584ac-f6a5-4dba-becb-51ba58935a4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "print(path_grav)\n",
    "print(path_topo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79d1fcfe-bf5b-437b-8b92-187d7e5b3e7d",
   "metadata": {},
   "source": [
    "Use Pandas to read the gravity data from the CSV file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3b1f231-6a82-4ee9-adac-963adfa4b061",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5a444869-3928-4d74-9240-b0b155b04d43",
   "metadata": {},
   "source": [
    "Use xarray to read the topography data from the netCDF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc6c22e5-662c-47e4-ac9d-a431d5c81ad9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "58426cbb-9c60-47b9-a448-a9a0e959547c",
   "metadata": {},
   "source": [
    "Use pygmt to plot the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca57463-bbbd-4ede-82c3-1603bc7a273a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "pygmt.makecpt(cmap=\"viridis\", series=[data.gravity.min(), data.gravity.max()])\n",
    "fig.plot(\n",
    "    x=data.longitude,\n",
    "    y=data.latitude,\n",
    "    color=data.gravity,\n",
    "    cmap=True,\n",
    "    style=\"c4p\",\n",
    "    projection=\"M15c\", \n",
    "    frame=True,\n",
    ")\n",
    "fig.colorbar(frame='af+l\"observed gravity [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d11271d0-6cc7-4ac3-b567-9e4673ba0cf1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "pygmt.makecpt(cmap=\"earth\", series=[topography.values.min(), topography.values.max()])\n",
    "fig.grdimage(topography, shading=True, projection=\"M15c\", frame=True)\n",
    "fig.colorbar(frame='af+l\"topography [m]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62e3b36c-1a23-4163-8bc3-bb100807fccf",
   "metadata": {},
   "source": [
    "## Gravity disturbance\n",
    "\n",
    "We can use [Boule](https://www.fatiando.org/boule) for computing the normal gravity of the WGS84 reference ellipsoid on any point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b87abd97-99bd-4172-ac7c-c60ccef667c0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7d15e31e-b04f-4762-b309-700b3aae2a2b",
   "metadata": {},
   "source": [
    "And compute the gravity disturbance as the difference between the observed gravity and the normal gravity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8baba02-be0e-4656-b7bb-d41b19232636",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1cf7e04-4dea-4064-9774-a8bf2f7ef86e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "maxabs = vd.maxabs(data.disturbance)\n",
    "pygmt.makecpt(cmap=\"polar\", series=[-maxabs, maxabs])\n",
    "fig.plot(\n",
    "    x=data.longitude,\n",
    "    y=data.latitude,\n",
    "    color=data.disturbance,\n",
    "    cmap=True,\n",
    "    style=\"c4p\",\n",
    "    projection=\"M15c\", \n",
    "    frame=True,\n",
    ")\n",
    "\n",
    "fig.colorbar(frame='af+l\"gravity disturbance [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22095adc-e2f5-4cbf-bd04-ed0cd97bf7f8",
   "metadata": {},
   "source": [
    "## Project the data\n",
    "\n",
    "From here on, we will work in Cartesian coordinates. We need to project our data so that we can do topographic correction and remove trends."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97a3a040-1c3c-49ed-8333-24774d2c1466",
   "metadata": {},
   "source": [
    "Define the Mercator projeciton using `pyproj`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c964b39-1202-41e2-8890-400e06045fa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "projection = pyproj.Proj(proj=\"merc\", lat_ts=data.latitude.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fe562d5-29a4-4389-abdb-a26edb4dc8ab",
   "metadata": {},
   "source": [
    "And use it to project the gravity data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913b9ee8-a825-4a2d-a73a-e3993b2cbbc7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "12031e53-9183-4831-9e5a-c14c644a8fa7",
   "metadata": {},
   "source": [
    "The topography grid is a bit trickier but we can use Verde to do the projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6bfdf31-a373-4311-80c4-a987b31cabd6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6b140e48-f725-4f41-9f4b-261ceb642232",
   "metadata": {},
   "source": [
    "## Topographic correction\n",
    "\n",
    "We can use [Harmonica](https://www.fatiando.org/harmonica) for forward modelling the gravitational effect of the terrain through rectangular prisms.\n",
    "\n",
    "First, define a layer of prisms from the topography grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a936d7eb-865b-47dd-bc19-9a64b9686255",
   "metadata": {},
   "outputs": [],
   "source": [
    "topo_prisms = hm.prism_layer(\n",
    "    coordinates=,\n",
    "    surface=,\n",
    "    reference=,\n",
    "    properties={\"density\": }\n",
    ")\n",
    "topo_prisms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b65a62da-cc21-471e-8e96-1a9341b8623c",
   "metadata": {},
   "source": [
    "Second, forward model the gravitational effect of the terrain at the data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b221dac-736f-41b9-89da-099af0ed27f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "598f5d79-b55e-4b7c-984d-c9df170c5443",
   "metadata": {},
   "source": [
    "Calculate the Bouguer disturbance (topography-free)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0de1433-5b55-43fb-a5eb-efbfb5a37ccc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "39e485ff-5a03-4507-ad59-97a4ecd4ed73",
   "metadata": {},
   "source": [
    "Plot it with pygmt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d47db642-5830-4a4e-b67a-99e8fe061e8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "maxabs = vd.maxabs(data.bouguer)\n",
    "pygmt.makecpt(cmap=\"polar\", series=[-maxabs, maxabs])\n",
    "fig.plot(\n",
    "    x=data.longitude,\n",
    "    y=data.latitude,\n",
    "    color=data.bouguer,\n",
    "    cmap=True,\n",
    "    style=\"c4p\",\n",
    "    projection=\"M15c\", \n",
    "    frame=True,\n",
    ")\n",
    "fig.colorbar(frame='af+l\"Bouguer disturbance [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "705be8cf-af67-45e8-ad5b-08e58fe90d5e",
   "metadata": {},
   "source": [
    "## Regional-residual separation\n",
    "\n",
    "We can use [Verde](https://www.fatiando.org/verde) to remove a second degree trend from the Bouguer disturbance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ff6631c-1a14-4a13-947c-d2011873f6b9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e9bbb1-9fb8-450c-8e51-4e6b47934413",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "maxabs = np.quantile(np.abs(data.residuals), 0.99)\n",
    "pygmt.makecpt(cmap=\"polar\", series=[-maxabs, maxabs])\n",
    "fig.plot(\n",
    "    x=data.longitude,\n",
    "    y=data.latitude,\n",
    "    color=data.residuals,\n",
    "    cmap=True,\n",
    "    style=\"c5p\",\n",
    "    projection=\"M15c\", \n",
    "    frame=True,\n",
    ")\n",
    "fig.colorbar(frame='af+l\"residual field [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35dc8dcb-d53b-45e7-8c27-93a629575fab",
   "metadata": {},
   "source": [
    "## Grid the residuals with Equivalent Sources\n",
    "\n",
    "Use the equivalent sources implementation in Harmonica to grid the residuals and upward-continue them to the same height (all in one step)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b177b8-85be-4000-9e22-4e81398a3f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae412802-c7de-4636-aa73-cc1c8dce4d58",
   "metadata": {},
   "source": [
    "Use the source model to forward model the grid at a uniform height. We can use the projection to generate a grid in **geographic coordinates** instead of Cartesian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44a837ef-e3c6-4071-90c5-e7e77a7f7a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "region = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99855692-3e69-430e-88f2-1c414fe2d80a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "grid = eql.grid(\n",
    "    upward=,\n",
    "    region=,\n",
    "    spacing=,\n",
    "    projection=,\n",
    "    data_names=[\"residuals\"],\n",
    "    dims=(\"latitude\", \"longitude\"),\n",
    ")\n",
    "\n",
    "grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "113e1df9-55c2-4057-8d39-17db7b3e210c",
   "metadata": {},
   "source": [
    "Plot the gridded residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "327a4f1d-c2e4-4c66-8c8e-be567377066e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "scale = np.quantile(np.abs(grid.residuals), 0.995)\n",
    "pygmt.makecpt(cmap=\"polar\", series=[-scale, scale], no_bg=True)\n",
    "fig.grdimage(\n",
    "    grid.residuals,\n",
    "    shading=\"+a45+nt0.15\",\n",
    "    cmap=True,\n",
    "    projection=\"M15c\",\n",
    "    frame=True,\n",
    ")\n",
    "fig.colorbar(frame='af+l\"residual field [mGal]\"')\n",
    "fig.plot(\n",
    "    x=data.longitude,\n",
    "    y=data.latitude,\n",
    "    style=\"c0.05c\",\n",
    "    color=\"gray\",\n",
    ")\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:2021-gfz]",
   "language": "python",
   "name": "conda-env-2021-gfz-py"
  },
  "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}