Creating a map using x/y data

I have X, Y, Z data where X and Y are lon and lat coordinates and Z is value of interest. All variables are 12500*1800 (too large to attach here).

I would like to plot these variables on map axes so I can add shape files with additional info., aswell as a scale bar.

clf
axesm('MapProjection','mercator') %define axis

surfm(y,x,z) % plot X/Y data

S = shaperead('study_site.shp','UseGeoCoords',true); %add shapefile on top of X/Y data
info = shapeinfo('study_site.shp');
crs = info.CoordinateReferenceSystem
geoshow(S,'FaceColor','white'); 

Apart from needing to adjust the axes limits, this looks good. But, I am struggling to add a scale bar. When I try to do this, the plot is just replaced with a scale bar:

hold on

scaleruler on

How do I do this properly?

Also, how do I add lat and lon information to the frame?