Nearest linestring given a point, a tolerance and a set of linestrings?

Discussion:

Stefan Keller

2013-01-23 02:24:43 UTC

What is the nearest interpolated point on a linestring given a single
point, a tolerance (say 100m) and a set of linestrings (i.e. a table
called myways with linestring geometry column); return way/linestring
id and distance?

It's a use case similar to snapping or to a navigation system where
you have a point and a set of lines.

I found hints to ST_Line_Interpolate_Point and ST_Line_Locate_Point
but theses posts/blogs are rather outdated and I'm not sure if KNN
index would do the job better?

- S.

M***@dlr.de

2013-01-23 08:06:13 UTC

Permalink

Hey! What about ST_ClosestPoint (http://postgis.refractions.net/documentation/manual-2.0/ST_ClosestPoint.html)?

Greets, Michi

-----Ursprüngliche Nachricht-----
Gesendet: Mittwoch, 23. Januar 2013 03:25
An: PostGIS Users Discussion
Betreff: [postgis-users] Nearest linestring given a point, a tolerance and a set
of linestrings?
What is the nearest interpolated point on a linestring given a single point, a
tolerance (say 100m) and a set of linestrings (i.e. a table called myways with
linestring geometry column); return way/linestring id and distance?
It's a use case similar to snapping or to a navigation system where you have a
point and a set of lines.
I found hints to ST_Line_Interpolate_Point and ST_Line_Locate_Point but
theses posts/blogs are rather outdated and I'm not sure if KNN index would
do the job better?
- S.
_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Stefan Keller

2013-01-23 11:04:34 UTC

Permalink

Hi Michi,

Post by M***@dlr.de
Hey! What about ST_ClosestPoint (http://postgis.refractions.net/documentation/manual-2.0/ST_ClosestPoint.html)?
Greets, Michi

That's it. but I got unexpected malfunctioning of the query which uses
ST_Distance:

-- Solution with KNN order index => OK:
SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom)) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
WHERE tags @> hstore('railway', 'rail')
ORDER BY ST_ClosestPoint(osm_line.way, mypt.geom) <-> mypt.geom
LIMIT 1

-- Solution with ST_Distance => wrong result. linestring is far away :-<
SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom)) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
WHERE tags @> hstore('railway', 'rail')
ORDER BY ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom))
LIMIT 1

You can reproduce that in http://labs.geometa.info/postgisterminal .
Simply paste it into the Query textarea.
Here's 'mypoint'

SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(8.8211
47.2221)', 4326), 900913)) as geom, 'XXX' AS label
LIMIT 1

-S.

Post by M***@dlr.de

_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

=?utf-8?Q?Nicklas_Av=E9n?=

2013-01-23 11:59:43 UTC

Permalink

Hallo

Your problem is:
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom))

there is two things with this, on warning and one error.

The error is that it should be
ST_Distance( mypt.geom,ST_ClosestPoint(osm_line.way, mypt.geom))
to get what you want. The way you wrote it you ask for thedistance from the interpolated point on the line to the line which will be 0 in all cases.
The notice or warning is that you should instead write
ST_Distance( osm_line.way, mypt.geom)

It is the same thing and will be twice as fast since it is only doing the calculation once instead of doing the exactly same thing twice.

ST_ClosestPoint is just one of the points that ST_Distance uses to get its result. To get both points use ST_ShortestLine.

So, from your query it doesn't look like you want any interpolated opint at all, just the distance (wich is calculated from a interpolated point behind the scenes)

Then you can use ST_DWithin to speed up since ST_Dwithin can use spatial indexes to do a first bounding box search.

Second and very important
Using KNN distance as you do will not always give you the right answer since the <-> operator only works on the center of the bounding box, not the geoemetry itself. So another geometry with the bbox-center much more far away can come much closer than the one you have found. In this case it happened to be the same and give correct result.

Your query rewritten would be something like:

SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way, mypt.geom) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
WHERE tags @> hstore('railway', 'rail') AND ST_DWithin(osm_line.way, mypt.geom,100)
ORDER BY ST_Distance(osm_line.way, mypt.geom)
LIMIT 1;

A slightly other approach I wrote about long ago here:
http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

HTH

Nicklas Avén

2013-01-23 Stefan Keller wrote:

Hi Michi,

Post by Stefan Keller

