I would do it like this: first use a Capture Attribute node set to Vector and Face to get the normals of the faces.
Plug the Attribute output into a Vector Math node set to Dot Product, as the second vector you use 0/0/1, that is the pointing upwards in Z direction which would be perfectly horizontal.
Now the dot product returns 1, when both vectors point in the same direction (in this case upwards) and lower values when they differ from each other. So now with a Math node set to Less Than you can set a threshold how much deviation from the perfectly flat you will accept. Values less than the threshold will result in a Boolean value of 1, which you plug into a Delete Geometry node. This helps you see which areas now remain as flat.
Now some of the remaining mesh islands might be flat, but too small and you do not want to distribute points there. So you take Face Area node and a Mesh Island node's Island Index output and plug them into an Accumulate Field node set to Face. Then plug the Total output into another Less Than node and the result in another Delete Geometry node, to delete mesh islands that are smaller than the threshold you've set in the Less Than.
The remaining mesh you can now plug into a Distribute Points on Faces node. The resulting points you combine with the original landscape mesh through a Join Geometry node and that's it. You can later decide how flat it should be or how large the islands.
By the way, the face area threshold cannot really distinguish if an area is large enough in length and width, so a very long thin strip of flat terrain will be possible, too.