Additional User-contributed Functions

This is an area to put utility functions or wrappers around PostGIS Distance /Spatial Reference support functions PL/PGSQL/SQL Functions

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

 CREATE FUNCTION st_createarc(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(double precision, double precision, double precision, double precision, double precision, integer)
 RETURNS geometry AS
 $$
   SELECT ST_Translate( ST_Rotate( ST_Scale( ST_Buffer(ST_Point(0,0), 0.5, $6), $3, $4), $5), $1, $2)
$$
 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)';