Thanks for the reply Thomas, I think I did it in a slightly different way.

Here is my code:

```
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from mpl_toolkits.mplot3d import axes3d
# Normalize vector
def normalize(X):
return X / (1e-16 + np.sqrt((np.array(X) ** 2).sum(axis=-1)))[..., np.newaxis]
# get data
X, Y, Z = axes3d.get_test_data(0.02)
fig, ax = plt.subplots()
cs = plt.contour(X, Y, Z, levels=14, colors='None')
plt.contourf(X, Y, Z, levels=14, cmap="Blues")
# Get contour segments
segs = cs.allsegs
# number of contour levels
num_levels = len(cs.allsegs)
# LightSource direction 225deg
ls = np.array([-np.sqrt(2)/2, -np.sqrt(2)/2])
# Loop through levels
for lvl in range(num_levels):
# Loop through elements on each level
for el in range(len(cs.allsegs[lvl])):
# Get segment from level=lvl and element=el
xn = cs.allsegs[lvl][el][:, 0]
yn = cs.allsegs[lvl][el][:, 1]
points = np.array([xn, yn]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-2], points[1:-1],points[2:]], axis=1)
# Compute normal vector from segments
dy = segments[:, 1, 1] - segments[:, 0, 1]
dx = segments[:, 1, 0] - segments[:, 0, 0]
N = normalize(np.column_stack([-dy,dx]))
# Compute cosine of angle theta between normal and lightsource
costheta = -np.dot(N, ls)
# LineWidths based on cos(theta)
lwidths = np.abs(0.7 * costheta)
lc = LineCollection(segments, linewidths=lwidths, cmap="Greys", norm=plt.Normalize(-1, 1))
# Color based on cosine
lc.set_array(costheta)
ax.add_collection(lc)
plt.savefig('tanaka.png', dpi=600)
```

And here is the resulting image:

It works pretty well, perhaps two points of improvement might be:

- getting the contour line collection without having to plot the contour
- the linecollection isn’t a smooth curve so if you zoom in you might see discontinuities between each segment.

Thanks again.