PostGIS Ellipse Issue

Question

I am trying to create ellipses from fields in my posgres table. I found this custom Ellipse function and have been trying to tweak it to meet my needs.

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'; 

My inputs into this are a lat and lon. Smajor/Sminor measurements in meters. I receive my rotation in degrees but have a written a function that has converted them to radians per the requirements as stated in some of the other forums I've seen with this similar problem. I have tweaked the above method to this:

CREATE OR REPLACE FUNCTION Ellipse2(double precision, double precision, double precision, double precision, double precision, integer) 
RETURNS geometry AS 
$$ 
SELECT ST_Translate( ST_SetSRID( ST_Rotate( ST_Scale( ST_Buffer( ST_SetSRID( ST_Point(0,0), 3395), 0.5, $6), $3, $4), $5), 4326), $1, $2) 
$$ 
LANGUAGE 'sql'; 

I create the point (0,0) and then set it's SRID to 3395(world meters based system) in hopes that the follow on functions use meters when it extrapolates the smajor/sminor sides. I'm getting ellipses that are massive though (i.e. continent size not 600 meter size) I'm pretty sure my issue is with projection based units. Can someone point out to me what I am doing wrong?

Answer 1

The SRID has been set already on your point.

You need to transform the SRID, not set it, after you rotate it:

SELECT ST_Translate( ST_Transform( ST_Rotate( ST_Scale( ST_Buffer( ST_SetSRID( ST_Point(0,0), 3395), 0.5, $6), $3, $4), $5), 4326), $1, $2) 

Answer 2

The function above doesn't work perfectly because the scale and projection may distort the ellipse.

Normally the coordinates are managed with WGS84 that is not a metric system so you have to convert the points in a metric system and then translate it again.

First of all, you have to choose a metric system for your region. Here in Italy the best map system is RDN2008 that is a metric system.

So we start with:

CREATE OR REPLACE FUNCTION wgs84toRdn2008(double precision, double precision)
RETURNS geometry AS 
$$ 
SELECT ST_Transform(ST_setsrid(ST_MakePoint( $1, $2),4326),6875)
$$ 
LANGUAGE 'sql'; 

CREATE OR REPLACE FUNCTION Rdn2008toWGS84(double precision, double precision)
RETURNS geometry AS 
$$ 
SELECT ST_Transform(ST_setsrid(ST_MakePoint( $1, $2),6875),4326)
$$ 
LANGUAGE 'sql';

For a better use in the big Store Procedure, I have created this two function to extract only the X or Y.

CREATE OR REPLACE FUNCTION wgs84ToRdn2008X(double precision, double precision)
RETURNS double precision AS 
$$ 
SELECT ST_X(wgs84toRdn2008($1,$2))
$$ 
LANGUAGE 'sql';

CREATE OR REPLACE FUNCTION wgs84ToRdn2008Y(double precision, double precision)
RETURNS double precision AS 
$$ 
SELECT ST_Y(wgs84toRdn2008($1,$2))
$$ 
LANGUAGE 'sql';

Now, the final EllipseFunction:

CREATE OR REPLACE FUNCTION WGS84Ellipse(x double precision, y double precision, max double precision, min double precision, rotationAngle integer) 
RETURNS geometry AS 
$BODY$
DECLARE
  rdn2008dX double precision;
  rdn2008dY double precision;
  rotationA1 double precision;
  rotationA2 double precision;
  rotationB1 double precision;
  rotationB2 double precision;
  newX double precision;
  newY double precision;
  calcAngle integer;  
  tmpPoint geometry;
  ellipsePoints geometry;  
BEGIN
    rdn2008dX := wgs84ToRdn2008X(x,y);      
    rdn2008dY := wgs84ToRdn2008Y(x,y);
    calcAngle := 0; 

    rotationAngle := 360 - rotationAngle;

    create temporary table geo_temp_ellipse ("geometry" geometry not null) on commit drop;

    LOOP 
        EXIT WHEN calcAngle > 360 ; 
        -- Calculating rotation factor --
        rotationA1 := cos(radians(calcAngle)) * cos(radians(rotationAngle));
        rotationA2 := sin(radians(calcAngle)) * sin(radians(rotationAngle));    

        rotationB1 := cos(radians(calcAngle)) * sin(radians(rotationAngle));        
        rotationB2 := sin(radians(calcAngle)) * cos(radians(rotationAngle));        
        -- Calculating x and y on ellipse arc
        newX := ((max * rotationA1) - (min * rotationA2) + rdn2008dX);      
        newY := ((max * rotationB1) + (min * rotationB2) + rdn2008dY);              
        tmpPoint := rdn2008toWGS84(newX,newY);

        INSERT INTO geo_temp_ellipse VALUES (tmpPoint);                 

        calcAngle := calcAngle + 10 ; 
    END LOOP ; 

     ellipsePoints := ST_ConvexHull(ST_Collect(ARRAY(SELECT geometry from geo_temp_ellipse)));      
     return ellipsePoints;  
END; 
$BODY$
LANGUAGE PLPGSQL;

The parameters are: x, y => center of ellipse in WGS84 max, min => small and big diameters of ellipse rotation => rotation (degrees) from the vertical

Hope this procedure may help someone!