Error: shape mismatch when plotting certain DEM Geotiffs in 3D - using 'plot_surface' from Matplotlib

I’m trying to plot various DEM geotiff rasters in Python 3.7 using plot_surface in matplotlib but encountered an error that I can’t solve. I have 2 DEM rasters, the small one plots correctly, the other doesn’t.

Here is the plot of the small raster:

Here is the error I get with the larger raster:

line 30: 
surf = ax.plot_surface(X,Y,Z, rstride=1, zorder = 0, cstride=1, cmap=plt.cm.RdYlBu_r,  linewidth=0, antialiased=False)

..axes3d.py", line 1504, in plot_surface
    X, Y, Z = np.broadcast_arrays(X, Y, Z)

..stride_tricks.py", line 258, in broadcast_arrays
    shape = _broadcast_shape(*args)

..stride_tricks.py", line 189, in _broadcast_shape
    b = np.broadcast(*args[:32])

ValueError: shape mismatch: objects cannot be broadcast to a single shape

The two geotiffs are downloadable here.

Here is the script:

import numpy as np
import gdal
import geopandas 
from geopandas import GeoSeries
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from shapely.geometry import shape


fig, ax = plt.subplots(figsize=(8,5), subplot_kw={'projection': '3d'})
plt.ioff()

demFile = r"C:Tempsmall_DEM_8bit.tif"     # works
demFile = r"C:TempDEM_8bit.tif"           # fails

ds = gdal.Open(demFile)
gt  = ds.GetGeoTransform()
dem = ds.ReadAsArray()

xres = gt[1]
yres = gt[5]
X = np.arange(gt[0], gt[0] + dem.shape[1]*xres, xres)
Y = np.arange(gt[3], gt[3] + dem.shape[0]*yres, yres)

X, Y = np.meshgrid(X, Y)
Z = dem

print("X: {}. Y: {}. Z: {}".format(len(X),len(Y),len(Z)))

surf = ax.plot_surface(X,Y,Z, rstride=1, zorder = 0, cstride=1, cmap=plt.cm.RdYlBu_r,  linewidth=0, antialiased=False)
print("plotted surface")
ax.set_zlim(0, 100)
ax.view_init(30,-105)
ax.set_zlim3d(-1.01, 200)

fig.colorbar(surf, shrink=0.4, aspect=20)
fig.show()

Both rasters are 8 bit unsigned geotiffs, both produced by gdal with the same extract process. Obviously the second raster is larger but still only 681KB and the dimensions of the X,Y,Z match each other for both the small and large tif. I tried transposing the Z (np.transpose(Z)) but it didn’t work either.

Please can someone put me out of my misery as to why I get the error?