OSM and speed bump map in navigators

Since starting from the beginning is not interesting, I will start from the end.
We did it the same. We got speed bumps from OpenStreetMap data, crossed them with the terrible commercial program of Navitel, made a web viewer of these speed bumps, and an interface for adding them to newbies at http://latlon.org/tc/ . And they even wrote a small press release , a link to which can be sent to friends and acquaintances, motorists.
But for Habr, you can tell something special: how it all works inside, and how it was done.

The main rule of OpenStreetMap is Have fun. Everything described was done not for the sake of money or some other goodies, and not even for your own convenience.

At one of the Minsk mini-OSMovok (mini, because nowhere, except for the IRC channel, they are not announced) I was awarded with the task: “Do you want to do a useful thing? get from osma lezhekopy in navitel. ” Then everything was somehow forgotten, and after a couple of months I still remembered, thanks to the heroic work of comrade stas_symba from the forum w3bsit3-dns.com.ru. He has been collecting and posting updates on speed bumps around Belarus for several months alone. “But it’s easier to collect such data immediately in OSM,” I thought, and I sat down at the export.

Format

')
Since I have neither a navigator, nor Navitel, nor a car, I had to do everything blindly, reading the forums and thinking with brains. For a start and an understanding of the scale of the disaster - in short about the fact that we managed to get about the format of data from the forums.

The file cops.txt (and / or speedcam.txt) is the usual csv, in the first line is the header, in all subsequent lines the data. Title:

idx,x,y,type,speed,dirtype,direction
Where:

• idx - unique object number;
• x, y - coordinates, in the projection longitude-latitude EPSG: 4326;
• speed - the allowed speed for the object;
• dirtype - type of direction of action - in 1 or 2 sides;
• direction - action direction;

Tool selection

Tools for solving the problem are usually chosen from what is at hand. On hand there were a lot of scripts for threading OSM XML , and the PostGIS base. Despite the fact that it’s more familiar to work with streaming scripts, the last field of the speedcam file hinted that you would have to work with the geometry of objects, which is what PostGIS is famous for.

To start, the usual Mapnik database was imported on the home server with osm2pgsql . The first request was written quickly and simply:
select * from planet_osm_point where traffic_calming is not NULL;
In response to this, I returned a lot of fields, including geometry in WKB format, encoded in hex. Not good at all.

I had to get into the maual of PostGIS, in which the functions ST_X and ST_Y were found. It seemed that is necessary. Rewrite:
gis=> select st_x(way), st_y(way) from planet_osm_point where traffic_calming is not NULL limit 1;
st_x | st_y
------------------+------------------
3085590.21426068 | 7159526.18035388
(1 )


Suddenly? It is expected, but it is necessary to make an explanation.

Two projections are widely used in web-cartography: for transmission and display to the user - latitude-longitude, it’s EPSG: 4326 , and “Google’s”, it’s also “Mercator on a sphere”, it’s EPSG: 3857, it’s also EPSG: 900913, it is EPSG: 3785 (why are there so many codes? a long history of disputes between large corporations, amateurs and registrars, worthy of a separate post). The projection is good because the transition from it to 4326 in any language takes up two lines of mathematics on the strength . It is she who is metric, and not in angular coordinates, used when displaying maps on the screen. And to facilitate calculations, it is stored in the database Mapnik.

Well okay. We re-project.
gis=> select ST_X(transform(p.way,94326)) as X, ST_Y(transform(p.way,94326)) as Y from planet_osm_point p where traffic_calming is not NULL limit 1;
x | y
------------+------------
27.7183285 | 53.9438334
(1 )

What is next? Speed.

Nobody in OSM puts speed on speed bumps. We need to figure out where to get it. We ask the first person that came across (it turned out to be wildMan ), we get the answer - according to the traffic rules, the allowed speed is 60 in the city, 90 outside the city. And there may be a speed limit on the road itself, where it will be entered into the corresponding tag - maxspeed. Well, at the same time it is possible to get a sign of one-sidedness from the oneway tag - in vain, is it possible, in the format, the field is reserved for it?

We merge with the table of polygons - it contains administrative boundaries, and a table of lines - lines of roads are stored in it.

