[postgis-users] ST_LineToCurve for less than 4 points

[Simon Greener]

Finding the centre and radius of a circle (for projected data) can be done as follows:

CREATE OR REPLACE FUNCTION Find_Circle(p_pt1 geometry, p_pt2 geometry, p_pt3 geometry)
   RETURNS geometry AS
$BODY$
DECLARE
    v_Centre geometry;
    v_radius Numeric;
    v_CX     Numeric;
    v_CY     Numeric;
    v_dA     Numeric;
    v_dB     Numeric;
    v_dC     Numeric;
    v_dD     Numeric;
    v_dE     Numeric;
    v_dF     Numeric;
    v_dG     Numeric;
BEGIN
    IF ( p_pt1 is null OR p_pt2 is null OR p_pt3 is null ) Then
       RAISE EXCEPTION 'All supplied points must be not null.';
       Return Null;
    End If;
    If ( ST_GeometryType(p_pt1) <> 'ST_Point' OR
         ST_GeometryType(p_pt1) <> 'ST_Point' OR
         ST_GeometryType(p_pt1) <> 'ST_Point' ) Then
       RAISE EXCEPTION 'All supplied geometries must be points.';
       Return Null;
    End If;
    v_dA := ST_X(p_pt2) - ST_X(p_pt1);
    v_dB := ST_Y(p_pt2) - ST_Y(p_pt1);
    v_dC := ST_X(p_pt3) - ST_X(p_pt1);
    v_dD := ST_Y(p_pt3) - ST_Y(p_pt1);
    v_dE := v_dA * (ST_X(p_pt1) + ST_X(p_pt2)) + v_dB * (ST_Y(p_pt1) + ST_Y(p_pt2));
    v_dF := v_dC * (ST_X(p_pt1) + ST_X(p_pt3)) + v_dD * (ST_Y(p_pt1) + ST_Y(p_pt3));
    v_dG := 2.0  * (v_dA * (ST_Y(p_pt3) - ST_Y(p_pt2)) - v_dB * (ST_X(p_pt3) - ST_X(p_pt2)));
    -- If v_dG is zero then the three points are collinear and no finite-radius
    -- circle through them exists.
    If ( v_dG = 0 ) Then
       RETURN NULL;
    Else
       v_CX := (v_dD * v_dE - v_dB * v_dF) / v_dG;
       v_CY := (v_dA * v_dF - v_dC * v_dE) / v_dG;
       v_Radius := sqrt(power(ST_X(p_pt1) - v_CX,2) + power(ST_Y(p_pt1) - v_CY,2) );
    End If;
    RETURN ST_SetSRID(ST_MakePoint(v_CX, v_CY, v_radius),ST_Srid(p_pt1));
END;
$BODY$
   LANGUAGE plpgsql VOLATILE STRICT;

select ST_AsText(Finv_dCircle(ST_MakePoint(525000, 5202000),
                             ST_MakePoint(525500, 5202000),
                             ST_MakePoint(525500, 5202500)));

select ST_AsEWKT(Finv_dCircle(ST_SetSrid(ST_MakePoint(525000, 5202000),32615),
                             ST_SetSrid(ST_MakePoint(525500, 5202000),32615),
                             ST_SetSrid(ST_MakePoint(525500, 5202500),32615)));

I have a whole bunch of COGO functions for PostGIS that I have never posted on my website.
Perhaps it is time to do this.

regards
Simon
On Wed, 29 Feb 2012 01:06:58 +1100, Stephen Woodbridge <woodbri at swoodbridge.com> wrote:

> On 2/28/2012 3:38 AM, Sandro Santilli wrote:
>> On Tue, Feb 28, 2012 at 09:36:34AM +0100, Denis Rouzaud wrote:
>>> Hi,
>>>
>>> Can someone tells me why ST_LineToCurve needs at least a linestring
>>> of 4 points to be run?
>>> If I have three points, I should be able to create a circularstring, right?
>>>
>>> Here the error I have:
>>> SELECT ST_LineToCurve(ST_GeomFromText('LINESTRING(554803.177682475
>>> 145390.853708235,554796.469135702 145401.404724093,554769.634948609
>>> 145404.042478058)'));
>>> ERROR:  pta_desegmentize needs implementation for npoints<  4
>>
>> If I read the message correctly it answers exactly your questions.
>> It "needs implementation". A patch is welcome.
>
> Is this the algorithm that is needed to solve this?
> http://paulbourke.net/geometry/circlefrom3/
>
> -Steve
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>


-- 
Holder of "2011 Oracle Spatial Excellence Award for Education and Research."
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, FME, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
   Email: simon at spatialdbadvisor.com
   Voice: +61 362 396397
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
GeoHash: r22em9r98wg
NAC:W80CK 7SWP3