The function doit is written to return what gscatter returns, but gscatter returns an array of line handles - one line handle per group plotted - so doit returns a non-scalar array whenever there is more than one group plotted.
The solution is to make doit return a scalar no matter what. For instance, doit can return a scalar cell array containing the line handles returned from gscatter. (Or, if you don't need those line handles once the lines are created, you can just have doit return a scalar number, say 0 or 1 or NaN or anything else, or have doit return a status, say, logical true indicating it completed successfully.)
tT=readtable("data1.csv", "VariableNamingRule","preserve");
tT.Properties.VariableNames(2)={'Site'}; % shorten to be more convenient to use
G=grpstats(tT,{'Site','Depth'},{'mean','median','std'},'DataVars',{'Date','d18O'});
hSc=rowfun(@doit,tT,'GroupingVariables',{'Site'},'InputVariables',{'Date','d18O','Site','Depth'},'OutputFormat','uniform');
disp(hSc)
function out=doit(x,y,s,d)
figure
h=gscatter(x,y,d);
out = {h}; % return a scalar cell array containing the lines
xlabel('Date'),ylabel('d18O')
title(sprintf('Site: %d',s(1))) % title the plot with the site value
end
