Digital Elevation Map (DEM) with texture

[banesullivan]

Description

Originally post on Slack by user Andrew Annex

hey all, checking out vtki to replace vispy code I have laying around, ran into a texture issue, I have a DEM as a numpy array and an grayscale image overlay I want on the texture but there is some sort of rotation or scaling issue that I can’t seem to solve by using np.rot90, I will post screen shots in just a sec

dem no texture

with texture, looks like from this view it needs to be rotated 180 degrees and not repeat

Example Data

I am using the texture map function…

dem, img # load from npy files
x = np.arange(0, dem.shape[0])
y = np.arange(0, dem.shape[1])
mx, my = np.meshgrid(x,y, indexing='ij')
grid = vtki.StructuredGrid(mx,my,dem)
grid.texture_map_to_plane(inplace=True)
img = data['img'] * 255
img = np.dstack((img,img,img)).astype(np.uint8)
print(dem.shape, img.shape, img.min(), img.max(), img.dtype,)
texture = vtki.numpy_to_texture(img.transpose((1,0,2)))
grid.plot(texture=texture, cmap='gray')

Archive.zip

okay that should be it… I am also running into some silly issue with the gray scale as I moved away from using the matplotlib/cmap code, it is probably just a nan issue

[banesullivan]

NaN values plot as solid grey by default - you could change the nan_color argument when plotting or change it in the vtki.rcParams['nan_color'] = 'yellow'

[banesullivan]

Andrew Annex:

I take that back, I did have masked arrays but there are no masked values, other dems will have them though so it is good to know
but yeah I don’t get the rotation issue, there should be a 1:1 correspondence with each dem pixel and each img pixel

doing: “grid.plot(texture=vtki.numpy_to_texture(np.rot90(img)))” gets me closer, but there is still that odd repeat at either end of the mesh
thanks for the quick responses though, I need to go but can check back in later

[banesullivan]

So real quick… are you sure you want to map your data onto the mesh as a texture? If your data values have:

a 1:1 correspondence with each dem pixel and each img pixel

then you might actually want to put the data on the cells (“pixels”) of your mesh and plot it as data values - texture mapping is meant for plotting images that don’t have a 1:1 correspondence with the mesh discretization.

Try this for example:

import numpy as np
import vtki
# Panel doesn't handle structured grids very well
vtki.rcParams['use_panel'] = False

# load the data
dem = np.load('dem_vtki.npy')
rimg = np.load('img_vtki.npy')

# A nice default view direction for this data
cpos = [-1,1,1]

# Make a mesh from the DEM
x = np.arange(0, dem.shape[0])
y = np.arange(0, dem.shape[1])
mx, my = np.meshgrid(x,y, indexing='ij')
grid = vtki.StructuredGrid(mx,my,dem)#.elevation()

# Add image data to the mesh with 1:1 correspondance
grid['dem'] = rimg.ravel(order='F')
grid.plot(cpos=cpos, cmap='gray')

Which produces:

[banesullivan]

Andrew Annex:

ahh, that works, I am not familiar with all of the terminology yet so I may have made that more confusing than it needed to bbe

so as a follow on question, say if I have the same image but with 4x the resolution, can I “drape” that texture over the lower resolution cells or would I need to interpolate the elevation data to match the dimensions?

[banesullivan]

That’s when texture mapping makes sense!

[banesullivan]

Solves problem originally posted by Andrew Annex - then interpolates the given image to a high resolution to showcase texture mapping a high res image to a low res mesh

import numpy as np
import vtki
# Panel doesn't handle structured grids very well
vtki.rcParams['use_panel'] = False

# load the data
dem = np.load('dem_vtki.npy')
rimg = np.load('img_vtki.npy')

x = np.arange(0, dem.shape[0])
y = np.arange(0, dem.shape[1])
mx, my = np.meshgrid(x,y, indexing='ij')

grid = vtki.StructuredGrid(mx,my,dem)
# Make sure to map texture on a flat plane
#   defaults to best fitting plane which is not 
#   what you need
e = grid.bounds
grid.texture_map_to_plane(origin=[e[0], e[2], 0],
                          point_u=[e[1], e[3], 0],
                          point_v=[e[0], e[3], 0],
                          inplace=True)

norm = lambda x: (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x))
img = norm(rimg) * 255
img = np.dstack((img,img,img)).astype(np.uint8)
img = np.flip(np.moveaxis(np.transpose(img), 0, -1), 0)
print(dem.shape, img.shape, img.min(), img.max(), img.dtype,)

texture = vtki.numpy_to_texture(img)
grid.plot(texture=texture, cpos=[-1,1,1])

##############################################

# Interpolate the image
from scipy.interpolate import griddata 

# Make a mesh from the DEM
xi = np.linspace(0, dem.shape[0], num=dem.shape[0]*4)
yi = np.linspace(0, dem.shape[1], num=dem.shape[1]*4)
mxi, myi = np.meshgrid(xi, yi)

result = griddata(grid.points[:,0:2], 
                  rimg.ravel(order='F'), 
                  tuple((mxi, myi)), 
                  method='linear')

# Make texture:
norm = lambda x: (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x))
result = norm(result) * 255

fake_image = result.reshape((mxi.shape[0], mxi.shape[1]), order='F')
# fake_image = np.swapaxes(fake_image, 0, 1)

# fake_image = np.flip(fake_image, 1)
fake_image = np.flip(fake_image, 0)
img = np.dstack((fake_image,fake_image,fake_image)).astype(np.uint8)


print(img.shape, img.min(), img.max(), img.dtype,)
import matplotlib.pyplot as plt
plt.imshow(img)
plt.show()

texture = vtki.numpy_to_texture(img)
grid.plot(texture=texture, cpos=[-1,1,1], )