ST_Ellipse.sql

-- PostGIS functions which will create an pollygon representing an ellipse using the "Pins and String" method
-- 
-- https://en.wikipedia.org/wiki/Ellipse#Pins-and-string_method
--
-- Usage: ST_Ellipse(LINESTRING, DISTANCE, NUMBER OF NODES);
-- Usage: ST_Ellipse(LINESTRING, DISTANCE); -- defaults to 8 nodes
--
-- Args: 
-- LINESTRING -- A two point linestring. Use ST_MakeLine(p1,p2) if you have two points
-- DISTANCE -- The sum of the distance from both foci to a single point on the ellipse
-- NUMBER OF NODES -- How many points to render in each quarter of the ellipse
--
-- Notes: 
--  * If the distance requested is less than the distance between the two nodes, an exception is raised
--  * If the start and end points are the same, then ST_Buffer is called instead
--  * Only tested with UTM-15N data. Will maybe work with anything as long as your geometry and distance use the same units.
--  * YMMV, no guarantees, etc.

CREATE OR REPLACE FUNCTION ST_Ellipse(g geometry, dist float8) RETURNS geometry AS
$$
BEGIN
	RETURN ST_Ellipse(g,dist,8);
END;
$$
LANGUAGE 'plpgsql';


CREATE OR REPLACE FUNCTION ST_Ellipse(g geometry,dist float8,nodes integer) RETURNS geometry AS
$$

DECLARE

-- Some info about our input linestring
srid integer;
startpoint geometry;
endpoint geometry;
radians float8;
length float8;

-- Some info about our centroid
center geometry;
centerX float8;
centerY float8;

-- Temp variables used in the loop as we create points
angle float8;
x float8;
y float8;

-- Some info about our ellipse
xRadius float8;
yRadius float8;

-- Some working variables for data we generate
newpoint geometry;
points geometry[];
poly geometry;

BEGIN
	-- Bail if we're going to have a bad time
	IF ST_Length(g) >= dist THEN
		RAISE EXCEPTION 'The length of this linestring (%) is shorter than or equal to the requested string distance (%)', ST_Length(g), dist;
	END IF;

	IF ST_Length(g) = 0 THEN
		return ST_BUFFER(g,dist,nodes);
	END IF;



	-- Calculate the center of the line string
	center	:= ST_Centroid(g);
	centerX := ST_X(center);
	centerY := ST_Y(center);
	srid	:= ST_Srid(center);
	length  := ST_Length(g);

	-- Calculate the radiuses
	xRadius := ((length/2) + ((dist - length)/2));
	yRadius := sqrt(( dist / 2 )^2 - ( ST_Length(g) / 2 )^2);

	RAISE NOTICE 'For length %, distance % I got x=%,y=%', length, dist, xRadius, yRadius;

	-- And the slope of the string, in radians
	startpoint = ST_StartPoint(g);
	endpoint = ST_EndPoint(g);
	radians := atan( (ST_Y(endpoint) - ST_Y(startpoint)) / (ST_X(endpoint) - ST_X(startpoint)) );


	-- How many nodes? We ask for nodes in a quarter segment to be more similar to st_buffer
	nodes := nodes * 4;

	-- Our point accumulator
	points = '{}'::geometry[];


	FOR node IN 1..nodes LOOP
		angle := node * ( 2 * Pi() / nodes );
		x := xRadius * cos(angle) + centerX;
		y := yRadius * sin(angle) + centerY;
		newpoint = ST_SetSrid(ST_Point(x,y),srid);

		points = array_append(points,newpoint);
	END LOOP;
	
	-- Close the linestring
	points = array_append(points,points[1]);

	
poly = ST_MakePolygon(ST_MakeLine(points));

-- Rotate our ellipse around the centroid
poly = ST_Rotate(poly, radians, center);		

RETURN poly;
END;
$$
LANGUAGE 'plpgsql';