adding a geotiff as a texture to unstructured grid

[craigmillernz]

I'm trying to add a geotiff (a topo map) to an unstructured grid of the topography DEM.

topo = pv.read(topo_file.vtk)
topo_warp  = topo.warp_by_scalar()
topo_clip = topo_warp.clip_box(bounds=clipbox, invert=False)

topo_map = pv.read_texture(topo_map.tif)
p = pv.BackgroundPlotter()
p.add_mesh(topo_clip, name='topo_surface', cmap='gist_gray', 
                      opacity=0.3,  texture=topo_map)

However it produces this error. I've clipped the topo_map geotiff to the same extent as the topo dem.

AssertionError: Input mesh does not have texture coordinates to support the texture.

topo_map
(vtkRenderingOpenGL2Python.vtkOpenGLTexture)0000018E1D9C3DC8

I also tried adding texture_map_to_plane to the topo_clip

topo_clip.texture_map_to_plane(inplace=True)

which worked but the geotiff was incorrectly scaled to the underlying topography surface.

Can you advise the proper way to add a geotiff to an unstructured grid such as a DEM?

[craigmillernz]

ok, my problem was initially with the wrong clip of the DEM and secondly with not doing an inplace=True on the texture_map_to_plane

[banesullivan]

Hi @craigmillernz, I'm glad to see that you figured this out on your own! To be thorough, I thought I'd add a bit more details on texture mapping for a GeoTIFF on a topography surface - also, I do this ALOT so I have plenty of examples.

As you found out, it's necessary to have texture coordinates ("texture mapping") matching the proper extents of the mesh/surface you'd like to drape the texture (GeoTIFF) on.

I commonly 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. Unfortunately, I don't have a good example to demo that particular case, but I do have an example where we explicitly set the texture extents onto a topography surface where the texture/GeoTIFF has a much larger extent than the topography surface.

The data files for this example are quite large (~300MB) so I have them on my Dropbox here

Here's the example:

import pyvista as pv
import numpy as np
import gdal

topo = pv.read('topo_clean.vtk')
topo

UnstructuredGrid	Information
N Cells	824278
N Points	413250
X Bounds	3.299e+05, 3.442e+05
Y Bounds	4.253e+06, 4.271e+06
Z Bounds	1.494e+03, 2.723e+03
N Arrays	0

# Load the GeoTIFF/texture
filename = 'Geologic_map_on_air_photo.tif'

def get_gcps(filename):
    """This helper function retrieves the Ground Control 
    Points of a GeoTIFF. Note that this requires gdal"""
    get_point = lambda gcp : np.array([gcp.GCPX, gcp.GCPY, gcp.GCPZ])
    # Load a raster
    ds = gdal.Open(filename)
    # Grab the Groung Control Points
    points = np.array([get_point(gcp) for gcp in ds.GetGCPs()])
    # 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

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

# Show GCPs in relation to topo surface with texture coordinates displayed
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')

# Read the GeoTIFF as a vtkTexture:
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.show(cpos='xy')

Please note that you do not need to use GDAL to get the ground control points (GCPs) - all you need to know are the coordinates of the BOTTOM LEFT CORNER, BOTTOM RIGHT CORNER, and TOP LEFT CORNER of the GeoTIFF/texture to feed to the texture_map_to_plane method to properly build the texture coordinates.

[bluetyson]

I actually get p = pv.Plotter() p.add_point_labels(np.array([origin, point_u, point_v,]), ['Origin', 'Point U', 'Point V'], point_size=5, color='red')

p.add_mesh(elevation) p.show(cpos='xy')

when trying this

TypeError Traceback (most recent call last)

in 1 # Show GCPs in relation to topo surface with texture coordinates displayed 2 p = pv.Plotter() ----> 3 p.add_point_labels(np.array([origin, point_u, point_v,]), 4 ['Origin', 'Point U', 'Point V'], 5 point_size=5, color='red') TypeError: add_point_labels() got an unexpected keyword argument 'color'

[banesullivan]

The API has changed since that was posted and now that argument is point_color to be more explicit about what is getting colored because the text and the box can also have color (text_color and shape_color). Please see the documentation for this method

[banesullivan]

I edited the original snippet to remove that

[bluetyson]

Thanks!

[bluetyson]

Ok, so what happens if you take your topography map and download the puppy texture example to it?

[banesullivan]

@bluetyson, can you please open a new issue for that.

I'd recommend just trying it. The puppy texture would just map over the topography via the UV texture coordinates produced by texture_map_to_plane for a playful and adorable topography surface

[bluetyson]

Will do. Obviously not an overly serious example, but I was just wondering. e.g. I have a silly case where I have a mesozoic basement surface - and was wondering, what if I put a 'mesozoic monster' on it. However, small resolution image, large scale regional area - it would just stretch out?

[banesullivan]

I have a mesozoic basement surface - and was wondering, what if I put a 'mesozoic monster' on it.

😂 Please share the data for this in a new issue, I'd have fun putting an example together for it

small resolution image, large scale regional area - it would just stretch out?

The image will stretch but should still look okay.

[GGDRriedel]

I'm trying to map a geotif to a dummy plane so I can see something over a geotif of OSM data.

Trying to create that plane:

topo_map = pv.read_texture('geotiff.tif)

#coordinates are known
x = np.arange(403926.1008, 418211.4231, 10)
y = np.arange(5470594.0881, 5484179.9408, 10)
x, y = np.meshgrid(x, y)
#height is local height minus 100 for viewing
z=np.ones(x.shape)+323-100
dummyplane= pv.StructuredGrid(x, y, z)

p = pvqt.BackgroundPlotter()
p.add_mesh(grid, cmap="seismic", clim=[-1*per,per])
p.add_mesh(pv.PolyData(path), color='orange')
p.add_mesh(dummyplane,texture=topo_map)
#p.add_floor()
p.show()

However I'm also getting the error. Input mesh does not have texture coordinates to support the texture.

Any way to customky set these coordinates for the mesh without using gdal and all of the other stuff?

[GGDRriedel]

Nevermind.

dummyplane.texture_map_to_plane(origin=[x[0],y[0],323-100],
                             point_u =[x[1],y[0],323-100],
                             point_v=[x[0],y[1],323-100],
                             inplace=True)

[GGDRriedel]

welp, this does not throw errors, however my texure doesn't show. using

p.add_mesh(dummyplane)
#p.add_floor()
p.show()

Without texture sows the plane.

Do the units(mesh cellseize) and pixes of the texture need to align?

[banesullivan]

You need to pass the texture to the add_mesh method. For further help, please open a new support ticket.