5

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.

  • 1
    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
4

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)),
                radians(($2+s.a)*-1)),
          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)),
quadrant AS
(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
    quadrant.ID,
    (ST_DumpPoints(quadrant.geom)).geom
FROM quadrant, YourPolygonTable
    WHERE ST_Intersects(YourPolygonTable.geom,quadrant.geom)
        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):

enter image description here

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
3

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.

  • 1
    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
  • 2
    try ST_XMax and ST_YMin, i'm using an older postgis version. – geogeek Dec 16 '13 at 10:11
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.