Plot 3D surface from .tiff - why does it look so strange?

Question

I am trying to visualize two .tiffs as 3D surface plots with matplotlib. The outcome somehow looks strange - can anyone help me why?

# import relevant moduls
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# Specify input files
dem_old = '../data/Nisqually_1987.tif'
dem_new = '../data/Nisqually_2017.tif'

# Creating first 3D figure
# set up a figure twice as wide as it is tall
fig = plt.figure(figsize=plt.figaspect(0.5))

# Old DEM plot
ax = fig.add_subplot(1, 2, 1, projection='3d')
dem = gdal.Open(dem_old)
dem_array = dem.ReadAsArray()

lin_x = np.linspace(0,1,dem_array.shape[0],endpoint=False)
lin_y = np.linspace(0,1,dem_array.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y,lin_x)
z = dem_array
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=15)
ax.set_title('3D Nisqually-Gletscher Oberfläche (1987)')
plt.xticks([])
plt.yticks([])

# New DEM 3D plot
ax = fig.add_subplot(1, 2, 2, projection='3d')
dem_2 = gdal.Open(dem_new)
dem_array_2 = dem_2.ReadAsArray()

lin_x_2 = np.linspace(0,1,dem_array_2.shape[0],endpoint=False)
lin_y_2 = np.linspace(0,1,dem_array_2.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y_2,lin_x_2)
z = dem_array_2
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_title('3D Nisqually-Gletscher Oberfläche (2017)')
plt.xticks([])
plt.yticks([])

# show plot
plt.savefig ('glaciers_3D.png')
plt.show()

Data was obtained from: https://www.sciencebase.gov/catalog/item/6140d0a3d34e1449c5d6000a

Answer 1

After checking the values in the tiff and as shown in your figure, I found the normal values are between [1000, 3000], but some values are a constant of -3.40282e+38, so just change these values to np.nan by doing this :

dem_array = dem.ReadAsArray()
dem_array[abs(dem_array) > 1e10] = np.nan

dem_array_2 = dem_2.ReadAsArray()
dem_array_2[abs(dem_array_2) > 1e10] = np.nan

And then you will get the figure :

But in the figure of 2017, I found some "holes" in the surface plot, So I made a linear interpolation for filling the "holes" with pandas :

import pandas as pd 

z = dem_array_2
df = pd.DataFrame.from_records(z)
df_padded  = df.interpolate(limit_area='inside')
z = df_padded.values

We get the new figure :