examples/user_guide/Analyzing_Meshes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Very often there will be specific regions of a set of results that are the most of interest, and in those cases it can be helpful to pull out a particular region's data for visualization in a clearer way.  The ``LineCrossSection`` and ``SurfaceCrossSection`` classes allow taking cross-sections of a mesh (or some other rasterizable object) along poly-lines drawn using the ``PolyDraw`` tool, letting you make an *x*/*z* or *x*/*z*/*t* plot from data laid out in *x*/*y*/*z* or *x*/*y*/*z*/*t* (respectively). These classes rasterize the data at a fixed resolution, and then sample the data at that resolution given the underlying polylines. The ``resolution`` may be defined in the units of the underlying coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "import geoviews as gv\n",
    "import cartopy.crs as ccrs\n",
    "from colorcet import cm_n\n",
    "\n",
    "from holoviews import opts\n",
    "\n",
    "from earthsim.analysis import LineCrossSection, SurfaceCrossSection\n",
    "from earthsim.io import read_3dm_mesh, read_mesh2d\n",
    "\n",
    "hv.extension('bokeh')\n",
    "\n",
    "opts.defaults(\n",
    "    opts.Curve(height=400, width=400, framewise=True),\n",
    "    opts.Image(width=500, height=400, colorbar=True, framewise=True),\n",
    "    opts.NdOverlay(legend_limit=0),\n",
    "    opts.Path(line_width=3, color='black'),\n",
    "    opts.RGB(width=500, height=400),\n",
    "    opts.VLine(color='black')\n",
    ")\n",
    "hv.output(holomap='scrubber')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: A static TriMesh\n",
    "\n",
    "As a first simple example, we will pass a static TriMesh to the ``LineCrossSection`` and display it. We'll pass a couple of sample polylines to start with, and if a live Python process is available we will be able to see the depth along the paths appear on the right, for those paths and any more that are subsequently drawn by the user."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../data/Chesapeake_and_Delaware_Bays.3dm'\n",
    "\n",
    "cb_paths = [[(-8523594.,  4588993.), (-8476533.,  4578872.),\n",
    "             (-8449931.,  4562126.), (-8435409.,  4539747.)],\n",
    "            [(-8477265.,  4544109.), (-8469725.,  4430387.)]]\n",
    "            \n",
    "tris, verts = read_3dm_mesh(filename)\n",
    "points = gv.operation.project_points(gv.Points(verts, vdims=['z']))\n",
    "trimesh = hv.TriMesh((tris, points))\n",
    "\n",
    "sector1 = LineCrossSection(trimesh, cb_paths)\n",
    "sector1.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each individual polyline (thick black connected line segments) on the left should result in one curve on the right, with a colored dot on that path corresponding to the color of the curve.  The location of the dot will change as you hover in the plot on the right, allowing you to see which value in the curve corresponds to which location along the polyline.\n",
    "\n",
    "## Example 2: A time-varying mesh\n",
    "\n",
    "``LineCrossSection`` also allows working with time-varying meshes. Here we will modify the static TriMesh data with a time-varying random offset to demonstrate that the ``LineCrossSection`` also works for data that is evolving temporally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sd_paths = [\n",
    "    [(-13037161.,   3843454.), (-13041435.,   3854275.), (-13045822.,   3857996.),\n",
    "     (-13048664.,   3857210.), (-13050429.,   3855575.), (-13048385.,   3849452.),\n",
    "     (-13048027.,   3847896.)],\n",
    "    [(-13042975.,   3850040.), (-13063919.,   3844636.)]\n",
    "]\n",
    "\n",
    "filename = '../data/SanDiego_Mesh/SanDiego.3dm'\n",
    "tris, verts = read_3dm_mesh(filename, skiprows=2)\n",
    "points = gv.operation.project_points(gv.Points(verts, vdims=['z'], crs=ccrs.UTM(11)))\n",
    "trimesh = gv.TriMesh((tris, points))\n",
    "\n",
    "filename2 = '../data/SanDiego_Mesh/SanDiego_ovl.dat'\n",
    "dfs = read_mesh2d(filename2)\n",
    "\n",
    "points = gv.operation.project_points(gv.Points((verts.x, verts.y), crs=ccrs.UTM(11)))\n",
    "\n",
    "def time_mesh(time):\n",
    "    depth_points = points.add_dimension('Velocity', 0, dfs[time].values[:, 0], vdim=True)\n",
    "    return gv.TriMesh((tris, depth_points), crs=ccrs.GOOGLE_MERCATOR)\n",
    "\n",
    "time_dim = hv.Dimension('Time', values=sorted(dfs.keys()), default=3600)\n",
    "meshes = hv.DynamicMap(time_mesh, kdims=time_dim)\n",
    "\n",
    "sector2 = LineCrossSection(meshes, sd_paths, resolution=100)\n",
    "sector2.view(cmap=cm_n.rainbow_r, shade=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the scrubber controls to animate both plots, showing both the map-based data and the curve cross sections over time.\n",
    "\n",
    "## Example 3: Sampling a Surface\n",
    "\n",
    "Instead of sampling individual Curve elements for each time as in example 2, we can take cross-sections across time, returning an Image plotting the sampled value across both distance and time. The ``SurfaceCrossSection`` class is a simple subclass of ``LineCrossSection`` that overrides the ``sample`` method to apply the sampling for all time values and return an ``Image``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sector3 = SurfaceCrossSection(meshes, sd_paths[:1], resolution=100)\n",
    "sector3.view(cmap=cm_n.rainbow_r, shade=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scrubber will now only animate the plot on the left, as time has been laid out vertically on the plot on the right. \n",
    "\n",
    "If you have a running Python process and want to draw your own line, first delete the existing one, because an overlay of non-transparent images will only show the one on top."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}