| Version 26 (modified by , 4 years ago) ( diff ) |
|---|
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
Generate 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 long 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';
- Create Ellipse (x,y,rx,ry,rotation,#of segments in 1/4 of ellipse)
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)
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
This was written before the existence of ST_OffSetCurve. Using ST_OffsetCurve(http://postgis.net/docs/ST_OffsetCurve.html) does the same 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)';
Note:
See TracWiki
for help on using the wiki.
