= Additional User-contributed Functions =
 

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

 * '''Generate an arc given two points on the arc, the centre point and direction'''

{{{
#!sql
 CREATE OR REPLACE FUNCTION st_createarc(startpoint geometry, endpoint geometry, arcenter geometry,direction text)
  RETURNS geometry AS
  $$
  DECLARE
    pointonarc geometry;
    thearc text;
    thedirection float;
    midpointrads float;
    originalsrid int;
    arcenterutm geometry;
    startpointutm geometry;
    endpointutm geometry;
  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 direction = 'cw' THEN thedirection := -1*midpointrads; ELSE thedirection := midpointrads; END IF;
        pointonarc := ST_Translate( ST_Rotate( ST_Translate( startpointutm, -1*ST_X(arcenterutm), -1*ST_Y(arcenterutm)), thedirection), ST_X(arcenterutm), ST_Y(arcenterutm));
        thearc := 'CIRCULARSTRING('||ST_X(startpointutm)||' '||ST_Y(startpointutm)||','||ST_X(pointonarc)||' '||ST_Y(pointonarc)||','||ST_X(endpointutm)||' '||ST_Y(endpointutm)||')';
        RETURN st_transform(st_setsrid(st_curvetoline(thearc),utmzone(arcenter)),st_srid(arcenter));
  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\'). 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.'; 





}}}

 * '''Find UTM (WGS84)  SRID for a point (in any SRID)'''

{{{
#!sql
 -- Function: utmzone(geometry)
 -- DROP FUNCTION utmzone(geometry);

 CREATE OR REPLACE FUNCTION utmzone(geometry)
   RETURNS integer AS
 $BODY$
 DECLARE
     geomgeog geometry;
     zone int;
     pref int;

 BEGIN
     geomgeog:= ST_Transform($1,4326);

     IF (ST_Y(geomgeog))>0 THEN
 	pref:=32600;
     ELSE
	pref:=32700;
     END IF;

     zone:=floor((ST_X(geomgeog)+180)/6)+1;

     RETURN zone+pref;
 END;
 $BODY$ LANGUAGE 'plpgsql' IMMUTABLE
   COST 100;
}}}

 * '''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(geom,x,y)'''
{{{
#!sql
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';
}}}

 * '''Create Ellipse (x,y,rx,ry,rotation,#of segments in 1/4 of ellipse)'''
{{{
#!sql
CREATE OR REPLACE FUNCTION 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';
}}}

 * '''Return a multilinestring consisting of the interior and exterior rings of a polygon/multipolygon'''

''This already exists in PostGIS: see 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';

}}}

 * '''Shift a straight line along dist units along its perpendicular'''
{{{
#!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)';
}}}