= User-contributed Functions = This is an area to put utility functions or wrappers around PostGIS. == Construct Ellipse == {{{ #!sql 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 an arc given two points on the arc, the centre point and direction or size == {{{ #!sql 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 Decimal Degrees == This is from Simon Greener's article [http://www.spatialdbadvisor.com/postgis_tips_tricks/115/dms2dd-for-postgis/ DMS2DD for PostGIS] and is an adaptation of his Oracle function described in [http://www.spatialdbadvisor.com/oracle_spatial_tips_tricks/114/converting-google-earth-formatted-longitude-latitude-points-to-decimal-degrees 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 == {{{ #!sql 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 [https://postgis.net/docs/manual-3.0/ST_Rotate.html ST_Rotate].'' {{{ #!sql CREATE OR REPLACE FUNCTION RotateAtPoint( 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 [https://postgis.net/docs/manual-3.0/ST_Boundary.html ST_Boundary(geometry)]'' {{{ #!sql 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 {{{ #!sql 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)'; }}}