You may try code below that works around StatsPlots to produce:
Plots / StatsPlots code
using DataFrames, StatsPlots; gr(dpi=600, size=(600,400))
stream = ["Glide", "Pool", "Riffle"]
colors = [:red, :blue, :green]
reach = ["Reach A", "Reach B", "Reach C"]
shapes = [:circle, :rect, :star5]
df = DataFrame(x=1:20, y=rand(20), reach=rand(reach,20), stream=rand(stream,20))
dshapes = Dict(reach .=> shapes)
dcolors = Dict(stream .=> colors)
layout = @layout [a b{0.1w}]
p1 = @df df scatter(:x, :y, color=getindex.((dcolors,), :stream), shape=getindex.((dshapes,), :reach), ms=4)
v = [[0] for _ in 1:length(reach)];
p2 = scatter(v, v, lims=(-2,-2), fs=:none, mc=:white, label=label=" " .* permutedims(reach), shape=permutedims(shapes), legend=:top, legendfontsize=3, ms=3, msw=0.5)
v = [[0] for _ in 1:length(stream)];
scatter!(v, v, lims=(-2,-2), framestyle=:none, mc=permutedims(colors), label=" " .* permutedims(stream), ms=4, msw=0, fg_legend = :transparent)
plot(p1, p2, layout=layout)