33

The PostGIS documentation states that ST_PointOnSurface returns "a POINT guaranteed to lie on the surface". It seems like this function could be trivially implemented to give results that satisfy the documentation but provide little real-world utility, though I'm certain that PostGIS provides a non-trivial implementation.

This introduction to PostGIS provides a nice comparison and contrast of ST_Centroid with ST_PointOnSurface and says that "[ST_PointOnSurface] is substantially more computationally expensive than the centroid operation".

Is there a more thorough explanation of how ST_PointOnSurface is calculated? I've been using ST_Centroid, but have encountered some edge cases in my data where the centroid is outside the geometry. I believe that ST_PointOnSurface is the correct substitute, but the function name and documentation leave room for uncertainty.

Further, is the computational expense of ST_PointOnSurface incurred even if the centroid does lie within the geometry already?

mjobrien
  • 433
  • 4
  • 5
  • It exists precisely since the centroid of non-convex polygons is not always included in it. It has nothing to do with heights and DEMs if that's the confusing part of the name. Implementation details are best checked in the code, but I believe you'd get a better answer on GIS.se. – lynxlynxlynx Nov 04 '13 at 18:16
  • Good point on GIS.se. Is there a way to migrate this question there? I understand why both functions exist. I find the name confusing because there are infinitely many points on the surface of the polygon goemetries I'm working with. However, only a small subset of those points serves my purpose. I want to know I'm getting a point that makes sense for how I want to use it. –  Nov 04 '13 at 20:16

2 Answers2

41

Based on a few experiments, I think ST_PointOnSurface() works roughly like this, if the geometry is a polygon:

  1. Trace an east-west ray, lying half-way between the northern and southern extents of the polygon.
  2. Find the longest segment of the ray that intersects the polygon.
  3. Return the point that is half-way along said segment.

That may not make sense, so here's a sketch of a polygon with a ray dividing it into a northern and southern parts:

             _
            / \             <-- northern extent
           /   \
          /     \
         /       \
        /         \      __
       /           \    /  \
      /_ _ _ P _ _ _\  / _ _\  P = point-on-surface
     /               \/      \
    /                         \
   /            C              \   C = centroid
  /                             \
 /                              /
/______________________________/  <-- southern extent

Thus, ST_PointOnSurface() and ST_Centroid() are usually different points, even on convex polygons.

The only reason for the "surface" in the name, I think, is that if the geometry has 3D lines then the result will be simply be one of the vertexes.

I would agree that more explanation (and better naming) would have been useful and hope a GEOS programmer might shed some more light on the matter.

Martin F
  • 8,948
  • 36
  • 58
-1

The Algorithm mentioned by Martin F is not guaranteed to work for every polygon. I am not sure if this is the algorithm used by said function or not, however consider the example below: enter image description here

In this example, the longest line segment that the bisector crosses is between the points (-0.147,2.5) and (3.207,2.5), and the midpoint of this segment is the point (1.53,2.5), which is certainly not on the surface of the polygon.

Edit: I do realise this should've been a reply to Mark F's answer, however I didn't know how to include or reference an image in a reply.

anis messaoudi
  • 1
  • 1
  • 3
    As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. – Community Dec 07 '22 at 10:00
  • 2
    Note that the algorithm, as i describe it, has step 2: "Find the longest segment of the ray that intersects the polygon." [emphasis added]. On your example, that would be on the right-hand portion of your U-shaped polygon. – Martin F Aug 22 '23 at 21:02