select ST_X(transform(p.way,94326)) as X,
ST_Y(transform(p.way,94326)) as Y, '102' as TYPE,
case when l.maxspeed is not NULL then l.maxspeed else case when t.admin_level = '8' then '60' else '90' end end as SPEED,
case when l.oneway = 'yes' then '1' else '2' end as DirType

from (planet_osm_point p
join planet_osm_line l on (l.highway is not NULL and ST_DWithin(p.way,l.way,.1)))
LEFT OUTER JOIN planet_osm_polygon t on (t.admin_level = '8' and ST_Within(p.way, t.way))
where p.traffic_calming is not NULL;


Using the standard SQL join, we compare all three tables — we find each speed bump a couple of lines in no more than ten centimeters from it, and the city polygon surrounding it, if possible. Next, we try to fill in the speed field with the maxspeed line, and if it isn’t filled in, we’ll see if the dot has fallen into the administrative boundary of the city: if yes, set 60, if not - 90.

idx - the sequence number - was tempted to take from osm_id. But then it became clear that a policeman who was lying at the junction of two lines would meet in both lines in both directions, and his id would be the same for different directions, and therefore the best option seemed to be to number them all in order.

Fear and lack of documentation

Since Navitel is a commercial product, it has neither the source code, nor the normal demo version, the last figure left without comment is the direction. No hints about how it is considered could not be found, except that it is clearly a figure in degrees.
“Well, since all the maps are in the merkatorovskaya projection, and it preserves the angles and directions , it means that they are counted in it!”

To calculate the azimuth, PostGIS requires two points. Obviously, they should lie on the same road, on which the speed bump is installed, at some distance on both sides of it.
The ST_Line_Locate_Point method allowed finding the length along the line from its beginning to the specified point, and ST_Line_Interpolate_Point - on the contrary, finding the point along the length along the line.

But psql draws in response to requests for me a tablet in pseudographics, when I need CSV! Is it possible to write a wrapper?
It turned out that everything is much simpler. The COPY function allows you to deliver the result of the request immediately to CSV - just what you need. Here is such a little cat joy.

I expel the first result, I find the owner of Navitel, forcing to check ... and oops. Directions are not parallel to the roads. The thought “really ...” arises, and the paws themselves write reprojection in 4326 when calculating the azimuth.

Yes, the glorious developers of Navitel consider a direction in a projection that does not preserve directions. Honor them, praise and honor.

Quasiitog

As a result, the request was such a monster, extracting from the data of OpenStreetMap what it originally seemed to be not in them - speed bumps in the Navitel format. You can write a request to get speed control cameras from them as homework. :)
CREATE TEMP SEQUENCE idx ;
COPY (

select nextval('idx') as idx,
ST_X(transform(p.way,94326)) as X,
ST_Y(transform(p.way,94326)) as Y, '102' as TYPE,
case when l.maxspeed is not NULL then l.maxspeed else case when t.admin_level = '8' then '60' else '90' end end as SPEED,
case when l.oneway = 'yes' then '1' else '2' end as DirType,
floor(ST_Azimuth(
transform(ST_Line_Interpolate_Point(l.way,
case when (ST_Line_Locate_Point(l.way,p.way)-0.01 < 0) then 0 else ST_Line_Locate_Point(l.way,p.way)-0.0000001 end),94326),
transform(ST_Line_Interpolate_Point(l.way,
case when (ST_Line_Locate_Point(l.way,p.way)+0.01 > 1) then 1 else ST_Line_Locate_Point(l.way,p.way)+0.0000001 end)
,94326))/(2*pi())*360) as Direction
from (planet_osm_point p
join planet_osm_line l on (l.highway is not NULL and ST_DWithin(p.way,l.way,.1)))
LEFT OUTER JOIN planet_osm_polygon t on (t.admin_level = '8' and ST_Within(p.way, t.way))
where p.traffic_calming is not NULL
) TO STDOUT WITH CSV HEADER;


But Kotyara would not be Kotyara if his cat ears did not stick out everywhere. Therefore, as a warm-up, in addition to the every-minute unloading of updated speed bumps and cameras for the CIS , the rendering of all this disgrace into tiles was organized so that those who do not have a navigator with Navitel can admire it.

Closer to the weekend, comrade andrewsh (I think he also has something to tell, if you ask him :) he wrote a convenient interface for reporting new speed bumps without registration - latlon.org/tc , and it was decided to put all the services to the people. I hope you will like it.

Have fun! ;)