add graph reduction docs / scripts

graph_reduction.txt

@@ -0,0 +1,63 @@
+Reducing a rail network to a minimal graph, using PostGIS
+=========================================================
+
+You will need:
+- osmosis (shell commands below are run from the root of the osmosis distribution)
+- PostgreSQL + PostGIS
+- Python
+- psycopg2 (Python PostgreSQL binding)
+
+1) Extract railway-related nodes and ways from your OSM data source
+
+./bin/osmosis --read-pbf wales.osm.pbf --tf accept-ways railway=* --tf reject-relations --used-node --write-xml railway-ways.osm
+./bin/osmosis --read-pbf wales.osm.pbf --tf accept-nodes railway=* --tf reject-ways --tf reject-relations --write-xml railway-nodes.osm
+./bin/osmosis --read-xml railway-ways.osm --read-xml railway-nodes.oxm --merge --write-xml railways_full.osm
+
+2) Create a PostGIS-enabled database
+
+createdb -Upostgres traingraph -T template_postgis
+psql -Upostgres traingraph -f script/pgsimple_schema_0.6.sql 
+
+3) Import data into the database
+
+./bin/osmosis --read-xml ~/Sites/trainhack/railways_wales_fixed.osm --write-pgsimp user="postgres" database="traingraph"
+
+4) Create table of stations
+
+(at the PostgreSQL command line)
+
+SELECT nodes.id, name_tags.v AS name, nodes.geom
+INTO stations
+FROM nodes
+INNER JOIN node_tags AS rail_tags ON (nodes.id=rail_tags.node_id AND rail_tags.k='railway' AND rail_tags.v='station')
+INNER JOIN node_tags AS name_tags ON (nodes.id=name_tags.node_id AND name_tags.k='name');
+
+ALTER TABLE stations ADD PRIMARY KEY (id);
+
+5) Create table of graph edges
+
+SELECT DISTINCT
+	node1.id AS node1_id,
+	node2.id AS node2_id,
+	way_tags.v AS rail_type,
+	ST_MakeLine(node1.geom, node2.geom) AS linestring
+INTO rail_segments
+FROM way_nodes AS way_node1
+INNER JOIN way_nodes AS way_node2 ON (way_node1.way_id = way_node2.way_id AND way_node1.sequence_id + 1 = way_node2.sequence_id)
+INNER JOIN way_tags ON (
+	way_node1.way_id = way_tags.way_id AND way_tags.k = 'railway'
+	AND way_tags.v = 'rail'  -- tweak for other railway types as desired...
+)
+INNER JOIN nodes AS node1 ON (way_node1.node_id = node1.id)
+INNER JOIN nodes AS node2 ON (way_node2.node_id = node2.id);
+
+ALTER TABLE rail_segments ADD FOREIGN KEY (node1_id) REFERENCES nodes(id);
+ALTER TABLE rail_segments ADD FOREIGN KEY (node2_id) REFERENCES nodes(id);
+
+6) Run reduce_graph.py to collapse edges into a minimal graph
+
+python reduce_graph.py
+
+7) Output as KML
+
+python graph_to_kml.py > railways.kml

graph_to_kml.py

@@ -0,0 +1,48 @@
+import psycopg2
+import random
+from xml.sax import saxutils
+
+conn = psycopg2.connect("dbname=traingraph user=postgres")
+cur = conn.cursor()
+
+COLOURS = ['ff0000', '00ff00', '0000ff', 'ffff00', 'ff00ff', '00ffff']
+
+print '''<kml xmlns="http://www.opengis.net/kml/2.2">
+	<Document>
+		<name>traingraph</name>
+		<Style id="dot-icon">
+			<IconStyle>
+				<scale>0.25</scale>
+				<Icon>
+					<href>http://sleeper.demozoo.org/images/dot.png</href>
+				</Icon>
+			</IconStyle>
+		</Style>
+'''
+
+# extract railway lines
+cur.execute('''
+	-- SELECT ST_ASKML(ST_MAKELINE(node1.geom, node2.geom))
+	SELECT ST_ASKML(linestring)
+	FROM rail_segments
+	INNER JOIN nodes AS node1 ON (node1_id = node1.id)
+	INNER JOIN nodes AS node2 ON (node2_id = node2.id)
+''')
+for (kml,) in cur:
+	print '<Placemark>%s<Style><LineStyle><width>2</width><color>ff%s</color></LineStyle></Style></Placemark>' % (kml, random.choice(COLOURS))
+
+# extract stations with wgs84 lng/lat
+cur.execute('''
+	SELECT name, ST_ASKML(geom)
+	FROM stations
+''')
+for (name, kml) in cur:
+	print '<Placemark>'
+	if name:
+		print '<name>%s</name>' % saxutils.escape(name)
+	print '<styleUrl>#dot-icon</styleUrl>'
+	print kml
+	print '</Placemark>'
+
+print '</Document>'
+print '</kml>'

