[postgis-users] Midpoint function?

Rodrigo Martín LÓPEZ GREGORIO rodrigomartin at lopezgregorio.com.ar
Fri Jun 22 06:27:00 PDT 2007


Ok, here is the function. It takes two parameters. The first is a LINESTRING
or a MULTILINESTRING and the second is the location, similar as
line_interpolate_point function. I think it works ok. Obviously if the
MULTILINESTRING have gaps between constituent LINESTRINGs and the location
is right in a "discontinuity" point it will return just one of the point (a
kind of mathematical "left limit" I think although I didn't make several
test about that).
I hope it will be useful for someone. Here goes the code:

-- Function: multiline_interpolate_point(amultils geometry,location float8)

-- DROP FUNCTION multiline_interpolate_point(amultils geometry,location
float8);

CREATE OR REPLACE FUNCTION multiline_interpolate_point(amultils
geometry,location float8)
  RETURNS geometry AS
$BODY$
DECLARE
    accumlength float8;
    totallength float8;
    nearestpoint geometry;
    i integer;

BEGIN
    IF ((location < 0) or (location > 1)) THEN
    RAISE EXCEPTION 'location must be between 0 and 1 --> %', location;
    END IF;
    IF (GeometryType(amultils) != 'LINESTRING') and (GeometryType(amultils)
!= 'MULTILINESTRING') THEN
    RAISE EXCEPTION 'the geometry must be a LINESTRING or a
MULTILINESTRING';
    END IF;

    IF (GeometryType(amultils) = 'LINESTRING') THEN
    nearestpoint:=line_interpolate_point(amultils,location);
    RETURN nearestpoint;
    END IF;
    accumlength := 0;
    totallength := Length(amultils);
    FOR i IN 1 .. NumGeometries(amultils) LOOP
        IF ((accumlength + Length(GeometryN(amultils,i))/totallength) >=
location) THEN
            IF
((location-accumlength)/(Length(GeometryN(amultils,i))/totallength) > 1)
THEN

nearestpoint:=line_interpolate_point(GeometryN(amultils,i),1);
            ELSE

nearestpoint:=line_interpolate_point(GeometryN(amultils,i),(location-accumlength)/(Length(GeometryN(amultils,i))/totallength));
            END IF;
            EXIT;
        ELSE

accumlength:=accumlength+Length(GeometryN(amultils,i))/totallength;
        END IF;
    END LOOP;
    RETURN nearestpoint;
END;
$BODY$
  LANGUAGE 'plpgsql' IMMUTABLE STRICT;
ALTER FUNCTION multiline_interpolate_point(amultils geometry,location
float8) OWNER TO postgres;

Rodrigo.

