BMO

From OpenStreetMap Wiki
Jump to navigation Jump to search
Bootstrap file-earmark-arrow-down.svg

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 :

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 ...