I'm fitting a gam with `mgcv` and plot the result with the default `plot.gam()` function. My model includes a 2D-smoother and I want to plot the result as a heatmap. Is there any way to add a colorbar for the heatmap?

I've previously looked into other GAM potting packages, but none of them provided the necessary visualisation. Please note, this is just a simplification for illustration purposes; the actual model (and reporting needs) is much more complicated

**edited:** I initially had swapped y and z in my tensor product, updated to reflect the correct version both in the code and the plot

```
df.gam<-gam(y~te(x,z), data=df, method='REML')
plot(df.gam, scheme=2, hcolors=heat.colors(999, rev =T), rug=F)

```

[![2D GAM heatmap](https://i.stack.imgur.com/BMIZM.png)](https://i.stack.imgur.com/BMIZM.png)

sample data:

```
structure(list(x = c(3, 17, 37, 9, 4, 11, 20.5, 11.5, 16, 17, 
18, 15, 13, 29.5, 13.5, 25, 15, 13, 20, 20.5, 17, 11, 11, 5, 
16, 13, 3.5, 16, 16, 5, 20.5, 2, 20, 9, 23.5, 18, 3.5, 16, 23, 
3, 37, 24, 5, 2, 9, 3, 8, 10.5, 37, 3, 9, 11, 10.5, 9, 5.5, 8, 
22, 15.5, 18, 15, 3.5, 4.5, 20, 22, 4, 8, 18, 19, 26, 9, 5, 18, 
10.5, 30, 15, 13, 27, 19, 5.5, 18, 11.5, 23.5, 2, 25, 30, 17, 
18, 5, 16.5, 9, 2, 2, 23, 21, 15.5, 13, 3, 24, 17, 4.5), z = c(144, 
59, 66, 99, 136, 46, 76, 87, 54, 59, 46, 96, 38, 101, 84, 64, 
92, 56, 69, 76, 93, 109, 46, 124, 54, 98, 131, 89, 69, 124, 105, 
120, 69, 99, 84, 75, 129, 69, 74, 112, 66, 78, 118, 120, 103, 
116, 98, 57, 66, 116, 108, 95, 57, 41, 20, 89, 61, 61, 82, 52, 
129, 119, 69, 61, 136, 98, 94, 70, 77, 108, 118, 94, 105, 52, 
52, 38, 73, 59, 110, 97, 87, 84, 119, 64, 68, 93, 94, 9, 96, 
103, 119, 119, 74, 52, 95, 56, 112, 78, 93, 119), y = c(96.535, 
113.54, 108.17, 104.755, 94.36, 110.74, 112.83, 110.525, 103.645, 
117.875, 105.035, 109.62, 105.24, 119.485, 107.52, 107.925, 107.875, 
108.015, 115.455, 114.69, 116.715, 103.725, 110.395, 100.42, 
108.79, 110.94, 99.13, 110.935, 112.94, 100.785, 110.035, 102.95, 
108.42, 109.385, 119.09, 110.93, 99.885, 109.96, 116.575, 100.91, 
114.615, 113.87, 103.08, 101.15, 98.68, 101.825, 105.36, 110.045, 
118.575, 108.45, 99.21, 109.19, 107.175, 103.14, 94.855, 108.15, 
109.345, 110.935, 112.395, 111.13, 95.185, 100.335, 112.105, 
111.595, 100.365, 108.75, 116.695, 110.745, 112.455, 104.92, 
102.13, 110.905, 107.365, 113.785, 105.595, 107.65, 114.325, 
108.195, 96.72, 112.65, 103.81, 115.93, 101.41, 115.455, 108.58, 
118.705, 116.465, 96.89, 108.655, 107.225, 101.79, 102.235, 112.08, 
109.455, 111.945, 104.11, 94.775, 110.745, 112.44, 102.525)), row.names = c(NA, 
-100L), class = "data.frame")

```