.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "tutorial/03_figures/c_geological-map.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_tutorial_03_figures_c_geological-map.py>`
        to download the full example code or to run this example in your browser via Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_tutorial_03_figures_c_geological-map.py:


Geological Map on Topography
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Texture mapping for a GeoTIFF on a topography surface.

To overlay an image/map from a GeoTIFF on to a topography surface, it's necessary to have texture coordinates ("texture mapping") matching the proper extends of the mesh/surface you'd like to drape the texture (GeoTIFF) on.

We can do this by using the spatial reference of the GeoTIFF itself, as this allows you to preserve the entire mesh that the texture is being draped on without having to clip out the parts where you don't have imagery. In this example, we explicitly set the texture extents onto a topography surface where the texture/GeoTIFF has a much larger extent than the topography surface.

Originally posted here: https://github.com/pyvista/pyvista-support/issues/14

.. GENERATED FROM PYTHON SOURCE LINES 13-22

.. code-block:: Python


    import os
    import tempfile

    import numpy as np
    import pyvista as pv
    from pyvista import examples
    import requests








.. GENERATED FROM PYTHON SOURCE LINES 24-28

.. code-block:: Python

    path = examples.download_file("topo_clean.vtk")
    topo = pv.read(path)
    topo






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">

    <table style='width: 100%;'>
    <tr><th>UnstructuredGrid</th><th>Information</th></tr>
    <tr><td>N Cells</td><td>824278</td></tr>
    <tr><td>N Points</td><td>413250</td></tr>
    <tr><td>X Bounds</td><td>3.299e+05, 3.442e+05</td></tr>
    <tr><td>Y Bounds</td><td>4.253e+06, 4.271e+06</td></tr>
    <tr><td>Z Bounds</td><td>1.494e+03, 2.723e+03</td></tr>
    <tr><td>N Arrays</td><td>0</td></tr>
    </table>


    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 29-31

Load the GeoTIFF/texture (this could take a minute to download)
https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0

.. GENERATED FROM PYTHON SOURCE LINES 31-38

.. code-block:: Python

    url = "https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0"

    response = requests.get(url)
    filename = os.path.join(tempfile.gettempdir(), "downsampled_Geologic_map_on_air_photo.tif")
    open(filename, "wb").write(response.content)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    8175934



.. GENERATED FROM PYTHON SOURCE LINES 39-43

In the block below, we can use the ``get_gcps`` function to get the
Ground Control Points of the raster, however this depends on GDAL. For this
tutorial, we are going to hard code the GCPs to avoid having users install
GDAL.

.. GENERATED FROM PYTHON SOURCE LINES 43-64

.. code-block:: Python



    def get_gcps(filename):
        """This helper function retrieves the Ground Control
        Points of a GeoTIFF. Note that this requires gdal"""
        import rasterio

        get_point = lambda gcp: np.array([gcp.x, gcp.y, gcp.z])
        # Load a raster
        src = rasterio.open(filename)
        # Grab the Groung Control Points
        points = np.array([get_point(gcp) for gcp in src.gcps[0]])
        # Now Grab the three corners of their bounding box
        # -- This guarantees we grab the right points
        bounds = pv.PolyData(points).bounds
        origin = [bounds[0], bounds[2], bounds[4]]  # BOTTOM LEFT CORNER
        point_u = [bounds[1], bounds[2], bounds[4]]  # BOTTOM RIGHT CORNER
        point_v = [bounds[0], bounds[3], bounds[4]]  # TOP LEFT CORNER
        return origin, point_u, point_v









.. GENERATED FROM PYTHON SOURCE LINES 65-74

.. code-block:: Python


    # Fetch the GCPs
    # origin, point_u, point_v = get_gcps(filename)

    # Hard code GCPs
    origin = [310967.75148705335, 4238841.045453942, 0.0]
    point_u = [358682.9364281533, 4238841.045453942, 0.0]
    point_v = [310967.75148705335, 4276281.98755258, 0.0]








.. GENERATED FROM PYTHON SOURCE LINES 75-79

.. code-block:: Python


    # Use the GCPs to map the texture coordinates onto the topography surface
    topo.texture_map_to_plane(origin, point_u, point_v, inplace=True)






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table style='width: 100%;'><tr><th>Header</th><th>Data Arrays</th></tr><tr><td>
    <table style='width: 100%;'>
    <tr><th>UnstructuredGrid</th><th>Information</th></tr>
    <tr><td>N Cells</td><td>824278</td></tr>
    <tr><td>N Points</td><td>413250</td></tr>
    <tr><td>X Bounds</td><td>3.299e+05, 3.442e+05</td></tr>
    <tr><td>Y Bounds</td><td>4.253e+06, 4.271e+06</td></tr>
    <tr><td>Z Bounds</td><td>1.494e+03, 2.723e+03</td></tr>
    <tr><td>N Arrays</td><td>1</td></tr>
    </table>

    </td><td>
    <table style='width: 100%;'>
    <tr><th>Name</th><th>Field</th><th>Type</th><th>N Comp</th><th>Min</th><th>Max</th></tr>
    <tr><td>Texture Coordinates</td><td>Points</td><td>float32</td><td>2</td><td>3.737e-01</td><td>8.576e-01</td></tr>
    </table>

    </td></tr> </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 80-81

Show GCPs in relation to topo surface with texture coordinates displayed

.. GENERATED FROM PYTHON SOURCE LINES 81-98

.. code-block:: Python

    p = pv.Plotter()
    p.add_point_labels(
        np.array(
            [
                origin,
                point_u,
                point_v,
            ]
        ),
        ["Origin", "Point U", "Point V"],
        point_size=5,
    )

    p.add_mesh(topo)
    p.show(cpos="xy")





.. image-sg:: /tutorial/03_figures/images/sphx_glr_c_geological-map_001.png
   :alt: c geological map
   :srcset: /tutorial/03_figures/images/sphx_glr_c_geological-map_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 99-100

Read the GeoTIFF as a ``Texture`` in PyVista:

.. GENERATED FROM PYTHON SOURCE LINES 100-113

.. code-block:: Python

    texture = pv.read_texture(filename)

    # Now plot the topo surface with the texture draped over it
    # And make window size large for a high-res screenshot
    p = pv.Plotter(window_size=np.array([1024, 768]) * 3)
    p.add_mesh(topo, texture=texture)
    p.camera_position = [
        (337461.4124956896, 4257141.430658634, 2738.4956020899253),
        (339000.40935731295, 4260394.940646875, 1724.0720826501868),
        (0.10526647627366331, 0.2502863297360612, 0.962432190920575),
    ]
    p.show()




.. image-sg:: /tutorial/03_figures/images/sphx_glr_c_geological-map_002.png
   :alt: c geological map
   :srcset: /tutorial/03_figures/images/sphx_glr_c_geological-map_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 114-121

.. raw:: html

    <center>
      <a target="_blank" href="https://colab.research.google.com/github/pyvista/pyvista-tutorial/blob/gh-pages/notebooks/tutorial/03_figures/c_geological-map.ipynb">
        <img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/ width="150px">
      </a>
    </center>


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 15.007 seconds)


.. _sphx_glr_download_tutorial_03_figures_c_geological-map.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/pyvista/pyvista-tutorial/gh-pages?urlpath=lab/tree/notebooks/tutorial/03_figures/c_geological-map.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: c_geological-map.ipynb <c_geological-map.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: c_geological-map.py <c_geological-map.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_