pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2 - Compute Curie depth\n",
    "\n",
    "Bouligand *et al.* (2009) formulated an expression for the radial power assuming a fractal model for crustal magnetisation (see [Ex1-Plot-power-spectra](./Ex1-Plot-power-spectrum.ipynb) for analytic solution), which depends on 4 parameters:\n",
    "\n",
    "- $\\beta$ - a fractal parameter\n",
    "- $z_t$ - top of magnetic sources\n",
    "- $\\Delta z$ - thickness of the magnetic layer\n",
    "- $C$ - a field constant\n",
    "\n",
    "These parameters can be fitted to the radial power spectrum computed from FFT to determine Curie depth ($z_\\mathrm{curie} = z_t + \\Delta z$).\n",
    "\n",
    "### Contents\n",
    "\n",
    "- [Plot radial power spectrum](#Plot-radial-power-spectrum)\n",
    "- [Optimise parameters](#Optimise-parameters)\n",
    "- [Varying window sizes](#Varying-window-sizes)\n",
    "- [Add prior constraints](#Add-prior-constraints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import pycurious"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load x,y,anomaly\n",
    "mag_data = np.loadtxt(\"../../data/test_mag_data.txt\")\n",
    "\n",
    "nx, ny = 305, 305\n",
    "\n",
    "x = mag_data[:,0]\n",
    "y = mag_data[:,1]\n",
    "d = mag_data[:,2].reshape(ny,nx)\n",
    "\n",
    "xmin, xmax = x.min(), x.max()\n",
    "ymin, ymax = y.min(), y.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot radial power spectrum\n",
    "\n",
    "The radial power spectrum is computed from a square window of the magnetic anomaly. Methods to select window sizes and compute the Fast Fourier Transform (FFT) belong to the `CurieGrid` object. We apply the default `np.hanning` taper to the power spectrum as in [Ex1-Plot-power-spectrum](#./Ex1-Plot-power-spectrum.ipynb).\n",
    "\n",
    "By default the FFT of the magnetic anomaly is raised to the power 2:\n",
    "\n",
    "```python\n",
    "grid.radial_spectrum(subgrid, taper=None, power=2.0, **kwargs)\n",
    "```\n",
    "\n",
    "which is compatible with Bouligand _et al._ (2009) computation of Curie depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = pycurious.CurieOptimise(d, xmin, xmax, ymin, ymax)\n",
    "\n",
    "# pick the centroid\n",
    "xpt = 0.5*(xmin + xmax)\n",
    "ypt = 0.5*(ymin + ymax)\n",
    "\n",
    "window_size = 200e3\n",
    "subgrid = grid.subgrid(window_size, xpt, ypt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid)\n",
    "\n",
    "# plot radial power spectrum\n",
    "plt.plot(k, Phi, '-o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimise parameters\n",
    "\n",
    "The `optimise` routine determines the Curie parameters ($\\beta, z_t, \\Delta z, C$) by minimising the misfit between the analytic power spectum and the computed power spectrum. Optimisation is handles by the nonlinear least-squares algorithm implemented in `scipy.optimize.minimize`.\n",
    "\n",
    "```python\n",
    "optimise(window, xc, yc, beta=3.0, zt=1.0, dz=10.0, C=5.0)\n",
    "```\n",
    "\n",
    "These are sane defaults to initiate the optimisation, however, values close to the final solution will result in faster convergence and stability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta, zt, dz, C = grid.optimise(window_size, xpt, ypt)\n",
    "print(\"beta = {:.2f}\\n  zt = {:.2f}\\n  dz = {:.2f}\\n   C = {:.2f}\".format(beta,zt,dz,C))\n",
    "\n",
    "Phi2 = pycurious.bouligand2009(k, beta, zt, dz, C)\n",
    "\n",
    "# plot radial power spectrum - compare analytic against computed\n",
    "plt.plot(k, Phi, '-o')\n",
    "plt.plot(k, Phi2, linewidth=2)\n",
    "\n",
    "print(\"\\nCurie depth is {:.2f} km\".format(zt + dz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Varying window sizes\n",
    "\n",
    "Here, we observe the effect decreasing the window size has on the optimisation of Curie parameters. The number of points in the power spectrum decrease with smaller window sizes, particularly at low wavenumbers (long wavelengths). This significantly reduces the quality of Curie depth determinations for window sizes less than 200 km.\n",
    "\n",
    "As a rule of thumb, the window size should be > 4 times the maximum possible Curie depth that could be resolved from a study area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# window size in kilometres\n",
    "window_size_range = [304e3, 200e3, 100e3, 50e3]\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "for i, window_size in enumerate(window_size_range):\n",
    "    \n",
    "    # compute radial power spectrum\n",
    "    subgrid = grid.subgrid(window_size, xpt, ypt)\n",
    "    k, Phi, sigma_Phi = grid.radial_spectrum(subgrid)\n",
    "    ax1.plot(k, Phi, '-o', color='C{}'.format(i))\n",
    "    \n",
    "    # optimise Curie variables\n",
    "    beta, zt, dz, C = grid.optimise(window_size, xpt, ypt)\n",
    "    print(\"beta = {:.2f}\\n  zt = {:.2f}\\n  dz = {:.2f}\\n   C = {:.2f}\\n\".format(beta,zt,dz,C))\n",
    "    \n",
    "    Phi2 = pycurious.bouligand2009(k, beta, zt, dz, C)\n",
    "    ax1.plot(k, Phi2, color='C{}'.format(i), linewidth=2, label=\"window size = {}km\".format(window_size/1e3))\n",
    "\n",
    "ax1.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add prior constraints\n",
    "\n",
    "Spurious Curie depth values can often be obtained if windows of the magnetic anomaly do not contain long wavelengths. This often arises from selecting too small a small window size, but can also occur from uncorrected artefacts in the magnetic anomaly, or simply the absence of magnetic minerals close to the surface of the crust. Adding *a priori* information on $\\beta, z_t, \\Delta z, C$ can reduce the incidence of anomalous values if there is reasonable grounds to do so.\n",
    "\n",
    "> **IMPORTANT:** strict priors can unfairly bias the results - *only* add priors if you have good reason to do so. (This is **not recommended** for most use cases.)\n",
    "\n",
    "```python\n",
    "from scipy import stats\n",
    "\n",
    "# define normal distributions\n",
    "beta_p = stats.norm(3.0, 1.0)\n",
    "zt_p = stats.norm(0.0, 1.0)\n",
    "dz_p = stats.norm(10.0, 8.0)\n",
    "C_p = stats.norm(10.0, 10.0)\n",
    "\n",
    "# add to grid object\n",
    "grid.add_prior(beta=beta_p, zt=zt_p, dz=dz_p, C=C_p)\n",
    "```\n",
    "\n",
    "Priors can be accessed from a dictionary:\n",
    "\n",
    "```python\n",
    "prior = grid.prior_pdf['beta'] # stats.norm object\n",
    "prior = grid.prior['beta'] # stats.norm object arguments (3.0, 1.0)\n",
    "```\n",
    "\n",
    "A combination of priors can be added. `reset_priors()` removes all prior information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.reset_priors()\n",
    "\n",
    "from scipy import stats\n",
    "beta_p = stats.norm(3.0, 0.01)\n",
    "grid.add_prior(beta=beta_p)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111, xlabel=\"wavenumber (rad/km)\", ylabel=\"radial power spectrum\")\n",
    "\n",
    "for i, window_size in enumerate(window_size_range):\n",
    "    \n",
    "    # compute radial power spectrum\n",
    "    subgrid = grid.subgrid(window_size, xpt, ypt)\n",
    "    k, Phi, sigma_Phi = grid.radial_spectrum(subgrid)\n",
    "    ax1.plot(k, Phi, '-o', color='C{}'.format(i), label='{} km window size'.format(window_size*1e-3))\n",
    "    \n",
    "    # optimise Curie variables\n",
    "    beta, zt, dz, C = grid.optimise(window_size, xpt, ypt)\n",
    "    print(\"beta = {:.2f}\\n  zt = {:.2f}\\n  dz = {:.2f}\\n   C = {:.2f}\\n\".format(beta,zt,dz,C))\n",
    "    \n",
    "    Phi2 = pycurious.bouligand2009(k, beta, zt, dz, C)\n",
    "    ax1.plot(k, Phi2, color='C{}'.format(i), linewidth=2)\n",
    "    \n",
    "ax1.legend()"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}