reduce_graph.py

@@ -0,0 +1,97 @@
+import psycopg2
+
+conn = psycopg2.connect("dbname=traingraph user=postgres")
+nodes_qry = conn.cursor()
+cur = conn.cursor()
+
+cur.execute("""
+SELECT nodes.id
+FROM nodes
+INNER JOIN rail_segments ON (nodes.id = node1_id OR nodes.id = node2_id)
+GROUP BY nodes.id
+HAVING count(distinct rail_type) > 1
+""")
+mixed_mode_results = cur.fetchall()
+mixed_mode_nodes = set([res[0] for res in mixed_mode_results])
+
+nodes_qry.execute("""
+	SELECT nodes.id
+	FROM nodes
+	INNER JOIN rail_segments ON (nodes.id = node1_id OR nodes.id = node2_id)
+	GROUP BY nodes.id
+	HAVING COUNT(node1_id) = 2
+""")
+
+for ((node_id,)) in nodes_qry:
+	# find the segments adjoining this one
+	print node_id
+	if node_id in mixed_mode_nodes:
+		continue
+
+	cur.execute("""
+		SELECT node1_id, node2_id, rail_type, linestring
+		FROM rail_segments
+		WHERE (node1_id = %s OR node2_id = %s)
+	""", (node_id, node_id))
+	results = cur.fetchall()
+	if len(results) != 2:
+		# raise Exception('Number of segments adjoining node %s was not exactly 2' % node_id)
+		print "warning: number of segments adjoining node %s was not exactly 2" % node_id
+		continue
+
+	[(seg1_node1, seg1_node2, seg1_type, seg1_linestring), (seg2_node1, seg2_node2, seg2_type, seg2_linestring)] = results
+	if seg1_type != seg2_type:
+		raise Exception('node %s had two neighbouring rail types: %s and %s' % (node_id, seg1_type, seg2_type))
+
+	if seg1_node2 == node_id and seg2_node1 == node_id:
+		# new segment is seg1 + seg2
+		cur.execute("""
+			INSERT INTO rail_segments (node1_id, node2_id, rail_type, linestring)
+			VALUES (
+				%s, %s, %s,
+				ST_LINEMERGE(ST_COLLECT(%s, %s))
+			)
+		""", (seg1_node1, seg2_node2, seg1_type, seg1_linestring, seg2_linestring))
+	elif seg1_node1 == node_id and seg2_node2 == node_id:
+		# new segment is seg2 + seg1
+		cur.execute("""
+			INSERT INTO rail_segments (node1_id, node2_id, rail_type, linestring)
+			VALUES (
+				%s, %s, %s,
+				ST_LINEMERGE(ST_COLLECT(%s, %s))
+			)
+		""", (seg2_node1, seg1_node2, seg1_type, seg2_linestring, seg1_linestring))
+	elif seg1_node2 == node_id and seg2_node2 == node_id:
+		# new segment is seg1 + rev(seg2)
+		cur.execute("""
+			INSERT INTO rail_segments (node1_id, node2_id, rail_type, linestring)
+			VALUES (
+				%s, %s, %s,
+				ST_LINEMERGE(ST_COLLECT(%s, ST_REVERSE(%s)))
+			)
+		""", (seg1_node1, seg2_node1, seg1_type, seg1_linestring, seg2_linestring))
+	elif seg1_node1 == node_id and seg2_node1 == node_id:
+		# new segment is rev(seg2) + seg1
+		cur.execute("""
+			INSERT INTO rail_segments (node1_id, node2_id, rail_type, linestring)
+			VALUES (
+				%s, %s, %s,
+				ST_LINEMERGE(ST_COLLECT(ST_REVERSE(%s), %s))
+			)
+		""", (seg2_node2, seg1_node2, seg1_type, seg2_linestring, seg1_linestring))
+	else:
+		raise Exception("don't know how to reduce node %s!" % node_id)
+
+	cur.execute("""
+		DELETE FROM rail_segments
+		WHERE node1_id = %s AND node2_id = %s and rail_type = %s
+	""", (seg1_node1, seg1_node2, seg1_type))
+	cur.execute("""
+		DELETE FROM rail_segments
+		WHERE node1_id = %s AND node2_id = %s and rail_type = %s
+	""", (seg2_node1, seg2_node2, seg2_type))
+	conn.commit()
+
+nodes_qry.close()
+cur.close()
+conn.close()