@natemcintosh, some possible issues in your two attempts. For the first: too few input points for bi-cubic (default) spline computations and possibly an inadequate error tolerance for the linear input data; for the second attempt: the sorted data defines a linear feature in the 2D plane, probably a ill-conditioned problem for 2D splines (plus the error tolerance issue).
Produced a working example below (ps: used surface for a better display)
using Dierckx, Plots
x = 2*rand(100); y = rand(100);
z = (x .- 1).^2 .+ y.^2
spl = Spline2D(x, y, z; kx=3, ky=3, s=1e-4)
xg = LinRange(0,2,200); yg = LinRange(0,1,200);
zspl = evalgrid(spl, xg, yg)
surface(xg, yg, zspl'; legend=:false, xl="x", yl="y", zl="z")
scatter!(x, y, z; ms=4, mc=:red, ma=0.6, title="Dierckx Spline2D",titlefont=12)
# xmesh = ones(200)'.*xg
# ymesh = yg'.*ones(200)
# scatter(xmesh, ymesh, zspl; ms=0, mc=:cyan, legend=:false)