Post by M***@dlr.de
Hey! What about ST_ClosestPoint (http://postgis.refractions.net/documentation/manual-2.0/ST_ClosestPoint.html)?
Greets, Michi

That's it. but I got unexpected malfunctioning of the query which uses
SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom)) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
ORDER BY ST_ClosestPoint(osm_line.way, mypt.geom) <-> mypt.geom
LIMIT 1
-- Solution with ST_Distance => wrong result. linestring is far away :-<
SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom)) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
ORDER BY ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom))
LIMIT 1
You can reproduce that in http://labs.geometa.info/postgisterminal .
Simply paste it into the Query textarea.
Here's 'mypoint'
SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(8.8211
47.2221)', 4326), 900913)) as geom, 'XXX' AS label
LIMIT 1
-S.

Post by M***@dlr.de

-----UrsprÃŒngliche Nachricht-----
Gesendet: Mittwoch, 23. Januar 2013 03:25
An: PostGIS Users Discussion
Betreff: [postgis-users] Nearest linestring given a point, a tolerance and a set
of linestrings?
What is the nearest interpolated point on a linestring given a single point, a
tolerance (say 100m) and a set of linestrings (i.e. a table called myways with
linestring geometry column); return way/linestring id and distance?
It's a use case similar to snapping or to a navigation system where you have a
point and a set of lines.
I found hints to ST_Line_Interpolate_Point and ST_Line_Locate_Point but
theses posts/blogs are rather outdated and I'm not sure if KNN index would
do the job better?
- S.
_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Stefan Keller

2013-01-23 14:10:52 UTC

Permalink

Hi Nicklas

Post by Stefan Keller
Hallo
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom))

(...)

Post by Stefan Keller
So, from your query it doesn't look like you want any interpolated opint at
all, just the distance (wich is calculated from a interpolated point behind
the scenes)

Exactly.

(...)

Post by Stefan Keller
Second and very important
Using KNN distance as you do will not always give you the right answer since
the <-> operator only works on the center of the bounding box, not the
geoemetry itself. So another geometry with the bbox-center much more far
away can come much closer than the one you have found. In this case it
happened to be the same and give correct result.

Yes, but I think one could use <#> instead of <-> (if newest release
is available):

SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom)) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)', 4326),
900913) AS geom) AS mypt
WHERE tags @> hstore('railway', 'rail')
ORDER BY ST_ClosestPoint(osm_line.way, mypt.geom) <#> mypt.geom
LIMIT 1
Many thanks. That seems to be the final solution

SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way, mypt.geom) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)', 4326),
900913) AS geom) AS mypt
WHERE tags @> hstore('railway', 'rail')
AND ST_DWithin(osm_line.way, mypt.geom, 100) -- tolerance=100m
ORDER BY 2
LIMIT 1

Post by Stefan Keller
http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

Nice blog.
But I think that asks either for the closest point itself, whereas I
ultimately only want the closest line (or id) and the distance.
Thus, consequently it leads to the same solution, as above?

-S.

Post by Stefan Keller
Hallo
ST_Distance(osm_line.way,ST_ClosestPoint(osm_line.way, mypt.geom))
there is two things with this, on warning and one error.
The error is that it should be
ST_Distance( mypt.geom,ST_ClosestPoint(osm_line.way, mypt.geom))
to get what you want. The way you wrote it you ask for thedistance from the
interpolated point on the line to the line which will be 0 in all cases.
The notice or warning is that you should instead write
ST_Distance( osm_line.way, mypt.geom)
It is the same thing and will be twice as fast since it is only doing the
calculation once instead of doing the exactly same thing twice.
ST_ClosestPoint is just one of the points that ST_Distance uses to get its
result. To get both points use ST_ShortestLine.
So, from your query it doesn't look like you want any interpolated opint at
all, just the distance (wich is calculated from a interpolated point behind
the scenes)
Then you can use ST_DWithin to speed up since ST_Dwithin can use spatial
indexes to do a first bounding box search.
Second and very important
Using KNN distance as you do will not always give you the right answer since
the <-> operator only works on the center of the bounding box, not the
geoemetry itself. So another geometry with the bbox-center much more far
away can come much closer than the one you have found. In this case it
happened to be the same and give correct result.
SELECT
ST_AsText(osm_line.way) AS geom,
ST_Distance(osm_line.way, mypt.geom) AS label
FROM osm_line,
(SELECT ST_Transform(ST_GeomFromText('POINT(8.8211 47.2221)',
4326), 900913) AS geom) AS mypt
ORDER BY ST_Distance(osm_line.way, mypt.geom)
LIMIT 1;
http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/
HTH
Nicklas Avén
Hi Michi,

Post by Stefan Keller

