OpenStreetMap
Small fix to osm2pgsql commited
Posted by SimonPoole on 9 December 2012 in English. Last updated on 10 December 2012.One property that differentiates the noaddress layer on qa.poole.ch from other similar layers is that it checks for address information on nodes inside and on the building polygons. This is trivial, ST_Covers for example should return true if a node is in or on the polygon, but, a tiny bit suprising, didn’t work with a standard (srid 900913) osm2psql imported database.
For example
select St_Covers(p.way,n.way) from planet_osm_polygon p, planet_osm_point n where p.osm_id=26681434 and n.osm_id=1946037917;
returns false even though node 1946037917 is one of the nodes used to construct the building polygon 26681434.
On inspection of the data in the database the issue was glaringly obvious: the nodes in the planet_osm_point table had a lot more decimal places than the corresponding coordinates in the polygons. This is simply due to the polygons and ways being constructed from the node cache that scales the coordnate values so that they fit in a 32bit integer. In the process the decimal places get truncated (note that there is no loss of information since the main database schema stores the coordinates scaled as 32bit integers anyway, the extra decimal places are likely a result of reprojection for the import).
The fix is easy: simply scale the node coordinates the same way (and back) before inserting them into the database.
If you have a database imported in the standard spherical mercartor projection you can truncate the existing coordinates in planet_osm_point with the following SQL and don’t need to reimport (you’ll need to get the patched osm2pgsql version from the OSM svn too naturally):
update planet_osm_point set way = ST_SetSRID(ST_MakePoint(trunc(cast(St_X(way) as numeric),2),trunc(cast(St_Y(way) as numeric),2)),900913) ;
Discussion
Comment from FraMauro on 9 December 2012 at 22:21
First of all I would like to thank you for qa.poole.ch: it’s really quite useful! I already noticed the fact that all the housenumbers “on” the building were rightly identified.
I just wanted to signal another possibility for the housenumber localization; maybe you’ll find it interesting. Sometimes the housenumber is on the fence that surrounds the building.
Take for example this building: way 33381972 It is a school. Apparently no housenumber.
Yet the school is surrounded by a fence: way 39419103
The node 443720945 is part of the fence and shows the housenumber.
Thanks again for your map!
Comment from pnorman on 1 April 2013 at 21:43
I ran into a similar problem with addressmerge, except my addresses could be marginally outside the buildings. I solved it by adding a column in which I buffered the geometry with
UPDATE import_addresses SET buffered_geom = geometry(ST_Buffer(geography(geom),N));where N was the distance, on the order of 1-2 meters (real meters, not mercator meters). For some purposes you might want to reproject but I was working in WGS84.
For mercator you could just use ST_Buffer and multiply N by a factor, but I had to do the reprojections.
Another trick for doing the JOINs between addresses and buildings was to pre-filter with && and then do ST_DWithin(geography(…), geography(…), N)
In my code this was expressed as
a LEFT JOIN b ON a.id != b.id AND ST_Expand(a.geom, PRECOMPUTED_DEGREES) && b.geom AND ST_DWithin(geography(a.geom), geography(b.geom), DISTANCE)with PRECOMPUTED_DEGREES being DISTANCE * scale, where scale is chosen from the point closest to the pole in the region I’m considering.