Friday, December 9, 2011

How-To find the closest point of a MultiLineString from a Point in PostGIS

I had the same trouble before with MultiLineStrings, I realized that when a MultiLinestring can't be merged, all functions like ST_ClosestPoint and ST_Line_Locate_Point doesn't work.(you can know if a MultiLineString can't be merged using the ST_LineMerge function) I've made a pl/pgSQL function based in an old maillist but I added some performance tweaks, It only works with MultiLineStrings and LineStrings (but can be easily modified to work with Polygons). First it checks if the geometry only has 1 dimension, if it has, you can use the old ST_Line_Interpolate_Point and ST_Line_Locate_Point combination, if not, then you have to do the same for each LineString in the MultiLineString. Also I've added a ST_LineMerge for pre 1.5 compatibility :


CREATE OR REPLACE FUNCTION ST_MultiLine_Nearest_Point(amultiline geometry,apoint geometry)
  RETURNS geometry AS
$BODY$
DECLARE
    mindistance float8;
    adistance float8;
    nearestlinestring geometry;
    nearestpoint geometry;
    simplifiedline geometry;
    line geometry;
BEGIN
        simplifiedline:=ST_LineMerge(amultiline);
        IF ST_NumGeometries(simplifiedline) <= 1 THEN
            nearestpoint:=ST_Line_Interpolate_Point(simplifiedline, ST_Line_Locate_Point(simplifiedline,apoint) );
            RETURN nearestpoint;
      END IF;
--      *Change your mindistance according to your projection, it should be stupidly big*
        mindistance := 100000;
        FOR line IN SELECT (ST_Dump(simplifiedline)).geom as geom LOOP
                adistance:=ST_Distance(apoint,line);
            IF adistance < mindistance THEN
                mindistance:=adistance;
                nearestlinestring:=line;
            END IF;
        END LOOP;
        RETURN ST_Line_Interpolate_Point(nearestlinestring,ST_Line_Locate_Point(nearestlinestring,apoint));
    END;
    $BODY$
      LANGUAGE 'plpgsql' IMMUTABLE STRICT;

5 comments:

Nicklas Avén said...

Hmm. St_ClosestPoint should work on multilinestrings, and collections, also nested collections. Do you have an example when it doesn't work?
/Nicklas

Nicklas Avén said...
This comment has been removed by the author.
Nicklas Avén said...
This comment has been removed by the author.
Paco V. said...

yes, I do have an example:

geometry:

MULTILINESTRING((-99.171266262679 19.4251711795705,-99.1712950964062 19.42523852826,-99.1713225890486 19.4252821626056,-99.1713547755568 19.4253207378868,-99.1713795860115 19.4253359150378),(-99.1713795860115 19.4253359150378,-99.1714191485945 19.4253681664926,-99.1714721222225 19.4253909322215,-99.1715277780595 19.4254073741347,-99.1715814222398 19.4254143303281,-99.1716377486291 19.4254105360408,-99.1716913928094 19.4254023150847,-99.1717436958852 19.4253814465015,-99.1717966695132 19.4253542541011),(-99.1679064606144 19.4274357243262,-99.1679845799519 19.4273803916483,-99.1680549879385 19.4273241103911,-99.1681113143278 19.4272545492595,-99.1681522180153 19.4271818262267,-99.1681823928667 19.4270945585441,-99.168191110046 19.4270072908147,-99.1681803812099 19.4269257144165,-99.1681582529855 19.4268757567571,-99.1681448419404 19.4268435056018,-99.168286328466 19.4267859594068,-99.1684941996646 19.4266879411155,-99.1688965310167 19.4264956986176,-99.1693042267868 19.4263021911367,-99.1696783949443 19.4261244928874,-99.1699425925322 19.4259992819046,-99.1700619508333 19.4259417354105,-99.1701504637308 19.4259006307595,-99.1704193551845 19.4257741548446,-99.1704857398576 19.4257412710905,-99.1711220939462 19.4254434198618,-99.1712468166653 19.4254042122283,-99.1713795860115 19.4253359150378))


query:

SELECT ST_Line_Interpolate_Point(ST_LineMerge(the_geom), ST_Line_Locate_Point(ST_LineMerge(the_geom),ST_GeomFromEWKT('srid=4326;POINT(99.170133 19.425992)') ) ) FROM segments WHERE gid =35;

response error:

ERROR: line_locate_point: 1st arg isnt a line

postgis version:

postgis_full_version
----------------------------------------------------------------------------------------
POSTGIS="1.4.0" GEOS="3.1.0-CAPI-1.5.0" PROJ="Rel. 4.7.1, 23 September 2009" USE_STATS

also tried with a higher version, It didn't work either.

Nicklas Avén said...

But that is not ST_ClosestPoint

from your data, this is the result
:

POINT(-99.1679064606144 19.4274357243262)

from this query in PostGIS 1.5 (ST_ClosestPoint was introduced in 1.5)
select st_astext(st_closestpoint(ST_GeomFromEWKT('srid=4326;MULTILINESTRING((-99.171266262679 19.4251711795705,-99.1712950964062 19.42523852826,
-99.1713225890486 19.4252821626056,-99.1713547755568 19.4253207378868,-99.1713795860115 19.4253359150378),(-99.1713795860115 19.4253359150378,
-99.1714191485945 19.4253681664926,-99.1714721222225 19.4253909322215,-99.1715277780595 19.4254073741347,-99.1715814222398 19.4254143303281,
-99.1716377486291 19.4254105360408,-99.1716913928094 19.4254023150847,-99.1717436958852 19.4253814465015,-99.1717966695132 19.4253542541011),
(-99.1679064606144 19.4274357243262,-99.1679845799519 19.4273803916483,-99.1680549879385 19.4273241103911,-99.1681113143278 19.4272545492595,
-99.1681522180153 19.4271818262267,-99.1681823928667 19.4270945585441,-99.168191110046 19.4270072908147,-99.1681803812099 19.4269257144165,
-99.1681582529855 19.4268757567571,-99.1681448419404 19.4268435056018,-99.168286328466 19.4267859594068,-99.1684941996646 19.4266879411155,
-99.1688965310167 19.4264956986176,-99.1693042267868 19.4263021911367,-99.1696783949443 19.4261244928874,-99.1699425925322 19.4259992819046,
-99.1700619508333 19.4259417354105,-99.1701504637308 19.4259006307595,-99.1704193551845 19.4257741548446,-99.1704857398576 19.4257412710905,
-99.1711220939462 19.4254434198618,-99.1712468166653 19.4254042122283,-99.1713795860115 19.4253359150378))'),
ST_GeomFromEWKT('srid=4326;POINT(99.170133 19.425992)')));