{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "eddf2e87-d54b-4d44-903c-24124afbd262",
   "metadata": {},
   "source": [
    "# Live demo: Processing gravity data with Fatiando a Terra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae028f3f-b962-49d0-b401-e93e2eefa897",
   "metadata": {},
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce93d2aa-cefa-4d3a-9518-ab7c806d5cfe",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pygmt\n",
    "import pyproj\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "\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": [
    "## Load Bushveld Igneous Complex gravity data (South Africa) and a DEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b8c430-73e0-4f98-a0cb-6cfa9d65e30e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "url = \"https://github.com/fatiando/2021-gsh/main/raw/notebook/data/bushveld_gravity.csv\"\n",
    "md5_hash = \"md5:45539f7945794911c6b5a2eb43391051\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1e584ac-f6a5-4dba-becb-51ba58935a4a",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3b1f231-6a82-4ee9-adac-963adfa4b061",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv(fname)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca57463-bbbd-4ede-82c3-1603bc7a273a",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Obtain the region to plot using Verde ([W, E, S, N])\n",
    "region_deg = vd.get_region((data.longitude, data.latitude))\n",
    "\n",
    "fig = pygmt.Figure()\n",
    "fig.basemap(projection=\"M15c\", region=region_deg, frame=True)\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",
    ")\n",
    "fig.colorbar(frame='af+l\"Observed Gravity [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "989eb172-c7bc-4d02-b779-a3f3943ae6db",
   "metadata": {},
   "source": [
    "Let's download a DEM for the same area:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75d2cfb6-3972-4c12-97d5-38c042d047b7",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "url = \"https://github.com/fatiando/transform21/raw/main/data/bushveld_topography.nc\"\n",
    "md5_hash = \"md5:62daf6a114dda89530e88942aa3b8c41\"\n",
    "\n",
    "fname = pooch.retrieve(url, known_hash=md5_hash, fname=\"bushveld_topography.nc\")\n",
    "fname"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bce4e958-d43c-4be8-b252-5fc4997c41f7",
   "metadata": {},
   "source": [
    "And use Xarray to load the netCDF file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01551e59-7f77-4996-8068-dd0b26ec01de",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "topography = xr.load_dataarray(fname)\n",
    "topography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d11271d0-6cc7-4ac3-b567-9e4673ba0cf1",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot topography using pygmt\n",
    "topo_region = vd.get_region((topography.longitude.values, topography.latitude.values))\n",
    "\n",
    "fig = pygmt.Figure()\n",
    "topo_region = vd.get_region((topography.longitude.values, topography.latitude.values))\n",
    "fig.basemap(projection=\"M15c\", region=topo_region, frame=True)\n",
    "\n",
    "vmin, vmax = topography.values.min(), topography.values.max()\n",
    "pygmt.makecpt(cmap=\"batlow\", series=[vmin, vmax])\n",
    "fig.grdimage(topography)\n",
    "\n",
    "fig.colorbar(frame='af+l\"Topography [m]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62e3b36c-1a23-4163-8bc3-bb100807fccf",
   "metadata": {},
   "source": [
    "## Compute gravity disturbance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b87abd97-99bd-4172-ac7c-c60ccef667c0",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8baba02-be0e-4656-b7bb-d41b19232636",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "data[\"disturbance\"] = data.gravity - normal_gravity\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1cf7e04-4dea-4064-9774-a8bf2f7ef86e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "fig.basemap(projection=\"M15c\", region=region_deg, frame=True)\n",
    "\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",
    ")\n",
    "\n",
    "fig.colorbar(frame='af+l\"Gravity disturbance [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f66973a9-b078-49e0-9997-d1f5843ed580",
   "metadata": {},
   "source": [
    "## Remove terrain correction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22095adc-e2f5-4cbf-bd04-ed0cd97bf7f8",
   "metadata": {},
   "source": [
    "### Project the data to plain coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c964b39-1202-41e2-8890-400e06045fa2",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "projection = pyproj.Proj(proj=\"merc\", lat_ts=data.latitude.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913b9ee8-a825-4a2d-a73a-e3993b2cbbc7",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dfef710-efd5-4add-9a95-89d23b0939df",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "data[\"easting\"] = easting\n",
    "data[\"northing\"] = northing\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12031e53-9183-4831-9e5a-c14c644a8fa7",
   "metadata": {},
   "source": [
    "### Project the topography to plain coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6bfdf31-a373-4311-80c4-a987b31cabd6",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e769f871-af99-4669-9232-ea6c1d42c231",
   "metadata": {},
   "source": [
    "### Compute gravitational effect of the layer of prisms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e67fbeba-e48e-4be0-a318-d18c246f9c3b",
   "metadata": {},
   "source": [
    "Create a model of the terrain with prisms\n",
    "\n",
    "![](images/prisms.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a936d7eb-865b-47dd-bc19-9a64b9686255",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b65a62da-cc21-471e-8e96-1a9341b8623c",
   "metadata": {},
   "source": [
    "Calculate the gravitational effect of the terrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b221dac-736f-41b9-89da-099af0ed27f1",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "598f5d79-b55e-4b7c-984d-c9df170c5443",
   "metadata": {},
   "source": [
    "Calculate the Bouguer disturbance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0de1433-5b55-43fb-a5eb-efbfb5a37ccc",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "data[\"bouguer\"] = data.disturbance - terrain_effect\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41559110-927f-4913-9678-b720e7ca501e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "fig.basemap(projection=\"M15c\", region=region_deg, frame=True)\n",
    "\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",
    ")\n",
    "fig.colorbar(frame='af+l\"Bouguer gravity disturbance [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "705be8cf-af67-45e8-ad5b-08e58fe90d5e",
   "metadata": {},
   "source": [
    "## Calculate residuals\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": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48dbc668-2196-48eb-9644-0701cfdecfdc",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "data[\"residuals\"] = residuals\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e9bbb1-9fb8-450c-8e51-4e6b47934413",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "fig.basemap(projection=\"M15c\", region=region_deg, frame=True)\n",
    "\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",
    ")\n",
    "fig.colorbar(frame='af+l\"Residuals [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35dc8dcb-d53b-45e7-8c27-93a629575fab",
   "metadata": {},
   "source": [
    "## Grid the residuals with Equivalent Sources\n",
    "\n",
    "We can use [Harmonica](https://www.fatiando.org/harmonica) to grid the residuals though the equivalent sources technique\n",
    "\n",
    "![](images/eql.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b177b8-85be-4000-9e22-4e81398a3f77",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99855692-3e69-430e-88f2-1c414fe2d80a",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "327a4f1d-c2e4-4c66-8c8e-be567377066e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = pygmt.Figure()\n",
    "fig.basemap(projection=\"M15c\", region=region_deg, frame=True)\n",
    "\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",
    ")\n",
    "fig.colorbar(frame='af+l\"Residuals [mGal]\"')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ab85f6b-46b6-453d-a953-3213f47ed610",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "source": [
    "![Bushveld geologic map](images/study-area.png)\n",
    "\n",
    "Susan J. Webb, R. Grant Cawthorn, Teresia Nguuri, David James; Gravity modeling of Bushveld Complex connectivity supported by Southern African Seismic Experiment results. South African Journal of Geology 2004;; 107 (1-2): 207–218. doi: https://doi.org/10.2113/107.1-2.207"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}