BMO
BMO (brest-metropole-oceane) stands for "Brest Métropole Océane". Brest is a city in Britany, France.
In the second part of 2009, BMO gave its data to OSM. A team composed of french OSMers, namely Steven Le Roux, Emilie Laffray, Pieren, Yann Coupin and François Van Der Biest managed to get the data into OSM.
Some History
Describe here how the Urban Community of Brest decided to give its data to OSM
Links to some pictures showing raw data :
- http://www.flickr.com/photos/stevenleroux/3815231303/in/set-72157620241394011/
- http://www.flickr.com/photos/stevenleroux/3946052456/sizes/l/in/set-72157620241394011/
- http://www.flickr.com/photos/stevenleroux/3946052754/sizes/l/in/set-72157620241394011/
Data import
Source SRS = "Lambert Nord France" http://www.spatialreference.org/ref/epsg/27561/
The first step is to convert the data to EPSG:4326 using for instance :
ogr2ogr -s_srs EPSG:27561 -t_srs EPSG:4326 bati_4326.shp bati.shp
NOTE: This is WRONG and leads to erroneous data (~5 meters shift) We should have used this instead :
ogr2ogr -s_srs "+init=IGNF:LAMB1 +wktext" -t_srs EPSG:4326 bati_4326.shp bati.shp
Cadastre Data has been fixed by a bot in april/may 2010 to adress this issue. Vdb 06:50, 18 May 2010 (UTC)
Cadastre data
The shapefiles were converted to OSM file format, using polyshp2osm.py, modified by Emilie Laffray :
python polyshp2osm.py -o 9999999 bati_4326.shp
... generated a poly_output.1.osm file
In polyshp2osm.py, we used some fixed tags, which are :
fixed_tags = {'building': 'yes', 'source': 'Brest Metropole Oceane, 2009'}
Plus, we chose to keep the TYPE attribute, in the bmo namespace. TYPE attribute can have two values : 'D' or 'L', which is the french for "Dur" (heavy) or "Léger" (light).
Then, we used Pieren's Merge software in order to make adjacent polygons use the same nodes:
java -Xmx1300m Merge poly_output.1.osm
... which generated a poly_output.1.osm.new file
The last step is the data upload using the same version of bulk_upload.py as the one used for Corine Land Cover data import (kudos to Yann Coupin):
python bulk_upload.py -i poly_output.1.osm.new -u BMO_2009 -p ***** -c 'Brest Metropole Oceane data import - see https://wiki.openstreetmap.org/wiki/BMO'
The import started on October the 14th, at 22:50, Paris time, and ended the day after, at 23:14
Road network
Differential import
I describe here the process of selecting roads from the original dataset, which are probably not yet on OSM.
The process is based on the idea that roads which are contained in a 10m buffer around existing OSM ways for more than 50% of their length should not be imported. The implementation provided here is pure SQL, using the powerful PostGIS spatial extension.
Some screenshots help get a clear picture of what's going on.
We first create a dedicated "osm" database :
sudo su postgres createdb -E UTF-8 osm createlang plpgsql osm psql -d osm -f /usr/share/postgresql-8.3-postgis/lwpostgis.sql psql -d osm -f /usr/share/postgresql-8.3-postgis/spatial_ref_sys.sql
Alternative depending os which Linux distribution :
psql -d osm -f /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql psql -d osm -f /usr/share/postgresql/8.4/contrib/postgis-1.5/spatial_ref_sys.sql
We then download from the OSM database the most recent ways in the bounding box of interest, using the XAPI:
wget http://www.informationfreeway.org/api/0.6/way[highway=*][bbox=-4.67,48.31,-4.31,48.49] -O /tmp/brest_ways.osm
Next, these data are imported in the newly created db:
osm2pgsql -d osm -H localhost -P 5432 -l /tmp/brest_ways.osm
The -l switch stands for "lon/lat" and should be equivalent to using --proj 4326 (epsg:4326). Unfortunately, at the time of writing, I found out that using the "--proj 4326" switch results in getting erroneous geometries.
If you want to see the data using QGis for instance, you're probably interested in adding a primary key to the planet_osm_line table:
ALTER TABLE planet_osm_line ADD PRIMARY KEY (osm_id);
We now focus on the shapefile we've been given by BMO.
It first needs conversion to EPSG:4326:
ogr2ogr -s_srs "+init=IGNF:LAMB1 +wktext" -t_srs EPSG:4326 voies_4326.shp graphe_voies.shp
Then, it is converted to SQL inserts:
shp2pgsql -I -s 4326 -g geom voies_4326.shp voies_cub > /tmp/voies_cub.sql
Finally, it's imported in the db:
sudo su postgres; psql -d osm -f /tmp/voies_cub.sql
Note that we had to set client_encoding to iso-8859-1 just for this import, because of the dbf file encoding.
Then, we create a new table which contains the buffered geometries around the OSM highways:
CREATE TABLE osm_ways_buffered ( id integer NOT NULL, CONSTRAINT osm_ways_buffered_pkey PRIMARY KEY (id) ) WITH (OIDS=FALSE); SELECT AddGeometryColumn('','osm_ways_buffered','the_geom','4326','GEOMETRY',2); INSERT INTO osm_ways_buffered(id, the_geom) SELECT c.osm_id, buffer(c.way,0.00018) FROM planet_osm_line AS c
A rough calculation shows that 0.00018 degree is about 10m at 45° latitude.
The next step is the creation of an intermediate table to hold the ratio between the length of lines inside the buffered geometries and the total length:
CREATE TABLE cub_ways_ratio ( gid integer NOT NULL, ratio float, CONSTRAINT cub_ways_ratio_pkey PRIMARY KEY (gid) ) WITH (OIDS=FALSE); INSERT INTO cub_ways_ratio (gid, ratio) (SELECT DISTINCT cub.gid, MAX(ST_LENGTH(ST_Intersection(buf.the_geom, cub.geom)) *100 / ST_LENGTH(cub.geom)) FROM osm_ways_buffered AS buf INNER JOIN voies_cub AS cub ON (buf.the_geom && cub.geom AND ST_intersects(buf.the_geom, cub.geom) = true) WHERE cub.geom IS NOT NULL GROUP BY cub.gid);
Finally, if we run this kind of SQL SELECT statement :
SELECT cub.gid, cub.libru as name, cub.geom FROM voies_cub cub LEFT JOIN cub_ways_ratio r ON cub.gid=r.gid WHERE r.ratio < 50
... we miss :
- all the roads which do not intersect at all the buffered geometries.
- and the small roads, whose length is comparable to the buffer width (and for which the ratio is often very high)
The final export to shapefile takes that into account :
pgsql2shp -f /tmp/to_import.shp -g geom osm "SELECT cub.gid, cub.libru as name, cub.geom FROM voies_cub cub LEFT JOIN cub_ways_ratio r ON cub.gid=r.gid WHERE r.ratio < 50 OR (r.ratio between 50 and 85 AND st_length(cub.geom) < 3*0.00018) UNION SELECT cub.gid, cub.libru as name, cub.geom FROM voies_cub cub WHERE NOT EXISTS(SELECT 1 FROM cub_ways_ratio WHERE cub_ways_ratio.gid=cub.gid) "
Global import
BMO - Global Import - Voting principle
Please vote here : BMO - Global Import - voting principle
POI
no yet tackled ...