examples/00-load/create-surface-draped.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline\nfrom pyvista import set_plot_theme\nset_plot_theme('document')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Drape 2D Surface From Line {#create_draped_surf_example}\n==========================\n\nDrape a surface (2D array) from a line in 3D space.\n\nThis is a common task to create a 2.5D image/sectional mesh of data like\nGPR or airborne EM profiles (geophysics applications). This example\nprovides a look into how to create a 2.5D sectional mesh from typical\ndata in those use cases.\n\nFor this example, we have an instrument path on the ground surface (the\nline) and a 2D array of the collected image under that line.\n\nOriginally posted in [this support\nissue](https://github.com/pyvista/pyvista-support/issues/135).\n\nSuppose you have some GPR data (or anything that produces a line of data\nwith values at depth). With these data, you\\'ll have a 2D image/array of\nyour data values and 3D coordinates of where that line/profile is in 3D\nspace (often where you collected the data on the surface of topography).\nAttached below are some example data for this: 1) XYZ coordinates of a\nGPR path and 2) a 2D array of data values produced from the GPR.\n\nthe data here are wacky (it\\'s difficult to get shareable data of decent\nquality), so ignore them but pay attention to the structure. The\ncoordinates we have are technically shifted up and we have some NaN\nfiller above the surface - its weird and just ignore it. You\\'ll\ntypically have a more uniform looking profile in 2D with the coordinates\nassociated to the top of each column in your 2D array.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pyvista as pv\nfrom pyvista import examples\n\n# Extract the data archive and load these files\n# 2D array of XYZ coordinates\npath = examples.download_gpr_path().points\n# 2D array of the data values from the imaging equipment\ndata = examples.download_gpr_data_array()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(15, 3))\nplt.pcolormesh(data, cmap=\"seismic\", clim=[-1, 1])\nplt.gca().invert_yaxis()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "View the path of the GPR profile from a top-down perspective. Since we\nhave the full coordinates (XY and Z), we can create a structured mesh\n\\\"draping\\\" down from those coordinates to hold the GPR image data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.scatter(path[:, 1], path[:, 0])\nplt.axis(\"image\")\nplt.xlabel(\"Northing\")\nplt.ylabel(\"Easting\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "assert len(path) in data.shape, \"Make sure coordinates are present for every trace.\"\n# If not, you'll need to interpolate the path\n\n# Grab the number of samples (in Z dir) and number of traces/soundings\nnsamples, ntraces = data.shape  # Might be opposite for your data, pay attention here\n\n# Define the Z spacing of your 2D section\nz_spacing = 0.12\n\n# Create structured points draping down from the path\npoints = np.repeat(path, nsamples, axis=0)\n# repeat the Z locations across\ntp = np.arange(0, z_spacing * nsamples, z_spacing)\ntp = path[:, 2][:, None] - tp\npoints[:, -1] = tp.ravel()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a StructuredGrid from the structured points\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = pv.StructuredGrid()\ngrid.points = points\ngrid.dimensions = nsamples, ntraces, 1\n\n# Add the data array - note the ordering\ngrid[\"values\"] = data.ravel(order=\"F\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And now we can plot it, process it, or do anything, because it is a\nPyVista mesh and the possibilities are endless with PyVista.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cpos = [\n    (1217002.366883762, 345363.80666238244, 3816.828857791056),\n    (1216322.4753436751, 344033.0310674846, 3331.052985309526),\n    (-0.17716571330686096, -0.25634368781817973, 0.9502106207279767),\n]\n\np = pv.Plotter()\np.add_mesh(grid, cmap=\"seismic\", clim=[-1, 1])\np.add_mesh(pv.PolyData(path), color='orange')\np.show(cpos=cpos)"
      ]
    }
  ],
  "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.11.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}