pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1 - Plot power spectra\n",
    "\n",
    "`PyCurious` offers convenience functions to extract the radial power spectrum. This has been covered in the Bouligand folder of jupyter notebooks and is largely duplicated here.\n",
    "\n",
    "In this notebook we plot the radial power spectrum using the analytical expression, and the spectrum computed from a synthetic magnetic anomaly. (This can be generated from `Bouligand_forward.py` in the `tests` directory.)\n",
    "\n",
    "### Contents\n",
    "\n",
    "- [Radial power spectrum](#Power-spectrum-from-FFT)\n",
    "- [Azimuthal power spectrum](#Azimuthal-power-spectrum)"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "## 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.\n",
    "\n",
    "`CurieGrid` achieves the following purposes:\n",
    "\n",
    "- Upward continuation\n",
    "- Reduction to the pole\n",
    "- Compute the radial power spectrum using FFT\n",
    "\n",
    "The shape of the radial power spectrum is heavily dependent on window size. Resolution of long wavelength features require large windows.\n",
    "\n",
    "> Suggestion: use a window size **> 4 times** the maximum Curie depth.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load magnetic anomaly - i.e. random fractal noise\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot random fractal noise\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "im1 = ax1.imshow(d)\n",
    "fig.colorbar(im1, label='nT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialise CurieGrid object\n",
    "grid = pycurious.CurieGrid(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 = 304e3\n",
    "subgrid = grid.subgrid(window_size, xpt, ypt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute radial power spectrum\n",
    "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid)\n",
    "\n",
    "# plot radial power spectrum\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111, xlabel=\"wavenumber (rad/km)\", ylabel=\"radial power spectrum\")\n",
    "ax1.plot(k, Phi, '-o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Choice of taper**\n",
    "\n",
    "The default taper is the hanning filter (see [`numpy.hanning`](#hanning) for more details), but other functions can be passed to taper, or simply set it to `None`. There is a significant offset in the power spectrum with tapering functions, however, it is the *slope* of the spectrum with wavenumber that is most important when it comes to determining Curie depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "k1, Phi1, sigma_Phi1 = grid.radial_spectrum(subgrid, taper=None)\n",
    "k2, Phi2, sigma_Phi2 = grid.radial_spectrum(subgrid, taper=np.hanning)\n",
    "k3, Phi3, sigma_Phi2 = grid.radial_spectrum(subgrid, taper=np.hamming)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111, xlabel=\"wavenumber (rad/km)\", ylabel=\"radial power spectrum\")\n",
    "ax1.plot(k1, Phi1, '-o', label='none')\n",
    "ax1.plot(k2, Phi2, '-o', label='hanning')\n",
    "ax1.plot(k3, Phi3, '-o', label='hamming')\n",
    "ax1.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Azimuthal power spectrum\n",
    "\n",
    "The azimuthal spectrum computes the FFT on a square window that is polarised at a range of radii. Subdividing the transforms into bins by azimuth is useful to explore linear trends in the magnetic anomaly that may align with a particular foliation or strike orientation.\n",
    "\n",
    "```python\n",
    "azimuthal_spectrum(subgrid, power=2.0, theta=5.0)\n",
    "```\n",
    "\n",
    "`theta` controls the bin size of each azimuth (in degrees). The FFT of the magnetic anomaly is raised to the `power=2`(default) which is compatible with Bouligand _et al._ (2009) computation of Curie depth. For Tanaka *et al.* (1999), however, we need to take the square root which is set with `power=0.5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k, Phi, theta = grid.azimuthal_spectrum(subgrid, power=0.5, theta=20.0)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111, ylabel=r'azimuth $\\theta$', xlabel='wavenumber (rad/km)')\n",
    "im1 = ax1.pcolor(k, theta, Phi)\n",
    "fig.colorbar(im1, label='radial power spectrum')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting each azimuth on the power spectrum-wavenumber axis highlights the different in slope for different polarisations of the magnetic anomaly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111, xlabel=\"wavenumber (rad/km)\", ylabel=\"radial power spectrum\")\n",
    "\n",
    "for i, angle in enumerate(theta):\n",
    "    ax1.plot(k, Phi[i], label=r'$\\theta = {}^\\circ$'.format(angle))\n",
    "\n",
    "ax1.legend(bbox_to_anchor=(1,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}