= 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(startpoint)); endpointutm := st_transform(endpoint,utmzone(endpoint)); midpointrads := abs(st_azimuth(arcenterutm,startpointutm) - st_azimuth(arcenterutm,endpointutm))/2; 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)'; }}}