pycurious/Examples/Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 5 - Mapping Curie depth from EMAG2\n",
    "\n",
    "Using global compilations of the magnetic anomaly, such as EMAG2, means there is no limitation on window sizes. The mapping module in *PyCurious* has several functions to simplify:\n",
    "\n",
    "1. Translating between Coordinate Reference Systems (CRS)\n",
    "2. Griding scattered point data\n",
    "3. Importing and exporting Geotiffs\n",
    "\n",
    "These can be accessed from `pycurious.mapping`\n",
    "\n",
    "This requires some extra dependencies:\n",
    "\n",
    "- pyproj\n",
    "- pyepsg\n",
    "- cartopy (for visualisation)\n",
    "\n",
    "which can be installed via `pip`\n",
    "\n",
    "```shell\n",
    "pip install [--user] pyproj pyepsg cartopy\n",
    "```\n",
    "\n",
    "### Contents\n",
    "\n",
    "- [Import EMAG2 data](#Import-EMAG2-data)\n",
    "- [Interpolate onto grid projection](#Interpolate-onto-grid-projection)\n",
    "- [Compute Curie depth](#Compute-Curie-depth)\n",
    "- [Compare to reference models](#Compare-to-reference-models)\n",
    "- [Export GeoTIFF](#Export-GeoTIFF)"
   ]
  },
  {
   "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\n",
    "from pycurious import mapping\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.io.shapereader as shpreader"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import EMAG2 data\n",
    "\n",
    "The [EMAG2 (v3)](https://www.ngdc.noaa.gov/geomag/emag2.html) Earth magnetic anomaly grid$^1$ has been reformatted from a CSV file to a compressed `.npz` NumPy archive to efficiently load within Python workflows. This dataset can be downloaded from [Zenodo](https://zenodo.org/record/3245551). We have added some download tools accessed under `pycurious.download` to ease the download process and checksum validation.\n",
    "\n",
    "### References\n",
    "\n",
    "$^1$ Brian Meyer, Richard Saltus, Arnaud Chulliat (2017): EMAG2: Earth Magnetic Anomaly Grid (2-arc-minute resolution) Version 3. National Centers for Environmental Information, NOAA. Model. doi:[10.7289/V5H70CVX](https://data.nodc.noaa.gov/cgi-bin/iso?id=gov.noaa.ngdc.mgg.geophysical_models:EMAG2_V3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_resources = [\n",
    "    # EMAG2 v3\n",
    "    {\n",
    "     \"local_file\":\"../../data/EMAG2_V3_20170530.npz\",\n",
    "     \"md5\":'c0898b6a91efb3f13783873a8b67380c',\n",
    "     \"url\":\"https://zenodo.org/record/3245551/files/EMAG2_V3_20170530.npz?download=1\",\n",
    "     \"expected_size\":\"500Mb\"\n",
    "    },\n",
    "    \n",
    "    # Li et al. 2017\n",
    "    {\n",
    "     \"local_file\":\"../../data/Li_et_al_2017.txt\",\n",
    "     \"md5\":'5f0ea0af3e27c6c21cd7b776c979aa80',\n",
    "     \"url\":\"https://media.nature.com/original/nature-assets/srep/2017/170321/srep45129/extref/srep45129-s1.txt?download=1\",\n",
    "     \"expected_size\":\"15Mb\"\n",
    "    },    \n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycurious.download import download_cached_file\n",
    "\n",
    "for resource in file_resources:\n",
    "    print (\"\\nDownloading {:s}\".format(resource[\"local_file\"]))\n",
    "    download_cached_file(resource[\"url\"], resource[\"local_file\"], resource[\"md5\"], resource[\"expected_size\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once loaded, we trim the data to a region of interest - being careful to allow some buffer to accommodate the maximum window size.\n",
    "\n",
    "**EMAG2_V3_20170530.npz** contains the following columns:\n",
    "\n",
    "1. LON ; Geographic Longitude WGS84 (decimal degrees)\n",
    "2. LAT ; Geographic Latitude WGS84 (decimal degrees)\n",
    "3. SeaLevel ; Magnetic Anomaly Value at Sea Level(nT)\n",
    "4. UpCont ; Magnetic Anomaly Value at continuous 4km altitude (nT)\n",
    "5. Code ; Data Source Code (see table below)\n",
    "6. Error ; Error estimate (nT)\n",
    "\n",
    "Code 888 is assigned in certain cells on grid edges where the data source is ambiguous and assigned an error of -888 nT\n",
    "\n",
    "Code 999 is assigned in cells where no data is reported with the anomaly value assigned 99999 nT and an error of -999 nT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emag2_file = file_resources[0][\"local_file\"]\n",
    "with np.load(emag2_file) as f:\n",
    "    mag_data = f['data']\n",
    "\n",
    "# filter NaNs\n",
    "mag_data = mag_data[mag_data[:,3] != 99999.]\n",
    "mag_data = mag_data[mag_data[:,4] != -888.]\n",
    "mag_data = mag_data[mag_data[:,4] != -999.]\n",
    "\n",
    "# print min/max\n",
    "mincols = mag_data.min(axis=0)\n",
    "maxcols = mag_data.max(axis=0)\n",
    "\n",
    "fmt = \"min/max {:5.2f} -> {:5.2f}\"\n",
    "for col in range(mag_data.shape[1]):\n",
    "    print(fmt.format(mincols[col], maxcols[col]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolate onto grid projection\n",
    "\n",
    "We choose an arbitrary section of the Earth's surface and a suitable projection. PyCurious handles transformation between different CRS using a EPSG reference code. EMAG2 is in lons and lats (in decimal degrees), which we want to convert into eastings and northings (in metres). For this example, we transform WGS84 coordinates (EPSG:4326) to the IRENET95 projection (EPSG:2157).\n",
    "\n",
    "> **IMPORTANT:** The power spectrum must be in eastings & northings to compute the power spectrum in rad/km."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grid extent - grid of the British Isles in Irenet 95 grid projection\n",
    "xmin = 0.0\n",
    "xmax = 1900000.0\n",
    "ymin = 0.0\n",
    "ymax = 1700000.0\n",
    "extent_grid = [xmin, xmax, ymin, ymax]\n",
    "\n",
    "# extent on the sphere (WGS84)\n",
    "extent_sphere = mapping.convert_extent(extent_grid, epsg_in=2157, epsg_out=4326)\n",
    "\n",
    "# map extent - should be narrower than extent_sphere\n",
    "extent_map = [-10.8, 2, 49.5, 59]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx, dy = 1e3, 1e3 # 1 km grid resolution\n",
    "nx, ny = int(round((xmax-xmin)/dx)), int(round((ymax-ymin)/dy))\n",
    "\n",
    "# interpolate onto grid extent (for computation)\n",
    "mag_grid = mapping.grid(mag_data[:,:2], mag_data[:,3], extent_grid, shape=(ny,nx), epsg_in=4326, epsg_out=2157)\n",
    "\n",
    "# interpolate onto extent on the sphere (for mapping)\n",
    "mag_sphere = mapping.grid(mag_data[:,:2], mag_data[:,3], extent_sphere, shape=(ny,nx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.UTM(30)\n",
    "\n",
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=proj)\n",
    "ax.set_extent(extent_map)\n",
    "ax.coastlines(resolution='50m', linewidth=1.5)\n",
    "ax.gridlines()\n",
    "\n",
    "im1 = ax.imshow(mag_sphere, extent=extent_sphere, transform=ccrs.PlateCarree(),\n",
    "                cmap='Spectral_r', vmin=-200, vmax=200, zorder=0)\n",
    "\n",
    "fig.colorbar(im1, label='nT')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Curie depth\n",
    "\n",
    "For this example we use the Bouligand *et al.*, 2009 approach we outlined in previous notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)\n",
    "\n",
    "window_size = 400e3\n",
    "\n",
    "# centroid spacing of 50 km\n",
    "xc_list, yc_list = grid.create_centroid_list(window_size, spacingX=50e3, spacingY=50e3)\n",
    "print(\"number of centroids = {}\".format(len(xc_list)))\n",
    "\n",
    "\n",
    "beta, zt, dz, C = grid.optimise_routine(window_size, xc_list, yc_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "curie_depth = zt + dz\n",
    "centroids = np.column_stack([xc_list, yc_list])\n",
    "\n",
    "curie_depth_sphere = mapping.grid(centroids, curie_depth, extent_sphere, (ny,nx), epsg_in=2157, epsg_out=4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "proj = ccrs.UTM(30)\n",
    "\n",
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=proj)\n",
    "ax.set_extent(extent_map)\n",
    "ax.coastlines(resolution='50m', linewidth=1.5)\n",
    "ax.gridlines()\n",
    "\n",
    "im1 = ax.imshow(curie_depth_sphere, extent=extent_sphere, transform=ccrs.PlateCarree(),\n",
    "                cmap='BrBG', vmin=10, vmax=40, zorder=0)\n",
    "\n",
    "fig.colorbar(im1, label='nT')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to reference models\n",
    "\n",
    "A global Curie depth reference model has been created by [Li *et al.*, 2017](https://www.nature.com/articles/srep45129), which is a useful resource to compare these results. Download the supplementary material from the open access [article](https://www.nature.com/articles/srep45129#supplementary-information) to the same directory as EMAG2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "li2017_cpd_file = file_resources[1]['local_file']\n",
    "\n",
    "li_cpd = np.loadtxt(li2017_cpd_file, skiprows=1)\n",
    "li_cpd = li_cpd[~np.isnan(li_cpd[:,2])]\n",
    "\n",
    "li_grid = mapping.grid(li_cpd[:,:2], li_cpd[:,2], extent_sphere, shape=(ny,nx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj = ccrs.UTM(30)\n",
    "\n",
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=proj)\n",
    "ax.set_extent(extent_map)\n",
    "ax.coastlines(resolution='50m', linewidth=1.5)\n",
    "ax.gridlines()\n",
    "\n",
    "im1 = ax.imshow(li_grid, extent=extent_sphere, transform=ccrs.PlateCarree(),\n",
    "                cmap='BrBG', zorder=0, vmin=10, vmax=40)\n",
    "\n",
    "fig.colorbar(im1, label='km')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export GeoTIFF\n",
    "\n",
    "The GeoTIFF format contains projection information, which makes it ideal for GIS applications such as [QGIS](https://www.qgis.org/en/site/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapping.export_geotiff(\"my-curie-depth.tiff\", curie_depth_sphere, extent_sphere, epsg=4326)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also import a GeoTIFF using the `mapping` module. the `import_geotiff` function prints a wealth of information on the reference spheriod and the projection, including EPSG code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GTiff, GTextent = mapping.import_geotiff(\"my-curie-depth.tiff\")\n",
    "\n",
    "print(GTextent)\n",
    "print(GTiff.shape)"
   ]
  }
 ],
 "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
}