Not sure why its going wrong.
Somehow removing the only two points outside the region changes things
f = Function[{x, y, z},
N[-2 (1 + Sqrt[1 - y]) + (3 + 2 Sqrt[1 - y]) y] < x <
N[-2 + 2 Sqrt[1 - y] + 3 y - 2 Sqrt[1 - y] y]];
datas = Select[data, f @@ # &];
Select[data, ! (f @@ #) &]
Grid@{ListContourPlot[#,
BoundaryStyle -> Directive[Red, Thick],
RegionFunction -> f] & /@ {data, datas}}
But adding the option MaxPlotPoints
imposes a regular grid on your irregular data which seems to solve your issue
Grid@{
{ListContourPlot[data, RegionFunction -> f,
BoundaryStyle -> Directive[Red, Thick]],
ListContourPlot[data, RegionFunction -> f,
BoundaryStyle -> Directive[Red, Thick], MaxPlotPoints -> 40]
},
{ListContourPlot[data, RegionFunction -> f, Mesh -> All],
ListContourPlot[data, RegionFunction -> f, Mesh -> All,
MaxPlotPoints -> 40]
}}