notebooks/Conversion_between_lon-lat-depth_and_x-y-z.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conversion between lon/lat/depth and x/y/z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve issue with autocomplete\n",
    "%config Completer.use_jedi = False\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`geotree` can read lon/lat/depth or x/y/z as inputs. Here is a list of relevant functions:\n",
    "- `add_lonlatdep` (depth should be in meters; positive depths specify points inside the Earth.)\n",
    "- `add_lonlatdep_query` (same as above except for queries)\n",
    "- `add_xyz` (in meters)\n",
    "- `add_xyz_q` (for queries, in meters)\n",
    "\n",
    "In this section, we show two functions in geotree: `lonlatdep2xyz_spherical` and `xyz2lonlatdep_spherical`. \n",
    "These are used internally to convert between lon/lat/dep and x/y/z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geotree import convert as geoconvert\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define a set of lons/lats/depths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npoints = 100\n",
    "lons = np.random.randint(-180, 180, npoints)\n",
    "lats = np.random.randint(-90, 90, npoints)\n",
    "depths = np.zeros(npoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## lons/lats/depths ---> x/y/z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we use geoconvert.lonlatdep2xyz_spherical to convert lons/lats/depths ---> x/y/z (in meters)\n",
    "# Note that We set depths to zeros, i.e., all points are on a sphere with a radius of 6371000 meters.\n",
    "x, y, z = geoconvert.lonlatdep2xyz_spherical(lons, \n",
    "                                             lats, \n",
    "                                             depths, \n",
    "                                             return_one_arr=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Left pabel: in geographic coordinate, right panel: x/y/z in meters (on a sphere with radius of 6371000m)\n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "\n",
    "plt.scatter(lons, lats,\n",
    "            c=\"k\", \n",
    "            marker=\"o\",\n",
    "            zorder=100)\n",
    "\n",
    "plt.xlabel(\"lons\", size=20)\n",
    "plt.ylabel(\"lats\", size=20)\n",
    "plt.xticks(size=14); plt.yticks(size=14)\n",
    "plt.xlim(-180, 180); plt.ylim(-90, 90)\n",
    "plt.grid()\n",
    "\n",
    "# ---\n",
    "ax = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "\n",
    "ax.scatter3D(x, y, z, c=\"k\", marker=\"o\");\n",
    "\n",
    "ax.set_xlabel('X (m)', size=16)\n",
    "ax.set_ylabel('Y (m)', size=16)\n",
    "ax.set_zlabel('Z (m)', size=16)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## x/y/z ---> lons/lats/depths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just as test, we now use geoconvert.xyz2lonlatdep_spherical to convert x/y/z back to lons/lats/depths\n",
    "lons_conv, lats_conv, depths_conv = geoconvert.xyz2lonlatdep_spherical(x, y, z, \n",
    "                                                                       return_one_arr=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and, we measure the L1 error between original lons/lats/depths and the ones computed above:\n",
    "print(max(abs(lons - lons_conv)))\n",
    "print(max(abs(lats - lats_conv)))\n",
    "print(max(abs(depths - depths_conv)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (vwm_py37)",
   "language": "python",
   "name": "vwm_py37"
  },
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}