You can (1) use ListDensityPlot with options InterpolationOrder -> 1 and MeshFunctions to create mesh lines, and (2) use the function smoothCP from this answer to smooth the mesh lines:
ClearAll[smoothCP]
smoothCP = Module[{pr = PlotRange @ #, nF = Nearest[Join @@ Cases[Normal @ #,
Line[a_] :> Transpose[#2 /@ Transpose@a], ∞]]},
# /. GraphicsComplex[a_, b__] :> GraphicsComplex[
If[Or @@ MemberQ @@@ Thread @ {pr, #}, #, First @ nF @ #] & /@ a, b] ] &;
ldp1 = ListDensityPlot[data, InterpolationOrder -> 0,
PlotRange -> All, AspectRatio -> 4/5, ImageSize -> 300];
ldp2 = ListDensityPlot[data, InterpolationOrder -> 1,
PlotRange -> All, AspectRatio -> 4/5, ImageSize -> 300, PlotStyle -> None,
MeshFunctions -> {Round@#3 &}, Mesh -> 2, MeshStyle -> Directive[Thick, Red]];
Row[{Show[ldp1, ldp2], Show[ldp1, smoothCP[ldp2, GaussianFilter[#, {5, 5}] &],
PlotRange -> All]}]
Is it possible to extract the contours as separate lists?
To extract the smoothed lines, you can use
lines = Cases[Normal[smoothCP[ldp2, GaussianFilter[#, {5, 5}] &]], _Line, Infinity]
Short /@ lines
{Line[{{0., 0.387454}, << 58 >>, {0.392532, 0.}}],
Line[{{0.504121, 0.483339}, << 146 >>, {0.504121, 0.483339}}]}
To get just the coordinates of the two lines, use Line[x_, ___]:>x instead of _Line above.