User-contributed Functions

This is an area to put utility functions or wrappers around PostGIS.

Construct Ellipse

  x double precision, y double precision, 
  rx double precision, ry double precision, 
  rotation double precision DEFAULT 0.0, 
  quadSeg integer DEFAULT 8)
 RETURNS geometry AS
   SELECT ST_Translate( ST_Rotate( ST_Scale( ST_Buffer(ST_Point(0,0), 0.5, quadSeg), rx, ry), rotation), x, y)
 LANGUAGE 'sql';

Construct an arc given two points on the arc, the centre point and direction or size

 CREATE FUNCTION ST_Arc(startpoint geometry, endpoint geometry, arcenter geometry, direction text) RETURNS geometry
    AS $$
    cwpointonarc geometry;
    ccpointonarc geometry;
    ccarc text;
    cwarc text;
    cwdirection float;
    ccdirection float;
    midpointrads float;
    arcenterutm geometry;
    startpointutm geometry;
    endpointutm geometry;
    thearc geometry;
    majorarc text;
    minorarc text;
        arcenterutm := st_transform(arcenter,utmzone(arcenter));
        startpointutm := st_transform(startpoint,utmzone(arcenter));
        endpointutm := st_transform(endpoint,utmzone(arcenter));

        midpointrads := abs(st_azimuth(arcenterutm,startpointutm) - st_azimuth(arcenterutm,endpointutm));
        IF midpointrads > pi() THEN midpointrads:= (midpointrads - pi())/2; ELSE midpointrads:= midpointrads/2; END IF;
        IF midpointrads > st_azimuth(arcenterutm,startpointutm) THEN midpointrads:= st_azimuth(arcenterutm,startpointutm)/2; END IF;
        IF midpointrads > st_azimuth(arcenterutm,endpointutm) THEN midpointrads:= st_azimuth(arcenterutm,endpointutm)/2; END IF;
        cwdirection := -1*midpointrads;
ccdirection := midpointrads;
cwpointonarc := ST_Translate( ST_Rotate( ST_Translate( startpointutm, -1*ST_X(arcenterutm), -1*ST_Y(arcenterutm)), cwdirection), ST_X(arcenterutm), ST_Y(arcenterutm));
ccpointonarc := ST_Translate( ST_Rotate( ST_Translate( startpointutm, -1*ST_X(arcenterutm), -1*ST_Y(arcenterutm)), ccdirection), ST_X(arcenterutm), ST_Y(arcenterutm));
        cwarc := 'CIRCULARSTRING('||ST_X(startpointutm)||' '||ST_Y(startpointutm)||','||ST_X(cwpointonarc)||' '||ST_Y(cwpointonarc)||','||ST_X(endpointutm)||' '||ST_Y(endpointutm)||')';
ccarc := 'CIRCULARSTRING('||ST_X(startpointutm)||' '||ST_Y(startpointutm)||','||ST_X(ccpointonarc)||' '||ST_Y(ccpointonarc)||','||ST_X(endpointutm)||' '||ST_Y(endpointutm)||')';
IF st_length(st_curvetoline(cwarc)) > st_length(st_curvetoline(ccarc)) THEN majorarc := cwarc; minorarc := ccarc; ELSE majorarc := ccarc; minorarc := cwarc; END IF;
IF direction = 'major' THEN RETURN st_transform(st_setsrid(st_curvetoline(majorarc),utmzone(arcenter)),st_srid(arcenter)); ELSE IF direction = 'minor' THEN RETURN st_transform(st_setsrid(st_curvetoline(minorarc),utmzone(arcenter)),st_srid(arcenter)); END IF; END IF;
        IF direction = 'cw' THEN RETURN st_transform(st_setsrid(st_curvetoline(cwarc),utmzone(arcenter)),st_srid(arcenter)); ELSE IF direction = 'cc' THEN RETURN st_transform(st_setsrid(st_curvetoline(ccarc),utmzone(arcenter)),st_srid(arcenter)); END IF; END IF;

  COMMENT ON FUNCTION st_createarc(geometry,geometry,geometry,text) IS 'Generates an arc based on starting point, ending point, centre of arc, and direction (clockwise: \'cw\', conter-clockwise: \'cc\', or size \'major\' or \'minor\'). All geometries must be of type POINT and this function returns a linestring of the arc based on the original SRID and 32 points per quarter circle. Requires the utmzone function.'; 

Convert Degree Minutes Seconds to Decimal Degrees

This is from Simon Greener's article DMS2DD for PostGIS and is an adaptation of his Oracle function described in Converting Google Earth Formatted Longitude/Latitude points to decimal degrees

Example use:

