How to find the boundary of a region from a dataset

The following function extractregion might be something to start with. First, it smoothes along the third entries in the data set by TotalVariationFilter. TTotalVariationFilter has the advantage over, e.g. GaussianFilter, that it tends to preserve edge jumps, which is really desired here. Afterwards, we interpolate the smoothed data set, plot a ContourPlot, extract the contour lines, and convert them to a single MeshRegion.

extractregion[data_, level_, smoothing_] := 
 Module[{a, n, g, contour},
  a = data;
  n = FirstPosition[Unitize[Subtract[a[[1 ;; -2, 1]], a[[2 ;; -1, 1]]]], 1][[1]];
  a[[All, 3]] = Flatten[TotalVariationFilter[Partition[a[[All, 3]], n], smoothing]];
  g = Interpolation[a];
  contour = 
   ContourPlot[
    g[x, y] == level, {x, Sequence @@ g[[1, 1]]}, {y, Sequence @@ g[[1, 2]]}];
  RegionUnion@Map[
    gc \[Function] MeshRegion[gc[[1]], Cases[gc, _Line, All]],
    Cases[contour, _GraphicsComplex, All]
    ]
  ]

This is how it behaves on the provide example data:

Show[
 ListDensityPlot[data, InterpolationOrder -> 0, PlotRange -> All, 
  AspectRatio -> 4/5],
 extractregion[data, 0.5, 3],
 extractregion[data, -0.5, 1.2]
 ]

It doesn't close the quarter of the disk in the lower left corner. I don't now an elegant way to do that at the moment.

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.