Issues with ST_INTERSECTS for node geometry in way geometry

Question

Hello

I am trying to work with the full history file. I've used Mazderminds scripts to load an excerpt into a Postgres database.

Some spatial queries are giving me the results I expect, however I've found that none of nodes intersect the ways.

For example using way id : 3565059, version 5 and node id 17543678, version 5 MUST intersect. However ST_INTERSECTS returns FALSE.

On the same query using ST_DWITHIN returns TRUE using parameter of 0.1 metres, but FALSE with any greater precision.

Has anyone else encountered the same issue or can shed some light on this problem ?

Answer 1

I believe the problem lies in the fact the point geometry in the database created by MaZdermind's scripts is rounded to 6 decimal places, while the vertices in the line geometry is not rounded. The below query on a full history extract on nepal shows this. In the results ST_INTERSECTS will return True on record 2. The remainder returns False, but returns True using a max resolution of 0.000001 m and ST_DWITHIN.

The actual precision of the node_point geometry on the OSM database can be manually verified through the URL(openstreetmap.org/node/...(way id)../history)and paging down to the corresponding version number. Equally the point can be confirmed and being part of a way by the same method.

 --Create a function to find the closest vertex in a linestring to a given point (http://gis.stackexchange.com/questions/111711/find-the-closest-line-vertex-from-the-linestring-using-postgis)

CREATE OR REPLACE FUNCTION ST_AsMultiPoint(geometry) RETURNS geometry AS  'SELECT ST_Union((d).geom) FROM ST_DumpPoints($1) AS d;'LANGUAGE sql IMMUTABLE STRICT COST 10;

--Select a subset of data from a full history node table

WITH pointsub as (Select * FROM hist_point LIMIT 600) ,

--Find records in the full history way table (hist_line) where the nodes subset is 0.000001 metres from the way (hist_line). Add condition on time stamp creation to make sure you have found points in the lines. Display the vertex in the linestring which is closest to that point, compare against point geom

results1 as (SELECT hist_line.id line_id,pointsub.id point_id,pointsub.version point_version,hist_line.version,ST_AsText(ST_ClosestPoint(ST_AsMultiPoint(hist_line.geom), pointsub.geom)) closest_vertex_in_line,ST_AsText(pointsub.geom) point_geom
FROM hist_line,pointsub WHERE ST_DWITHIN(pointsub.geom, hist_line.geom,0.000001)AND hist_line.valid_from = pointsub.valid_from),

--As there are many versions of the lines. select 1 instance of the line where the geometry relates to each unique point

results2 as (SELECT *,ROW_NUMBER() OVER(PARTITION by results1.point_id) as rn FROM results1)

SELECT results2.point_id,results2.point_version,results2.line_id,
results2.line_version,results2.point_geom,results2.closest_vertex_in_line 
FROM results2 WHERE rn = 1

alt text

Answer 2

1

I would regard this as much more likely to be a question about PostGIS processing. There are numerous examples of problems with queries of this kind on GIS Stackinfo & Stackinfo. Various generic techniques can be used to mitigate against such issue: such as forcing valdity st_makevalid and checkgeometry, checking validity st_isvalid, pre-snapping objects with st_snap.

Also remember that there is no guarantee that any given polygon OSM way or relation at a particular time will make always be well formed for PostGIS. I wrote something about this a while back.

More detail of what you are doing may help as well. I presume these queries are against hist_view_line etc. as described MaZdermind's write-up on github.

permanent link

answered 11 Nov '15, 10:37

SK53 ♦
28.1k●48●268●433
accept rate: 22%

Thanks for the response. As I believe the problem is related to the underlying data or the transformation I hoped to have more success on this forum than stackoverflow. But I do appreciate this forum is weighted towards contributions, rather than consumption.

The queries are from hist_point against hist_line. I am doing some analysis on rollbacks and looking to see which rolled back nodes haved effected which ways. Spatial queries within the tables work perfectly, the issue arises on queries between tables.

(12 Nov '15, 09:27) Sam75

1

Well the other place is against MaZdermind's repo. This is a fairly involved technical area, and relatively few people have a setup which can test the issue. Certainly, I'd only really start to worry if it wasnt resolved with the techniques I suggest.

(12 Nov '15, 12:06) SK53 ♦

Answer 3

0

I suspect line 106 here https://github.com/MaZderMind/osm-history-renderer/blob/master/importer/handler.hpp is the culprit. This sets the total number of decimal places used and setting that to 8 would cause exactly the results you see above.

Removing "<< std::setprecision(8)" and reimporting should fix the issue.

permanent link

answered 03 Dec '15, 16:38

SimonPoole ♦
44.7k●13●326●701
accept rate: 18%

Many thanks for looking at this Simon. Unfortunately I had already tried your suggestion, << std::setprecision(8)" can also be found in line 260 of handler file and line 48 of the geombuilder file. Changing the precision or removing any of these lines, does not change the output format, the point geom remains rounded to precision 8 (or 6 decimal places - I have not imported files that with geometries above 100 degrees longitude). The linestring vertices precision are as per the url of the node. Any advice from Jochen Topf?

(03 Dec '15, 21:09) Sam75

Did you reimport? The rounding is done when writing to the database (I would suggest testing on a small extract naturally). Jochen is not aware of any rounding in OSMIUM (but the tools in question are using a very old version in any case).

(03 Dec '15, 21:11) SimonPoole ♦

Correct, amended script(s) and dropped existing tables and reimported. For good measure I set up a fresh database. I am aware it uses the older version of OSMIUM. There were two warnings when unpacking osmium regarding the xml tag on the doxyfile, I was given the option to update or remove, I chose to update. Since Madzermind's application is based on the original OSMIUM I assume it is best not to play with that. When you fixed the bug in OSM2PGSQL was it as simple as resetting a precision parameter ?

(03 Dec '15, 22:13) Sam75

Do you think it could be something in the WKB writer : http://geos.osgeo.org/doxygen/WKBWriter_8h_source.html and its point class dependency http://geos.osgeo.org/doxygen/Point_8h_source.html.

(03 Dec '15, 23:48) Sam75

Looking at line 117 it writes out WKT and IMHO that should be correct if you remove the setPrecision, what -dows- look wrong there is the SRID spec, given that you should be writing lat lon, but maybe I missed something.

(04 Dec '15, 13:33) SimonPoole ♦

I don't have a reason or good enough understanding to disagree with your rationale, but removing lines of code with setprecision definitely does not work for me. Nonetheless, many thanks for your efforts, it's much appreciated.

(04 Dec '15, 20:56) Sam75

Not fully understanding your problem, but the default precision in C++ is (often? always?) 6. So removing std::setprecision() doesn't necessarily lead to a higher precision. Instead it might be necessary to increase the specified precision or even add std::setprecision() at some places. So in this specific case std::setprecision(8) was not meant to lower the precision but to raise it instead.

(04 Dec '15, 22:18) scai ♦

Thanks, I have also tried raising the precision, and tried std ::setprecision() just now, again no change. This line of code can't be causing the bug.

(04 Dec '15, 22:33) Sam75

1

IMHO one us will actually have to install the software and test .. will however take a couple of days.

(05 Dec '15, 00:18) SimonPoole ♦

3

setprecision(8) is not enough. You need 7 digits after the decimal point to get the full precision OSM supports, plus 3 digits before it, makes 10. So setprecision(10) should be correct. Removing the setprecision() will get you the default of 6 which is definitely not enough.

(27 Dec '15, 16:28) Jochen Topf

1

Many thanks Simon, Alex and Jochen. Increasing setprecision solves the problem. I had failed with earlier attempts as I had not remade the executable file after amending the 'setprecision' in the handler and geombuilder files. As Jochen says setprecision(10) allows for the full OSM precision.

(28 Dec '15, 22:41) Sam75

showing 5 of 11 show 6 more comments

Issues with ST_INTERSECTS for node geometry in way geometry

Follow this question

Related questions