On 6/22/07, Rodrigo Martín LÓPEZ GREGORIO <
rodrigomartin at lopezgregorio.com.ar> wrote:
>
> I don't know why I mixed my ideas. The function that I send in my previous
> mail find the nearest point of a multilinestring to a given point so it has
> nothing to do with the middle point in the multilinestring. Tomorrow I will
> modify that function and I'll send it again. Sorry about that.
>
> Rodrigo.
>
> On 6/21/07, Rodrigo Martín LÓPEZ GREGORIO <
> rodrigomartin at lopezgregorio.com.ar > wrote:
> >
> > I was testing the LineMerge() solution with different MULTILINESTRINGs
> > and it will only return a LINESTRING geometry if there is no gap between
> > constituent linestrings. I remember then that a few days ago I answered a
> > similar question but that time the linestrings had gaps and I tryed with
> > LineMerge with no luck, so that time I created a function that work with
> > those MULTILINESTRINGs. I know that the LineMerge solution works fine in
> > your case but I send again the code of the function I created. I know maybe
> > there could be a better solution but at least it works :P
> >
> > [...] So I wrote a stored function (preety simple by the way) that get
> > what you need. The function get two parameters... the first is a
> > multilinestring, the second, your point. The function loops over all
> > linestrings in the multilinestring, get the nearest one and then with that
> > linestring calculates the location of the nearest point of the linestring to
> > your point. To get this function working run this script and it will create
> > your function (take a look at the end of the function and chanche the user
> > name if necesary; mine was postgres):
> >
> > -- Function: multiline_locate_point(amultils geometry,apoint geometry)
> >
> > -- DROP FUNCTION multiline_locate_point(amultils geometry,apoint
> > geometry);
> >
> > CREATE OR REPLACE FUNCTION multiline_locate_point(amultils
> > geometry,apoint geometry)
> >   RETURNS geometry AS
> > $BODY$
> > DECLARE
> >     mindistance float8;
> >     nearestlinestring geometry;
> >     nearestpoint geometry;
> >     i integer;
> >
> > BEGIN
> >     mindistance := (distance(apoint,amultils)+100);
> >     FOR i IN 1 .. NumGeometries(amultils) LOOP
> >         if distance(apoint,GeometryN(amultils,i)) < mindistance THEN
> >             mindistance:=distance(apoint,GeometryN(amultils,i));
> >             nearestlinestring:=GeometryN(amultils,i);
> >         END IF;
> >     END LOOP;
> >
> > nearestpoint:=line_interpolate_point(nearestlinestring,line_locate_point(nearestlinestring,apoint));
> >     RETURN nearestpoint;
> > END;
> > $BODY$
> >   LANGUAGE 'plpgsql' IMMUTABLE STRICT;
> > ALTER FUNCTION multiline_locate_point(amultils geometry,apoint geometry)
> > OWNER TO postgres;
> >
> > Once the function is stored and available for use, you can make the
> > query like:
> >
> > SELECT *, multiline_locate_point(the_geom,PointFromText('POINT(517651
> > 2121421)', 42102))
> > FROM mainroad ORDER BY Distance(the_geom,PointFromText('POINT(517651
> > 2121421)', 42102)) LIMIT 1
> >
> > However, the query speed is a little bit slow if you do the query in
> > that way. You can use the next query instead wich first gets the nearest
> > multilinestring to your point and then call my function with that geometry
> > instead of calling the function for all geometries:
> >
> > SELECT *, multiline_locate_point(the
> > _geom,PointFromText('POINT(517651 2121421)', 42102)) FROM
> > (SELECT *, Distance(the_geom,PointFromText('POINT(517651 2121421)',
> > 42102)) as dist
> > FROM mainroad ORDER BY Distance(the_geom,PointFromText('POINT(517651
> > 2121421)', 42102)) LIMIT 1 ) as foo
> >
> > On my PC the second query takes 1/6 of the time that tooks the first one
> > so the speed is much better at least here ;)
> >
> > Rodrigo
> >
> > On 6/21/07, Pradeep B V <pradeepbv at gmail.com> wrote:
> > >
> > >  On 6/21/07, Rodrigo Martín LÓPEZ GREGORIO <rodrigomartin at lopezgregorio.com.ar
> > > > wrote:
> > > >
> > > > As Nicolas says you should use:
> > > >
> > > > select line_interpolate_point(LineMerge(the_geom),0.5).
> > > >
> > > > you can try this and see if the point is the one you want:
> > >
> > >
> > > It works fine.
> > >
> > > and here is the proof.
> > >
> > > select distance(startpoint(the_geom),
> > > line_interpolate_point(LineMerge(the_geom),0.5)) as stdist,
> > > distance(endpoint(the_geom), line_interpolate_point(LineMerge(the_geom),
> > > 0.5)) as endist from testdataset limit 1;
> > >
> > >                     stdist         |       enddist
> > >               ----------------------+----------------------
> > >  0.000737003466840363 | 0.000737003466840363
> > >
> > > Thanks Rodrigo/Nicolas.
> > >
> > >  - Pradeep B V
> > > www.btis.in
> > >
> > > _______________________________________________
> > > postgis-users mailing list
> > > postgis-users at postgis.refractions.net
> > > http://postgis.refractions.net/mailman/listinfo/postgis-users
> > >
> > >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20070622/bc5dc2af/attachment.html>


More information about the postgis-users mailing list