[postgis-users] Calculating center point for geographic polygon that crosses IDL

Hazuka, Brad bhazuka at usgs.gov
Tue Sep 1 08:55:55 PDT 2015


Hello,

I am fairly new to PostGIS and I am hoping some one can provide me a hand
with an issue relating to scenes crossing the International Date Line (IDL)
and finding the center point.

I currently maintain a PostgreSQL database that contains footprints for
satellite scenes. The scenes have global coverage and range in size from
~1500Km squared down to ~100Km squared. The footprints are stored in a
geography column as basic polygons as seen below.

SRID=4326;POLYGON((179.255436389612 49.0944921724897,-179.572711844818
49.0656217040858,-179.619180111394 48.3778311277411,179.224815891643
48.4060150740941,179.255436389612 49.0944921724897))

I need to be able to determine the center point of each scene in order to
provide some additional information to various applications using the data.
Based on forums and various other web resources, it appears the way to
calculate the center point is by converting geography to geometry and using
the ST_CENTROID function as seen below.

ST_ASEWKT(ST_CENTROID(ST_GEOMFROMEWKT(ST_ASEWKT(footprint))))

Using my polygon example from above, the function call would look like the
following.

ST_ASEWKT(ST_CENTROID(ST_GEOMFROMEWKT('SRID=4326;POLYGON((179.255436389612
49.0944921724897,-179.572711844818 49.0656217040858,-179.619180111394
48.3778311277411,179.224815891643 48.4060150740941,179.255436389612
49.0944921724897))')))

Here is where the problem lies. In the above example, the centroid returns
a center point value of the following:

 st_asewkt
 ----------------------------------------------------
 SRID=4326;POINT(-0.148022213748242 48.7359898569253)

Note that the centroid longitude is -0.148022213748242. The footprint I
provided crosses the IDL. The scene center longitude should be ~179.8181.

I was able to hack together an inelegant solution that will correct for the
IDL on a Cartesian system up to a point. I have provided the code below.
The code appears to work but will breakdown if the scene footprint is to
large and I am sure I have not accounted for everything.

CREATE FUNCTION test (in i_wkt varchar)
RETURNS TABLE
( centerlat FLOAT
, centerlon FLOAT
);
AS
$$
DECLARE
    l_ctrlat    FLOAT;
    l_ctrlon    FLOAT;
    l_lonmin    FLOAT;
    l_lonmax    FLOAT;
    l_lontmp    FLOAT;

BEGIN
    SELECT ST_Y(ST_CENTROID(ST_GEOMFROMEWKT(i_wkt)))
         , ST_X(ST_CENTROID(ST_GEOMFROMEWKT(i_wkt)))
         , ST_XMIN(ST_GEOMFROMEWKT(i_wkt))
         , ST_XMAX(ST_GEOMFROMEWKT(i_wkt))
      INTO l_ctrlat
         , l_ctrlon
         , l_lonmin
         , l_lonmax;

    -- See if the scene crossesthe IDL if so correct the negative value
    IF (abs(l_lonmin) + abs(l_lonmax)) > 180 THEN
        l_lontmp := 360 + l_lonmin;
        l_ctrlon := ((l_lonmax + l_lontmp)/2);
    END IF;

    -- Determine which side of the IDL the longitude belongs on
    IF l_ctrlon > 180 THEN
        l_ctrlon := l_lonctr - 360;
    END IF;

    RETURN QUERY
    SELECT l_ctrlat AS centerlat
         , l_ctrlon AS centerlon;

END;
$$
LANGUAGE plpgsql;


Is there a way to determine the center point for a geographic scene that
crosses the IDL using current PostGIS functions that will remain true for
any foot print size that does not require additional code? If so, what
combination of functions do I need to use?

Do I need to re-project the scene prior to determining the center point
with ST_CENTROID since I move the footprint to a Cartesian geometry type?

If no functions are available, is there a preferred method for determining
center point of a polygon that crosses the IDL that will better allow for
larger sizes of polygons in a geographic system?

I am currently running PostgreSQL 9.3.4 with PostGIS extension 2.1.0.

Any help you could provide would be appreciated.




Bradley J. Hazuka
Database Administrator
Innovate!, Inc.
Contractor to the U.S. Geological Survey (USGS)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20150901/6b2c9543/attachment.html>


More information about the postgis-users mailing list