SELECT round(dms2dd('43° 0''50.60"S'),9) as latitude,
       round(dms2dd('147°12''18.20"E'),9) as longitude;
    RETURNS numeric
       i               numeric;
       intDmsLen       numeric;          -- Length of original string
       strCompassPoint Char(1);
       strNorm         varchar(16) = ''; -- Will contain normalized string
       strDegMinSecB   varchar(100);
       blnGotSeparator integer;          -- Keeps track of separator sequences
       arrDegMinSec    varchar[];        -- TYPE stringarray is table of varchar(2048) ;
       dDeg            numeric := 0;
       dMin            numeric := 0;
       dSec            numeric := 0;
       strChr          Char(1);
       -- Remove leading and trailing spaces
       strDegMinSecB := REPLACE(strDegMinSec,' ','');
       -- assume no leading and trailing spaces?
       intDmsLen := Length(strDegMinSecB);

       blnGotSeparator := 0; -- Not in separator sequence right now

       -- Loop over string, replacing anything that is not a digit or a
       -- decimal separator with
       -- a single blank
       FOR i in 1..intDmsLen LOOP
          -- Get current character
          strChr := SubStr(strDegMinSecB, i, 1);
          -- either add character to normalized string or replace
          -- separator sequence with single blank         
          If strpos('0123456789,.', strChr) > 0 Then
             -- add character but replace comma with point
             If (strChr <> ',') Then
                strNorm := strNorm || strChr;
                strNorm := strNorm || '.';
             End If;
             blnGotSeparator := 0;
          ElsIf strpos('neswNESW',strChr) > 0 Then -- Extract Compass Point if present
            strCompassPoint := strChr;
             -- ensure only one separator is replaced with a blank -
             -- suppress the rest
             If blnGotSeparator = 0 Then
                strNorm := strNorm || ' ';
                blnGotSeparator := 0;
             End If;
          End If;
       End Loop;

       -- Split normalized string into array of max 3 components
       arrDegMinSec := string_to_array(strNorm, ' ');

       --convert specified components to double
       i := array_upper(arrDegMinSec,1);
       If i >= 1 Then
          dDeg := CAST(arrDegMinSec[1] AS numeric);
       End If;
       If i >= 2 Then
          dMin := CAST(arrDegMinSec[2] AS numeric);
       End If;
       If i >= 3 Then
          dSec := CAST(arrDegMinSec[3] AS numeric);
       End If;

       -- convert components to value
       return (CASE WHEN UPPER(strCompassPoint) IN ('S','W') 
                    THEN -1 
                    ELSE 1 
               (dDeg + dMin / 60 + dSec / 3600));

Find Length of Day for a given location and date

CREATE OR REPLACE FUNCTION daylength(geometry, date)
 RETURNS double precision AS 
SELECT 24*(ACOS(1-(1-TAN(radians(Y(transform($1,4326))))*TAN(.409088 * COS(.0172024 * extract(doy from $2)))))/PI());
 LANGUAGE 'sql';

Rotate Geometry around a Point

Superseded by ST_Rotate.

  geom geometry, pt_x double precision, pt_y double precision, angle double precision)
 RETURNS geometry AS 
SELECT ST_Translate( ST_Rotate( ST_Translate( $1, -1*$2, -1*$3), $4), $2, $3)
   LANGUAGE 'sql';

Construct a multilinestring consisting of the interior and exterior rings of a polygon/multipolygon

Superseded by ST_Boundary(geometry)

Example use: - where the_geom is a multipolygon
 SELECT fnpoly_to_rings(the_geom)
 FROM sometable

 CREATE OR REPLACE FUNCTION fnpoly_to_rings(geometry)
 RETURNS geometry AS
 SELECT ST_Collect(the_line)  as multiline
 FROM (SELECT ST_ExteriorRing(the_poly) as the_line
        FROM (SELECT ST_GeometryN($1, g.n) As the_poly
            generate_series(1, ST_NumGeometries($1)) As g(n) ) As polys
        UNION ALL
        SELECT ST_InteriorRingN(the_poly,
        generate_series(1, ST_NumInteriorRings(the_poly))) as the_line
   FROM (SELECT ST_GeometryN($1, g.n) As the_poly
                FROM generate_series(1, ST_NumGeometries($1)) As g(n) ) As polys
        ) As all_lines
COMMENT ON FUNCTION fnpoly_to_rings(geometry) IS 'Takes as argument a multipolygon or polygon and returns a multilinestring consisting of the interior and exterior rings of the polygon/multipolygon';

Offset a straight line

This was written before the existence of ST_OffsetCurve(, which does the same thing and can handle non-straight lines

  CREATE OR REPLACE FUNCTION upgis_lineshift(centerline geometry, dist double precision)
  RETURNS geometry AS
        delx float;
        dely float;
        x0 float;
        y0 float;
        x1 float;
        y1 float;
        az float;
        dir integer;
        line geometry;
        az := ST_Azimuth (ST_StartPoint(centerline), ST_EndPoint(centerline));
        dir := CASE WHEN az < pi() THEN -1 ELSE 1 END;
        delx := ABS(COS(az)) * dist * dir;
        dely := ABS(SIN(az)) * dist * dir;

        IF az > pi()/2 AND az < pi() OR az > 3 * pi()/2 THEN
                line := ST_Translate(centerline, delx, dely) ;
                line := ST_Translate(centerline, -delx, dely);
        END IF;

        RETURN line;
  COMMENT ON FUNCTION upgis_lineshift(geometry, double precision) IS 'Takes a 2D line string and shifts it dist units along
 the perpendicular defined by the straight line between the start and end point
 Convention: (right is positive and left is negative.  right being defined as to right of observer
 standing at start point and looking down the end point)';
Last modified 14 months ago Last modified on Aug 5, 2020, 3:54:59 PM