Topographic Map

This is very similar to the Applying Textures example except it is focused on plotting aerial imagery from a GeoTIFF on top of some topography mesh.

from __future__ import annotations

import matplotlib as mpl
import matplotlib.pyplot as plt
import pyvista as pv
from pyvista import examples

Start by loading the elevation data and a topographic map.

# Load the elevation data as a surface
topo = examples.download_crater_topo().resample(0.5, anti_aliasing=True)
elevation = topo.warp_by_scalar()
# Load the topographic map from a GeoTiff
topo_map = examples.download_crater_imagery()
topo_map.to_image().resample(0.5, anti_aliasing=True, inplace=True)
topo_map = topo_map.flip_y()  # flip to align to our dataset

elevation

StructuredGrid (0x720568c08580)
  N Cells:      418701
  N Points:     420000
  X Bounds:     1.810e+06, 1.831e+06
  Y Bounds:     5.640e+06, 5.658e+06
  Z Bounds:     7.353e+02, 2.768e+03
  Dimensions:   700, 600, 1
  N Arrays:     1

StructuredGrid 420,000 points 418,701 cells 8.0 MiB

Bounds⧉X [1.810e+06, 1.831e+06]Y [5.640e+06, 5.658e+06]Z [7.353e+02, 2.768e+03]

Griddims 700 x 600 x 1

• Point Data: (1)
  scalar1of1 active
  • scalar1of1⧉
    scalar
    float64
    [7.353e+02, 2.768e+03]
    active

Let’s inspect the imagery that we just loaded.

mpl.rcParams['figure.dpi'] = 500
plt.imshow(topo_map.to_array())

<matplotlib.image.AxesImage object at 0x720568bf1e80>

Once you have a topography mesh loaded as a surface mesh (we use a pyvista.StructuredGrid here) and an image loaded as a pyvista.Texture using pyvista.read_texture(), then you can map that imagery to the surface mesh as follows:

# Bounds of the aerial imagery - given to us
bounds = (1818000, 1824500, 5645000, 5652500, 0, 3000)
# Clip the elevation dataset to the map's extent
local = elevation.clip_box(bounds, invert=False)
# Apply texturing coordinates to associate the image to the surface
local.texture_map_to_plane(use_bounds=True, inplace=True)

UnstructuredGrid (0x720568c5ffa0)
  N Cells:    110369
  N Points:   57064
  X Bounds:   1.818e+06, 1.825e+06
  Y Bounds:   5.645e+06, 5.653e+06
  Z Bounds:   1.379e+03, 2.768e+03
  N Arrays:   2

UnstructuredGrid 57,064 points 110,369 cells 5.0 MiB

Bounds⧉X [1.818e+06, 1.825e+06]Y [5.645e+06, 5.653e+06]Z [1.379e+03, 2.768e+03]

• Point Data: (2)
  scalar1of1 active
  • scalar1of1⧉
    scalar
    float64
    [1.379e+03, 2.768e+03]
    active
  • Texture Coordinates⧉
    2 comp
    float32
    [0.000e+00, 1.000e+00]
    tcoords

Now display it. Note that the imagery is aligned as we expect.

local.plot(texture=topo_map, cpos='xy')

Static Scene

Interactive Scene

And here is a 3D perspective.

local.plot(texture=topo_map)

Static Scene

Interactive Scene

We could also display the entire region by extracting the surrounding region and plotting the texture mapped local topography and the outside area

# Extract surrounding region from elevation data
surrounding = elevation.clip_box(bounds, invert=True)

# Display with a shading technique
pl = pv.Plotter()
pl.add_mesh(local, texture=topo_map)
pl.add_mesh(surrounding, color='white')
pl.enable_eye_dome_lighting()
pl.camera_position = pv.CameraPosition(
    position=(1831100.0, 5642142.0, 8168.0),
    focal_point=(1820841.0, 5648745.0, 1104.0),
    viewup=(-0.435, 0.248, 0.865),
)
pl.show()

Static Scene

Interactive Scene

Tags: plot