Topographic Map

Download Python source code | Download Jupyter notebook

Plot aerial imagery from a GeoTIFF on top of a topography mesh.

Similar to the Applying Textures example.

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.25, anti_aliasing=True)
elevation = topo.warp_by_scalar()
# Load the topographic map from a GeoTiff
topo_map = examples.download_crater_imagery()
topo_map = pv.Texture(topo_map.to_image().resample(0.5, anti_aliasing=True))
topo_map = topo_map.flip_y()  # flip to align to our dataset

elevation

StructuredGrid (0x7f533d448340)
  N Cells:      104351
  N Points:     105000
  X Bounds:     1.810e+06, 1.831e+06
  Y Bounds:     5.640e+06, 5.658e+06
  Z Bounds:     7.353e+02, 2.773e+03
  Dimensions:   350, 300, 1
  N Arrays:     1

StructuredGrid 105,000 points 104,351 cells 2.0 MiB

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

Griddims 350 x 300 x 1

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

Let’s inspect the imagery that we just loaded.

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

<matplotlib.image.AxesImage object at 0x7f534170aba0>

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 (0x7f533fb1cd60)
  N Cells:    28187
  N Points:   15041
  X Bounds:   1.818e+06, 1.825e+06
  Y Bounds:   5.645e+06, 5.653e+06
  Z Bounds:   1.376e+03, 2.773e+03
  N Arrays:   2

UnstructuredGrid 15,041 points 28,187 cells 1.3 MiB

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

• Point Data: (2)
  scalar1of1 active
  • scalar1of1⧉
    scalar
    float64
    [1.376e+03, 2.773e+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