How is ST_PointOnSurface calculated?

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?

• 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. Nov 4 '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.
– mjobrien
Nov 4 '13 at 20:16

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.

• Looking at the libgeos code, I believe you are right. The horizontal bisector is found, then the midpoint of the widest intersection is used. Nov 6 '13 at 6:41
• Does this mean that this call always produce the same point always? As in, same if called again with the same geometry? Oct 21 '20 at 11:11
• @iamthadiyan - I would assume so. Jan 9 '21 at 0:30