Post by M***@dlr.de
Hey! What about ST_ClosestPoint
(http://postgis.refractions.net/documentation/manual-2.0/ST_ClosestPoint.html)?
Greets, Michi

Post by M***@dlr.de

_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

f***@gmx.de

2013-01-23 14:33:05 UTC

Permalink

Hi Postgis User-list,

i've downloaded and migrated some srtm3 rasterdata into my PostgreSQL/Postgis 2.0 Database successfully. Now i would like to calculate the mean elevation value of a given polygon area by intersecting my polygon layer with the raster grid.

Using a SQL-Statement, which was shown by Pierre Racine at the FOSS4G 2011:

SELECT polyID,
(ST_Intersection(geom, rast)).geom poly,
(ST_Intersection(geom, rast)).val elevation
FROM polygon, srtm WHERE ST_Intersects(geom, rast);

What i expected the polygon-result geometry would look like:
_ _ _ _
|_|_|_|_|
|_|_|_|_|
|_|_|_|_|
|_|_|_|_|

What the result actually looks like:
_ _ _ _
|___|_ |
|_ |_| |
|_|___|_|
|_|_|_|_|

The raster cells with the same elevation value got merged, why?
The mean value which based on the actual result cannot be correct. Same result different SQL-Statement without selecting the geometry:

Select polyid,(stats).mean from (Select a.polyid,st_summarystats(st_clip(b.rast,1,a.geom))as stats from polygon a, srtm b where ST_intersects(a.geom,b.rast )

I'm very grateful for any help!

Best regards,
Simon

Bborie Park

2013-01-23 16:05:48 UTC

Permalink

By the looks of the docs for ST_Intersection, what you're getting is
correct. The second paragraph indicates that the raster gets
vectorized using ST_DumpAsPolygon() which does merge neighboring cells
with the same value in one larger geometry.

http://postgis.net/docs/manual-2.0/RT_ST_Intersection.html

The query using ST_Clip is more appropriate.

-bborie

Post by f***@gmx.de
Hi Postgis User-list,
i've downloaded and migrated some srtm3 rasterdata into my PostgreSQL/Postgis 2.0 Database successfully. Now i would like to calculate the mean elevation value of a given polygon area by intersecting my polygon layer with the raster grid.
SELECT polyID,
(ST_Intersection(geom, rast)).geom poly,
(ST_Intersection(geom, rast)).val elevation
FROM polygon, srtm WHERE ST_Intersects(geom, rast);
_ _ _ _
|_|_|_|_|
|_|_|_|_|
|_|_|_|_|
|_|_|_|_|
_ _ _ _
|___|_ |
|_ |_| |
|_|___|_|
|_|_|_|_|
The raster cells with the same elevation value got merged, why?
Select polyid,(stats).mean from (Select a.polyid,st_summarystats(st_clip(b.rast,1,a.geom))as stats from polygon a, srtm b where ST_intersects(a.geom,b.rast )
I'm very grateful for any help!
Best regards,
Simon
_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Simon Appelt

2013-01-23 20:08:43 UTC

Permalink

Hi Bborie,

thanks for your fast answer! Wow, i really missed to read that :-(
But my second query ,which i posted below, provides the same Mean-Value
result as the first. If i understand you correct, the st_clip command
should end up with a different result .

Well i'll do some further testing tomorrow and keep the list up to date.

Simon

Post by Bborie Park
By the looks of the docs for ST_Intersection, what you're getting is
correct. The second paragraph indicates that the raster gets
vectorized using ST_DumpAsPolygon() which does merge neighboring cells
with the same value in one larger geometry.
http://postgis.net/docs/manual-2.0/RT_ST_Intersection.html
The query using ST_Clip is more appropriate.
-bborie

_______________________________________________
postgis-users mailing list
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Pierre Racine

2013-01-23 19:39:47 UTC

Permalink

Post by f***@gmx.de
The raster cells with the same elevation value got merged, why?

Because intersecting with every single pixel outline would be just too long.

Post by f***@gmx.de
The mean value which based on the actual result cannot be correct.

You actually want to compute a weighted mean, not simply a mean. The weighted mean give a correct result (actually more accurate than the clip option but requiring more geometric processing).

Have a look at the tutorial which is now a little bit old:

http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01

The other way is as Bborie said. ST_Clip() the raster and ST_SummaryStat().

Pierre

=?utf-8?Q?Nicklas_Av=E9n?=

2013-01-23 14:14:58 UTC

Permalink

Post by =?utf-8?Q?Nicklas_Av=E9n?=
Second and very important
Using KNN distance as you do will not always give you the right answer since
the <-> operator only works on the center of the bounding box, not the
geoemetry itself. So another geometry with the bbox-center much more far
away can come much closer than the one you have found. In this case it
happened to be the same and give correct result.

Yes, but I think one could use <#> instead of <-> (if newest release

No, <#> calculates against the whole bounding box, still not the real geometry.

/Nicklas