{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Geological Map on Topography\n\nTexture mapping for a GeoTIFF on a topography surface.\n\nTo overlay an image/map from a GeoTIFF on to a topography surface, it\\'s\nnecessary to have texture coordinates (\\\"texture mapping\\\") matching the\nproper extends of the mesh/surface you\\'d like to drape the texture\n(GeoTIFF) on.\n\nWe can do this by using the spatial reference of the GeoTIFF itself, as\nthis allows you to preserve the entire mesh that the texture is being\ndraped on without having to clip out the parts where you don\\'t have\nimagery. In this example, we explicitly set the texture extents onto a\ntopography surface where the texture/GeoTIFF has a much larger extent\nthan the topography surface.\n\nOriginally posted here:\n<https://github.com/pyvista/pyvista-support/issues/14>\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nimport tempfile\n\nimport numpy as np\nimport pyvista as pv\nimport requests\nfrom pyvista import examples\n\ntry:\n import rasterio\nexcept ImportError:\n rasterio = None"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"path = examples.download_file(\"topo_clean.vtk\")\ntopo = pv.read(path)\ntopo"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the GeoTIFF/texture (this could take a minute to download)\n<https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0>\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"url = \"https://dl.dropbox.com/s/bp9j3fl3wbi0fld/downsampled_Geologic_map_on_air_photo.tif?dl=0\"\n\nresponse = requests.get(url) # noqa: S113\nfilename = os.path.join(tempfile.gettempdir(), \"downsampled_Geologic_map_on_air_photo.tif\") # noqa: PTH118\nopen(filename, \"wb\").write(response.content) # noqa: SIM115, PTH123"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the block below, we can use the `get_gcps` function to get the Ground\nControl Points of the raster, however this depends on GDAL. For this\ntutorial, we are going to hard code the GCPs to avoid having users\ninstall GDAL.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def get_gcps(filename):\n \"\"\"\n Helper function retrieves the Ground Control\n Points of a GeoTIFF. Note that this requires gdal.\n \"\"\"\n if rasterio is None:\n msg = \"rasterio is required for this function\"\n raise ImportError(msg)\n\n def get_point(gcp):\n return np.array([gcp.x, gcp.y, gcp.z])\n\n # Load a raster\n src = rasterio.open(filename)\n # Grab the Groung Control Points\n points = np.array([get_point(gcp) for gcp in src.gcps[0]])\n # Now Grab the three corners of their bounding box\n # -- This guarantees we grab the right points\n bounds = pv.PolyData(points).bounds\n origin = [bounds[0], bounds[2], bounds[4]] # BOTTOM LEFT CORNER\n point_u = [bounds[1], bounds[2], bounds[4]] # BOTTOM RIGHT CORNER\n point_v = [bounds[0], bounds[3], bounds[4]] # TOP LEFT CORNER\n return origin, point_u, point_v"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Fetch the GCPs\n# origin, point_u, point_v = get_gcps(filename)\n\n# Hard code GCPs\norigin = [310967.75148705335, 4238841.045453942, 0.0]\npoint_u = [358682.9364281533, 4238841.045453942, 0.0]\npoint_v = [310967.75148705335, 4276281.98755258, 0.0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Use the GCPs to map the texture coordinates onto the topography surface\ntopo.texture_map_to_plane(origin, point_u, point_v, inplace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Show GCPs in relation to topo surface with texture coordinates displayed\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"pl = pv.Plotter()\npl.add_point_labels(\n np.array(\n [\n origin,\n point_u,\n point_v,\n ]\n ),\n [\"Origin\", \"Point U\", \"Point V\"],\n point_size=5,\n)\n\npl.add_mesh(topo)\npl.show(cpos=\"xy\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Read the GeoTIFF as a `Texture` in PyVista:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"texture = pv.read_texture(filename)\n\n# Now plot the topo surface with the texture draped over it\n# And make window size large for a high-res screenshot\npl = pv.Plotter(window_size=np.array([1024, 768]) * 3)\npl.add_mesh(topo, texture=texture)\npl.camera_position = [\n (337461.4124956896, 4257141.430658634, 2738.4956020899253),\n (339000.40935731295, 4260394.940646875, 1724.0720826501868),\n (0.10526647627366331, 0.2502863297360612, 0.962432190920575),\n]\npl.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{=html}\n<center>\n <a target=\"_blank\" href=\"https://colab.research.google.com/github/pyvista/pyvista-tutorial/blob/gh-pages/notebooks/tutorial/03_figures/c_geological-map.ipynb\">\n <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/ width=\"150px\">\n </a>\n</center>\n```\n"
]
}
],
"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.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}