# Find ordinal polygon points

How might I go about finding the following points of a `MULTIPOLYGON` in PostGIS?

• Northeasternmost
• Northwesternmost
• Southeasternmost
• Southwesternmost

Can't seem to find a `ST_` function to accomplish this.

• I've 2 questions. Are you searching for an existing point(=vertice) of the polygon, or can it be a new one (like in my query)? Do you want to know the northeasternmost point within a range (that represents NE by 22.5 to 67.5 degrees)? – Stefan Dec 16 '13 at 11:26
• Thanks @StefanB. 1. It doesn't have to be an existing point. It can be on a line between points. 2. Yes, this sounds fine. To give some context. I want to add an 'x' button to the top-right of the poly, mimicking a window, as a close button. – Jordan Arseno Dec 18 '13 at 4:10
• I edited my answer one more time. The problem with the geometries is solved and you can specify the value that stands for a distinct direction (22.5 for NE, 111.5 for SE, ...) – Stefan Dec 19 '13 at 15:04

The function below can find the points that are in the northeasternmost, southeasternpost, southwesternmost, northwesternmost position.

Additional you have to create the TYPE `ordinal`

``````CREATE TYPE ordinal AS (ID smallint, geom geometry(point, <yourEPSGnumber>));
``````

You have to replace `<yourEPSGnumber>` with your EPSG defined in your geometry Table (e.g. 4326,...).

So here is my approach:

1. Calculating the centroid of a polygon
2. Creating 2 lines that span the geodirection (22.5 to 67.5)
3. Splitting the polygon with the merged lines
4. Querying the polygon that represents the NE,SE,SW,NW quadrants
5. Calculating the most distant point from the quadrant to the centroid

! You have to replace `ID` with an ID, or name that identifies your polygon (same in the `TYPE` defined above). You have to replace `YourPolygonTable` with the name of your table, where your polygons are stored.

``````CREATE OR REPLACE FUNCTION ordinal(ID integer, quadrant double precision) RETURNS SETOF ordinal AS \$\$
WITH centroid AS
(SELECT
YourPolygonTable.ID,
(ST_Centroid(YourPolygonTable.geom)) AS vertex
FROM YourPolygonTable
WHERE ID=\$1),
newline AS
(SELECT
ST_SetSRID(ST_Translate(
ST_Rotate(
ST_MakeLine(
ST_MakePoint(0.0,500.0),
ST_MakePoint(0.0,0.0)),
ST_X(vertex), ST_Y(vertex)),
ST_SRID(vertex)) AS geom
FROM centroid, generate_series(0,45,45) AS s(a)
),
span AS
(SELECT
centroid.ID,
ST_LineMerge(ST_Union(newline.geom)) AS geom
FROM newline, centroid
GROUP BY ID),
multiobject AS
(SELECT
span.ID,
ST_Split(YourPolygonTable.geom,span.geom) AS geom,
generate_series(1,100) AS n
FROM span, YourPolygonTable
WHERE YourPolygonTable.ID=\$1),
objects AS
(SELECT
n,
ST_GeometryN(multiobject.geom,n) AS geom
FROM multiobject
WHERE n <= ST_NumGeometries(multiobject.geom)),
(SELECT
multiobject.ID,
objects.geom AS geom
FROM objects, multiobject
WHERE multiobject.n <= ST_NumGeometries(multiobject.geom)
ORDER BY LEAST(ST_Area(ST_GeometryN(objects.geom,multiobject.n)))
LIMIT 1),
points AS
(SELECT
AND YourPolygonTable.ID=\$1)

SELECT DISTINCT ON
(ST_Distance(centroid.vertex,points.geom))
points.ID,
points.geom
FROM centroid, points
ORDER BY ST_Distance(centroid.vertex,points.geom) DESC
LIMIT 1
\$\$ LANGUAGE 'sql';
``````

One problem is that the ST_GeometryN(geom, n) varies from to polygon to polygon. Setting n = 1 in the "object" subquery doesn't select the smallest part for every polygon.

In the subquery "newline" the one point is 500m away from the vertex (centroid). When you have very large polygons, you have to change this value.

You can query NE data with:

``````SELECT (ordinal(4,22.5)).*
``````

You can query SE data with:

``````SELECT (ordinal(4,111.5)).*
``````

...

The first parameter (\$1) is your polygon ID, the second one (\$2, see "newline" subquery) the start value for the quadrants. When you want to select the southeasternmost point you have to set the paramter to 111.5, and so on! The clause generate_series(0,45,45) in the subquery "newline" counts the 22.5 degrees with 45 degrees up (one time).

Displaying the results in QGIS with the database manager (for NE): The green point is result of the query from geogeek, the purple one is the result of the function above!

Maybe we have to discuss the approach with the centroid.

EDIT

I solved the problem with geometries (n; which is the smaller one...) by using the `LEAST` function. The subquery "quadrant" results the smallest polygon from the split.

• This seems impressive, though, I'm not a PostGIS pro and could use some help calling this. When I execute this, I get `ERROR: relation "accumulation" does not exist`. Must I find and replace all occurrences of "accumulation" with the name of my table? How can we make this table independent? – Jordan Arseno Dec 18 '13 at 4:12
• Replace with the name of my geometry* rather. – Jordan Arseno Dec 18 '13 at 4:21
• I've edited my answer. You have to replace "yourtable" with the name of your table, where you have stored your polygon(s). The same with the partnumber! Consider that you can't see the quadrant like in the picture. – Stefan Dec 18 '13 at 11:31
• Can you share the source for the polygon you used? – Jakub Kania Dec 19 '13 at 13:39
• How can I share this with you? This is a big polygon with many vertices. Besides, I edited my answer... – Stefan Dec 19 '13 at 15:07

You can make a PostGis query to get feature points then select the nearest point to northeastern corner of the extent, using the query below:

select ST_AsText(dp) , st_distance(dp , (select ST_POINT( ST_MinX(geom ) , ST_MaxY(geom ) , ST_SRID(geom)) from geo_table where id = 1) ) as dis from (SELECT ST_DumpPoints(geom) AS dp from geo_table where id = 1 ) As foo

order by dis
limit 1;

ST_MinX(geom ) , ST_MaxY(geom ) is the Northeasternmost of the extent of the feature, for the Northwesternmost you will use ST_MaxX(geom ) , ST_MaxY(geom ), so the other corners will follow the same logic.

• I tried the function ST_MinX(geom) and ST_MaxY(geom) because I was curios about. But my postgis didn't know this functions, seemingly they don't exist. – Stefan Dec 16 '13 at 10:08
• try `ST_XMax` and `ST_YMin`, i'm using an older postgis version. – geogeek Dec 16 '13 at 10:11

I don't have the privilege to comment yet.

A quick heads up. Depending on the end goal, geogeeks answer can fail if the geometry crosses the 180th meridian. Looking at the globe along the equator plane, with north axis pointing upwards, the westmost point is not the same as the "leftmost" point. In normal cases, the leftmost point will have the least X-coordinate (i.e. westmost), while the rightmost point will have the greatest X-coordinate (i.e. eastmost). On the other hand, for a geometry crossing the 180th meridian, the leftmost point will have a greater X-coordinate than the rightmost point. This is due to the X-coordinate being defined within [-180, 180]. This can be solved by converting to all-positive X-coordinates by adding 360 degrees to the negative X-values.