User-contributed Functions

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

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 $$
  DECLARE
    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;
  BEGIN
        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;
  END;

  $$
  LANGUAGE 'plpgsql' IMMUTABLE;
  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 lon/lat decimals

This is from Simon Greener's article ​DMS2DD for PostGIS

and was 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;

CREATE OR REPLACE FUNCTION DMS2DD(strDegMinSec varchar)
    RETURNS numeric
    AS
    $$
    DECLARE
       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);
    BEGIN
       -- 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;
             Else
                strNorm := strNorm || '.';
             End If;
             blnGotSeparator := 0;
          ElsIf strpos('neswNESW',strChr) > 0 Then -- Extract Compass Point if present
            strCompassPoint := strChr;
          Else
             -- 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 
                END 
               *
               (dDeg + dMin / 60 + dSec / 3600));
    End 
$$
    LANGUAGE 'plpgsql' IMMUTABLE;

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(geom,x,y)

CREATE OR REPLACE FUNCTION RotateAtPoint(the_geom geometry, pt_x double precision, pt_y double precision, rotrads double precision)
 RETURNS geometry AS 
$$
SELECT ST_Translate( ST_Rotate( ST_Translate( $1, -1*$2, -1*$3), $4), $2, $3)
$$
   LANGUAGE 'sql';

Construct Ellipse ( x,y, rx,ry [,rotation [,#of segments in ¼ of ellipse] ] )

CREATE OR REPLACE FUNCTION ST_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 a multilinestring consisting of the interior and exterior rings of a polygon/multipolygon

This already exists in PostGIS: see 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
         FROM
            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
 $$
  LANGUAGE 'sql' IMMUTABLE;
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(​http://postgis.net/docs/ST_OffsetCurve.html), 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
  $$
  DECLARE
        delx float;
        dely float;
        x0 float;
        y0 float;
        x1 float;
        y1 float;
        az float;
        dir integer;
        line geometry;
  BEGIN
        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) ;
        ELSE
                line := ST_Translate(centerline, -delx, dely);
        END IF;

        RETURN line;
  END;
  $$
  LANGUAGE 'plpgsql' IMMUTABLE;
  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)';