This is a static archive of our old OpenStreetMap Help Site. Please post any new questions and answers at community.osm.org.

st_transform on planet_osm_line - information returned

I'm using this query to get all the lat and lon coordinates on a way from planet_osm_line:

select st_astext(st_transform(way, 4326)) from planet_osm_line where osm_id = 482283890;

The result is: LINESTRING(11.6168586 46.5434487002502,11.6168345 46.5431350002503,11.6167899 46.5429099002503,11.6172158 46.5425886002504,11.6176661 46.5424470002504,11.6180585 46.5424763002504,11.6181788 46.5425033002504,11.618848 46.5427728002503,11.6193672 46.5429585002503,11.6204657 46.5429528002503,11.6214048 46.5432146002503,11.6222983 46.5436187002502,11.6224615 46.5436944002502)

What I see are the lon and lat coordinates for every node constructing the line, but the lat seems to contain extra information (last 6 digits, ie 002502). What is this information and how can I seperate it from the lat coordinate?

st_transform planet_osm_line

asked 07 Feb '20, 12:15

Ruud Brandsma
11●1●1●2
accept rate: 0%

One Answer:

These are floating point represenation errors (from the conversion from lat/lon to mercator coordinates during import, and back in your query); they are not useful data. You can accomplish a rounding with ST_SnapToGrid:

SELECT st_astext(st_snaptogrid(st_transform(way, 4326), .000001))
FROM planet_osm_line 
WHERE osm_id = 482283890;

answered 07 Feb '20, 12:28

Frederik Ramm ♦
82.5k●92●720●1273
accept rate: 23%

Thanks for your reply. Why are these floating point representation errors only occur on the latitude (precision of 13 digits) and not on the longitude (precision always 7 digits)?

(07 Feb '20, 12:39) Ruud Brandsma

This is down to the formula used to compute the spherical mercator projection; computing the latitude value requires trigonometric functions, whereas computing the longitude is just simple multiplication:

double lon2x(double lon) { import std.math; return PI*(lon/180.0) * 6378137.0; }
double lat2y(double lat) { import std.math; return log(tan(PI_4+PI*(lat/360.0)))*6378137.0; }
double x2lon(double x) { import std.math; return 180*(x/6378137)/PI; }
double y2lat(double y) { import std.math; return 180*(2*atan(exp(y/6378137))-PI_2)/PI; }

You could largely avoid this issue if you used the -l (the letter L) flag on osm2pgsql import, this would then refrain from projecting coordinates in the first place and make your st_transform redundant. But it would then require some tweaking to use this database for rendering.

(07 Feb '20, 12:45) Frederik Ramm ♦