source: subversion/applications/utils/import/ogr2osm/ogr2osm.py @ 25444

Revision 25444, 25.2 KB checked in by aharvey, 3 years ago (diff)
  • make indentation consistant
  • fix long options (needed a = suffix for ones that accept arguments)
  • recognise file extensions in upper case
  • add MapInfo? TAB/MID/MIF driver as an accepted file type
  • skip features with no geometry (without this check the script may crash)
Line 
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
3
4""" ogr2osm beta
5 
6 (c) Iván Sánchez Ortega, 2009
7 <ivan@sanchezortega.es>
8
9
10 This piece of crap^H^H^H^Hsoftware is supposed to take just about any vector file
11 as an input thanks to the magic of the OGR libraries, and then output a pretty OSM XML
12 file with that data.
13
14 The cool part is that it will detect way segments shared between several ways, so
15 it will build relations outof thin air. This simplifies the structure of boundaries, for
16 example.
17
18 It is also able to translate attributes to tags, though there is only one such translation
19 scheme by now. In order to translate your own datasets, you should have some basic
20 understanding of python programming. See the files in the translation/ directory.
21 
22 An outstanding issue is that elevation in 2.5D features (that can be generated by
23 reprojecting) is ignored completely.
24 
25 Usage: specify a filename to be converted (its extension will be changed to .osm), and the
26 the projection the source data is in. You can specify the source projection by using either
27 an EPSG code or a Proj.4 string.
28 
29 If the projection is not specified, ogr2osm will try to fetch it from the source data. If
30 there is no projection information in the source data, this will assume EPSG:4326 (WGS84
31 latitude-longitude).
32 
33 python ogr2osm.py [options] [filename]
34       
35 Options:       
36   -e, --epsg=...       EPSG code, forcing the source data projection
37   -p, --proj4=...      PROJ4 string, forcing the source data projection
38   -v, --verbose        Shows some seemingly random characters dancing in the screen
39                            for every feature that's being worked on.
40   -h, --help           Show this message
41   -d, --debug-tags     Outputs the tags for every feature parsed
42   -a, --attribute-stats Outputs a summary of the different tags / attributes encountered
43   -t, --translation=... Select the attribute-tags translation method.
44                            See the translations/ diredtory for valid values.
45       
46 (-e and -p are mutually exclusive. If both are specified, only the last one will be
47 taken into account)
48
49 For example, if the shapefile foobar.shp has projection EPSG:23030, do:
50
51 python ogr2osm.py foobar.shp -e 23030
52
53 This will do an in-the-fly reprojection from EPSG:23030 to EPSG:4326, and write a file
54 called "foobar.osm"
55 
56
57#####################################################################################
58#   "THE BEER-WARE LICENSE":                                                        #
59#   <ivan@sanchezortega.es> wrote this file. As long as you retain this notice you  #
60#   can do whatever you want with this stuff. If we meet some day, and you think    #
61#   this stuff is worth it, you can buy me a beer in return.                        #
62#####################################################################################
63"""
64
65
66import sys
67import os
68import getopt
69from SimpleXMLWriter import XMLWriter
70
71try:
72        from osgeo import ogr
73except:
74        import ogr
75
76try:
77        from osgeo import osr
78except:
79        import osr
80
81# Some needed constants
82from ogr import wkbPoint
83from ogr import wkbLineString 
84from ogr import wkbPolygon 
85from ogr import wkbMultiPoint   
86from ogr import wkbMultiLineString 
87from ogr import wkbMultiPolygon   
88from ogr import wkbGeometryCollection   
89
90from ogr import wkbUnknown
91from ogr import wkbNone
92
93from ogr import wkbPoint25D
94from ogr import wkbLineString25D
95from ogr import wkbPolygon25D
96from ogr import wkbMultiPoint25D
97from ogr import wkbMultiLineString25D
98from ogr import wkbMultiPolygon25D
99from ogr import wkbGeometryCollection25D
100
101
102
103# Default options
104sourceEPSG       = 4326
105sourceProj4      = None
106detectProjection = True
107useEPSG          = False
108useProj4         = False
109showProgress     = False
110debugTags        = False
111attributeStats   = False
112translationMethod = None
113
114# Fetch command line parameters: file and source projection
115try:
116        (opts, args) = getopt.getopt(sys.argv[1:], "e:p:hvdat:", ["epsg=","proj4=","help","verbose","debug-tags","attribute-stats","translation="])
117except getopt.GetoptError:
118        print __doc__
119        sys.exit(2)
120for opt, arg in opts:
121        if opt in ("-h", "--help"):
122                print __doc__
123                sys.exit()
124        elif opt in ("-p", "--proj4"):
125                sourceProj4 = arg
126                useProj4 = True
127                useEPSG = False
128                detectProjection = False
129        elif opt in ("-e", "--epsg"):
130                try:
131                        sourceEPSG = int(arg)
132                except:
133                        print "Error: EPSG code must be numeric (e.g. '4326' instead of 'epsg:4326')"
134                        sys.exit(1)
135                detectProjection = False
136                useEPSG = True
137                useProj4 = False
138        elif opt in ("-v", "--verbose"):
139                showProgress=True
140        elif opt in ("-d", "--debug-tags"):
141                debugTags=True
142        elif opt in ("-a", "--attribute-stats"):
143                attributeStats=True
144                attributeStatsTable = {}
145        elif opt in ("-t", "--translation"):
146                translationMethod = arg
147        else:
148                print "Unknown option " + opt
149
150print (opts,args)
151file = args[0]
152fileExtension = file.split('.')[-1].lower()
153
154
155# FIXME: really complete this table
156if fileExtension == 'shp':
157        driver = ogr.GetDriverByName('ESRI Shapefile');
158elif fileExtension == 'tab' or fileExtension == 'mid' or fileExtension == 'mif':
159        driver = ogr.GetDriverByName('MapInfo File');
160elif fileExtension == 'gpx':
161        driver = ogr.GetDriverByName('GPX');
162elif fileExtension == 'dgn':
163        driver = ogr.GetDriverByName('DGN');
164elif fileExtension == 'gml':
165        driver = ogr.GetDriverByName('GML');
166elif fileExtension == 'csv':
167        driver = ogr.GetDriverByName('CSV');
168elif fileExtension == 'sqlite':
169        driver = ogr.GetDriverByName('SQLite');
170elif fileExtension == 'kml':
171        driver = ogr.GetDriverByName('KML');
172#elif fileExtension == 'kmz':
173        #driver = ogr.GetDriverByName('KML');
174else:
175        print "Error: extension " + fileExtension + " is invalid or not implemented yet."
176
177
178# Strip directories from output file name
179slashPosition = file.rfind('/')
180if slashPosition != -1:
181        #print slashPosition
182        outputFile = file[slashPosition+1:]
183        #print outputFile
184        #print len(fileExtension)
185else:
186        outputFile = file
187       
188outputFile = outputFile[: -len(fileExtension) ] + 'osm'
189
190
191# 0 means read-only
192dataSource = driver.Open(file,0); 
193
194if dataSource is None:
195        print 'Could not open ' + file
196        sys.exit(1)
197
198
199print
200print "Preparing to convert file " + file + " (extension is " + fileExtension + ") into " + outputFile
201
202if detectProjection:
203        print "Will try to detect projection from source metadata, or fall back to EPSG:4326"
204elif useEPSG:
205        print "Will assume that source data is in EPSG:" + str(sourceEPSG)
206elif useProj4:
207        print "Will assume that source data has the Proj.4 string: " + sourceProj4
208
209if showProgress:
210        print "Verbose mode is on. Get ready to see lots of dots."
211
212if debugTags:
213        print "Tag debugging is on. Get ready to see lots of stuff."
214
215
216# Some variables to hold stuff...
217nodeIDsByXY  = {}
218nodeTags     = {}
219nodeCoords   = {}
220nodeRefs     = {}
221segmentNodes = {}
222segmentIDByNodes = {}
223segmentRefs  = {}
224areaRings    = {}
225areaTags     = {}
226lineSegments = {}
227lineTags     = {}
228
229# nodeIDsByXY holds a node ID, given a set of coordinates (useful for looking for duplicated nodes)
230# nodeTags holds the tag pairs given a node ID
231# nodeCoords holds the coordinates of a given node ID (redundant if nodeIDsByXY is properly iterated through)
232# nodeRefs holds up the IDs of any segment referencing (containing) a given node ID, as a dictionary
233# segmentNodes holds up the node IDs for a given segment ID
234# segmentIDByNodes holds up segment IDs for a given pair of node IDs (useful for looking for duplicated segments)
235# segmentRefs holds up the IDs of any ways or areas referencing (containing) a given segment ID, as a dictionary with segment IDs as keys, and a boolean as value (the bool is a flag indicating whether the segment is an existing segment, but reversed - will probably screw things up with oneway=yes stuff)
236# areaRings holds up the rings, as a list of segments, for a given area ID
237# areaTags holds up the tags for a given area ID
238# lineSegments and lineTags work pretty much as areaRings and areaTags (only that lineSegments is a list, and areaRings is a list of lists)
239
240
241# Stuff needed for locating translation methods
242if translationMethod:
243        try:
244                sys.path.append(os.getcwd() + "/translations")
245                module = __import__(translationMethod)
246                translateAttributes = module.translateAttributes
247                translateAttributes([])
248        except:
249                print "Could not load translation method " + translationMethod + ". Check the translations/ directory for valid values."
250                sys.exit(-1)
251        print "Successfully loaded " + translationMethod + " translation method."
252else:
253        # If no function has been defined, perform no translation: just copy everything.
254        translateAttributes = lambda(attrs): attrs
255       
256
257
258elementIdCounter = -1
259nodeCount = 0
260segmentCount = 0
261lineCount = 0
262areaCount = 0
263segmentJoinCount = 0
264
265
266print
267print "Parsing features"
268
269
270# Some aux stuff for parsing the features into the data arrays
271
272def addNode(x,y,tags = {}):
273        "Given x,y, returns the ID of an existing node there, or creates it and returns the new ID. Node will be updated with the optional tags."
274        global elementIdCounter, nodeCount, nodeCoords, nodeIDsByXY, nodeTags, nodeCoords
275       
276        if (x,y) in nodeIDsByXY:
277                # Node already exists, merge tags
278                #print
279                #print "Warning, node already exists"
280                nodeID = nodeIDsByXY[(x,y)]
281                try:
282                        nodeTags[nodeID].update(tags)
283                except:
284                        nodeTags[nodeID]=tags
285                return nodeID
286        else:
287                # Allocate a new node
288                nodeID = elementIdCounter
289                elementIdCounter = elementIdCounter - 1
290               
291                nodeTags[nodeID]=tags
292                nodeIDsByXY[(x,y)] = nodeID
293                nodeCoords[nodeID] = (x,y)
294                nodeCount = nodeCount +1       
295                return nodeID
296       
297       
298def lineStringToSegments(geometry,references):
299        "Given a LineString geometry, will create the appropiate segments. It will add the optional tags and will not check for duplicate segments. Needs a line or area ID for updating the segment references. Returns a list of segment IDs."
300        global elementIdCounter, segmentCount, segmentNodes, segmentTags, showProgress, nodeRefs, segmentRefs, segmentIDByNodes
301       
302        result = []
303       
304        (lastx,lasty,z) = geometry.GetPoint(0)
305        lastNodeID = addNode(lastx,lasty)
306       
307        for k in range(1,geometry.GetPointCount()):
308               
309                (newx,newy,z) = geometry.GetPoint(k)
310                newNodeID = addNode(newx,newy)
311               
312                if (lastNodeID, newNodeID) in segmentIDByNodes:
313                        if showProgress: sys.stdout.write(u"-")
314                        segmentID = segmentIDByNodes[(lastNodeID, newNodeID)]
315                        reversed = False
316                        #print
317                        #print "Duplicated segment"
318                elif (newNodeID, lastNodeID) in segmentIDByNodes:
319                        if showProgress: sys.stdout.write(u"_")
320                        segmentID = segmentIDByNodes[(newNodeID, lastNodeID)]
321                        reversed = True
322                        #print
323                        #print "Duplicated reverse segment"
324                else:
325                        if showProgress: sys.stdout.write('.')
326                        segmentID = elementIdCounter
327                       
328                        elementIdCounter = elementIdCounter - 1
329                        segmentCount = segmentCount +1
330                        segmentNodes[segmentID] = [ lastNodeID, newNodeID ]
331                        segmentIDByNodes[(lastNodeID, newNodeID)] = segmentID
332                        reversed = False
333                       
334                        try:
335                                nodeRefs[lastNodeID].update({segmentID:True})
336                        except:
337                                nodeRefs[lastNodeID]={segmentID:True}
338                        try:
339                                nodeRefs[newNodeID].update({segmentID:True})
340                        except:
341                                nodeRefs[newNodeID]={segmentID:True}
342               
343               
344                try:
345                        segmentRefs[segmentID].update({references:reversed})
346                except:
347                        segmentRefs[segmentID]={references:reversed}
348
349                result.append(segmentID)
350               
351                # FIXME
352                segmentRefs
353               
354                lastNodeID = newNodeID
355        return result
356
357
358
359
360
361
362
363# Let's dive into the OGR data source and fetch the features
364
365for i in range(dataSource.GetLayerCount()):
366        layer = dataSource.GetLayer(i)
367        layer.ResetReading()
368       
369        spatialRef = None
370        if detectProjection:
371                spatialRef = layer.GetSpatialRef()
372                if spatialRef != None:
373                        print "Detected projection metadata:"
374                        print spatialRef
375                else:
376                        print "No projection metadata, falling back to EPSG:4326"
377        elif useEPSG:
378                spatialRef = osr.SpatialReference()
379                spatialRef.ImportFromEPSG(sourceEPSG)
380        elif useProj4:
381                spatialRef = osr.SpatialReference()
382                spatialRef.ImportFromProj4(sourceProj4)
383               
384       
385       
386        if spatialRef == None:  # No source proj specified yet? Then default to do no reprojection.
387                # Some python magic: skip reprojection altogether by using a dummy lamdba funcion. Otherwise, the lambda will be a call to the OGR reprojection stuff.
388                reproject = lambda(geometry): None
389        else:
390                destSpatialRef = osr.SpatialReference()
391                destSpatialRef.ImportFromEPSG(4326)     # Destionation projection will *always* be EPSG:4326, WGS84 lat-lon
392                coordTrans = osr.CoordinateTransformation(spatialRef,destSpatialRef)
393                reproject = lambda(geometry): geometry.Transform(coordTrans)
394               
395
396        featureDefinition = layer.GetLayerDefn()
397       
398        fieldNames = []
399        fieldCount = featureDefinition.GetFieldCount();
400       
401        for j in range(fieldCount):
402                #print featureDefinition.GetFieldDefn(j).GetNameRef()
403                fieldNames.append (featureDefinition.GetFieldDefn(j).GetNameRef())
404                if attributeStats:
405                        attributeStatsTable.update({featureDefinition.GetFieldDefn(j).GetNameRef():{} })
406       
407        print
408        print fieldNames
409        print "Got layer field definitions"
410       
411        #print "Feature definition: " + str(featureDefinition);
412       
413        for j in range(layer.GetFeatureCount()):
414                feature = layer.GetNextFeature()
415                geometry = feature.GetGeometryRef()
416               
417                if geometry == None:
418                        continue
419               
420                fields = {}
421               
422                for k in range(fieldCount-1):
423                        #fields[ fieldNames[k] ] = feature.GetRawFieldRef(k)
424                        fields[ fieldNames[k] ] = feature.GetFieldAsString(k)
425                        if attributeStats:
426                                try:
427                                        attributeStatsTable[ fieldNames[k] ][ feature.GetFieldAsString(k) ] = attributeStatsTable[ fieldNames[k] ][ feature.GetFieldAsString(k) ] + 1
428                                except:
429                                        attributeStatsTable[ fieldNames[k] ].update({ feature.GetFieldAsString(k) : 1})
430
431               
432                # Translate attributes into tags, as defined per the selected translation method
433                tags = translateAttributes(fields)
434               
435                if debugTags:
436                        print
437                        print tags
438               
439                # Do the reprojection (or pass if no reprojection is neccesary, see the lambda function definition)
440                reproject(geometry)
441               
442                # Now we got the fields for this feature. Now, let's convert the geometry.
443                # Points will get converted into nodes.
444                # LineStrings will get converted into a set of ways, each having only two nodes
445                # Polygons will be converted into relations
446                # Later, we'll fix the topology and simplify the ways. If a relation can be simplified into a way (i.e. only has one member), it will be. Adjacent segments will be merged if they share tags and direction.
447               
448                # We'll split a geometry into subGeometries or "elementary" geometries: points, linestrings, and polygons. This will take care of OGRMultiLineStrings, OGRGeometryCollections and the like
449               
450                geometryType = geometry.GetGeometryType()
451               
452                subGeometries = []
453               
454                if geometryType == wkbPoint or geometryType == wkbLineString or geometryType == wkbPolygon:
455                        subGeometries = [geometry]
456                elif geometryType == wkbMultiPoint or geometryType == wkbMultiLineString or geometryType == wkbMultiPolygon or geometryType == wkbGeometryCollection:
457                        if showProgress: sys.stdout.write('M')
458                        for k in range(geometry.GetGeometryCount()):
459                                subGeometries.append(geometry.GetGeometryRef(k))
460                               
461                elif geometryType == wkbPoint25D or geometryType == wkbLineString25D or geometryType == wkbPolygon25D:
462                        if showProgress: sys.stdout.write('z')
463                        subGeometries = [geometry]
464                elif geometryType == wkbMultiPoint25D or geometryType == wkbMultiLineString25D or geometryType == wkbMultiPolygon25D or geometryType == wkbGeometryCollection25D:
465                        if showProgress: sys.stdout.write('Mz')
466                        for k in range(geometry.GetGeometryCount()):
467                                subGeometries.append(geometry.GetGeometryRef(k))
468                               
469                elif geometryType == wkbUnknown:
470                        print "Geometry type is wkbUnknown, feature will be ignored\n"
471                elif geometryType == wkbNone:
472                        print "Geometry type is wkbNone, feature will be ignored\n"
473                else:
474                        print "Unknown or unimplemented geometry type :" + str(geometryType) + ", feature will be ignored\n"
475               
476               
477                for geometry in subGeometries:
478                        if geometry.GetDimension() == 0:
479                                # 0-D = point
480                                if showProgress: sys.stdout.write(',')
481                                x = geometry.GetX()
482                                y = geometry.GetY()
483                               
484                                nodeID = addNode(x,y,tags)
485                                # TODO: tags
486                               
487                        elif geometry.GetDimension() == 1:
488                                # 1-D = linestring
489                                if showProgress: sys.stdout.write('|')
490                               
491                                lineID = elementIdCounter
492                                elementIdCounter = elementIdCounter - 1
493                                lineSegments[lineID] = lineStringToSegments(geometry,lineID)
494                                lineTags[lineID] = tags
495                                lineCount = lineCount + 1
496                               
497                        elif geometry.GetDimension() == 2:
498                                # FIXME
499                                # 2-D = area
500                               
501                                if showProgress: sys.stdout.write('O')
502                                areaID = elementIdCounter
503                                elementIdCounter = elementIdCounter - 1
504                                rings = []
505                               
506                                for k in range(0,geometry.GetGeometryCount()):
507                                        if showProgress: sys.stdout.write('r')
508                                        rings.append(lineStringToSegments(geometry.GetGeometryRef(k), areaID))
509                               
510                                areaRings[areaID] = rings
511                                areaTags[areaID]  = tags
512                                areaCount = areaCount + 1
513                                # TODO: tags
514                                # The ring 0 will be the outer hull, any other rings will be inner hulls.
515       
516       
517       
518print
519print "Nodes: " + str(nodeCount)
520print "Way segments: " + str(segmentCount)
521print "Lines: " + str(lineCount)
522print "Areas: " + str(areaCount)
523       
524print
525print "Joining segments"
526
527       
528# OK, all features should be parsed in the arrays by now
529# Let's start to do some topological magic
530
531# We'll iterate through all the lines and areas, then iterate through all the nodes contained there
532# We'll then fetch all segments referencing that node. If a pair of segments share the same references (i.e. they are part of the same line or area), they will be joined as one and de-referenced from that node. Joining segments mean than the concept of segment changes at this point, becoming linestrings or ways.
533# There are some edge cases in which the algorithm may not prove optimal: if a line (or area ring) crosses itself, then the node will have more than two segments referenced to the line (or area), and does NOT check for the optimal one. As a result, lines that cross themselves may be (incorrectly) split into two and merged via a relation. In other words, the order of the points in a line (or ring) may not be kept if the line crosses itself.
534# The algorithm will not check if the node has been de-referenced: instead, it will check for the first and last node of the segments involved - if the segments have already been joined, the check will fail.
535
536
537
538
539def simplifyNode(nodeID):
540        global nodeRefs, segmentNodes, segmentRefs, showProgress, lineSegments, areaRings, segmentJoinCount
541        #for (nodeID, segments) in nodeRefs.items():
542        segments = nodeRefs[nodeID]
543       
544        segmentsJoined = 0
545        #print
546        #print "Node ID: " + str(nodeID)
547        #print "Node references to: " + str(segments)
548       
549        # We have to try all pairs of segments somehow
550        for segmentID1 in segments.copy():
551                for segmentID2 in segments.copy():      # We'll be changing the references, so make sure we iterate through the original list
552                        if segmentID1 != segmentID2:
553                                #print str(segmentID1) + " vs " + str(segmentID2)
554                                try:
555                                        if segmentNodes[segmentID1][-1] == segmentNodes[segmentID2][0] == nodeID and segmentRefs[segmentID1] == segmentRefs[segmentID2] :
556                                               
557                                                #print "Segment " + str(segmentID1) + ": " + str(segmentNodes[segmentID1])
558                                                #print "Segment " + str(segmentID2) + ": " + str(segmentNodes[segmentID2])
559                                               
560                                                #if showProgress: sys.stdout.write('=')
561                                                segmentNodes[segmentID1].extend( segmentNodes[segmentID2][1:] ) # Voila! Joined!
562                                                for nodeShifted in segmentNodes[segmentID2][1:]:        # Replace node references
563                                                        #print "deleting reference from node " + str(nodeShifted) + " to segment " + str(segmentID2) + "; updating to " + str(segmentID1)
564                                                        del nodeRefs[nodeShifted][segmentID2]
565                                                        nodeRefs[nodeShifted].update({segmentID1:True})
566                                               
567                                                # TODO: Check for potential clashes between the references? As in "way X has these segments in the wrong direction". The trivial case for this looks like a topology error, anyway.
568                                                # Anyway, delete all references to the second segment - we're 100% sure that the line or area references the first one 'cause we've checked before joining the segments
569                                                for segmentRef in segmentRefs[segmentID2]:
570                                                        try:
571                                                                lineSegments[segmentRef].remove(segmentID2)
572                                                        except:
573                                                                for ring in areaRings[segmentRef]:
574                                                                        try:
575                                                                                ring.remove(segmentID2)
576                                                                        except:
577                                                                                pass
578                                                       
579                                                       
580                                                del segmentRefs[segmentID2]
581                                               
582                                                del segmentNodes[segmentID2]
583                                                segmentJoinCount = segmentJoinCount +1
584                                                segmentsJoined = segmentsJoined + 1
585                                except:
586                                        pass    # This is due to the node no longer referencing to a segment because we just de-referenced it in a previous pass of the loop; this will be quite common
587       
588        # FIXME: if segmentsJoined > 1, this should mark the node for further testing - It's very likely to be a self-intersection.
589       
590        if showProgress: sys.stdout.write(str(segmentsJoined))
591       
592print
593print "Simplifying line segments"
594for line in lineSegments.values():
595        #print line
596        for segmentID in line:  # No need to check the last segment, it could not be simplyfied
597                #print segmentID
598                #print segmentNodes[segmentID]
599                for nodeID in segmentNodes[segmentID]:
600                        simplifyNode(nodeID)
601                        #simplifyNode(segmentNodes[segmentID][0])       # last node in segment 
602
603print
604print "Simplifying area segments"
605for area in areaRings.values():
606        for ring in area:
607                for segmentID in ring:
608                        for nodeID in segmentNodes[segmentID]:
609                                simplifyNode(nodeID)    # last node in segment 
610
611
612# That *should* do it... but a second pass through all the nodes will really fix things up. I wonder why some nodes are left out of the previous pass
613print
614print "Simplifying remaining nodes"
615for node in nodeRefs.keys():
616        simplifyNode(node)
617
618
619print
620print "Nodes: " + str(nodeCount)
621print "Original way segments: " + str(segmentCount)
622print "Segment join operations: " + str(segmentJoinCount)
623print "Lines: " + str(lineCount)
624print "Areas: " + str(areaCount)
625
626#print nodeRefs
627#print segmentNodes
628#print lineSegments
629#print areaRings
630#print segmentRefs
631
632
633
634print
635print "Generating OSM XML..."
636print "Generating nodes."
637
638
639#w = XMLWriter(sys.stdout)
640w = XMLWriter(open(outputFile,'w'))
641
642w.start("osm", version='0.6', generator='ogr2osm')
643
644# First, the nodes
645for (nodeID,(x,y)) in nodeCoords.items():
646        w.start("node", visible="true", id=str(nodeID), lat=str(y), lon=str(x))
647        for (tagKey,tagValue) in nodeTags[nodeID].items():
648                if tagValue:
649                        w.element("tag", k=tagKey, v=tagValue)
650        w.end("node")
651        if showProgress: sys.stdout.write('.')
652
653
654#print "Generated nodes. On to shared segments."
655
656# Now, the segments used by more than one line/area, as untagged ways
657
658
659#for (segmentID, segmentRef) in segmentRefs.items():
660        #if len(segmentRef) > 1:
661                #print "FIXME: output shared segment"
662                #outputtedSegments[segmentID] = True
663
664
665print
666print "Generated nodes. On to lines."
667
668# Next, the lines, either as ways or as relations
669
670outputtedSegments = {}
671
672for (lineID, lineSegment) in lineSegments.items():
673        if showProgress: sys.stdout.write(str(len(lineSegment)) + " ")
674        if len(lineSegment) == 1:       # The line will be a simple way
675                w.start('way', id=str(lineID), action='modify', visible='true')
676               
677                for nodeID in segmentNodes[ lineSegment[0] ]:
678                        w.element('nd',ref=str(nodeID))
679               
680                for (tagKey,tagValue) in lineTags[lineID].items():
681                        if tagValue:
682                                w.element("tag", k=tagKey, v=tagValue)
683               
684                w.end('way')
685                pass
686        else:   # The line will be a relationship
687                #print
688                #print "Line ID " + str(lineID) + " uses more than one segment: " + str(lineSegment)
689                for segmentID in lineSegment:
690                        if segmentID not in outputtedSegments:
691                                w.start('way', id=str(segmentID), action='modify', visible='true')
692                                for nodeID in segmentNodes[ segmentID ]:
693                                        w.element('nd',ref=str(nodeID))
694                                w.end('way')
695                w.start('relation', id=str(lineID), action='modify', visible='true')
696                for segmentID in lineSegment:
697                        w.element('member', type='way', ref=str(segmentID), role='')
698                for (tagKey,tagValue) in lineTags[lineID].items():
699                        if tagValue:
700                                w.element("tag", k=tagKey, v=tagValue)
701                w.end('relation')
702
703print
704print "Generated lines. On to areas."
705
706# And last, the areas, either as ways or as relations
707
708#print areaRings
709
710for (areaID, areaRing) in areaRings.items():
711        #sys.stdout.write(str(len(areaRings)))
712       
713        if len(areaRing) == 1 and len(areaRing[0]) == 1: # The area will be a simple way
714                w.start('way', id=str(areaID), action='modify', visible='true')
715               
716                for nodeID in segmentNodes[ areaRing[0][0] ]:
717                        w.element('nd',ref=str(nodeID))
718               
719                for (tagKey,tagValue) in areaTags[areaID].items():
720                        if tagValue:
721                                w.element("tag", k=tagKey, v=tagValue)
722               
723                w.end('way')           
724                if showProgress: sys.stdout.write('0 ')
725        else:
726                segmentsUsed = 0
727                segmentsUsedInRing = 0
728                #print "FIXME"
729               
730                for ring in areaRing:
731                        for segmentID in ring:
732                                if segmentID not in outputtedSegments:
733                                        w.start('way', id=str(segmentID), action='modify', visible='true')
734                                        for nodeID in segmentNodes[ segmentID ]:
735                                                w.element('nd',ref=str(nodeID))
736                                        w.end('way')
737                               
738               
739                w.start('relation', id=str(areaID), action='modify', visible='true')
740                w.element("tag", k='type', v='multipolygon')
741               
742                role = 'outer'
743                for ring in areaRing:
744                        for segmentID in ring:
745                                w.element('member', type='way', ref=str(segmentID), role=role)
746                                segmentsUsed = segmentsUsed + 1
747                                segmentsUsedInRing = segmentsUsedInRing + 1
748                        role = 'inner'
749                        #if showProgress: sys.stdout.write(str(segmentsUsedInRing)+'r')
750                        segmentsUsedInRing = 0
751       
752                for (tagKey,tagValue) in areaTags[areaID].items():
753                        if tagValue:
754                                w.element("tag", k=tagKey, v=tagValue)
755                w.end('relation')
756                if showProgress: sys.stdout.write(str(segmentsUsed) + " ")
757
758
759
760
761if attributeStats:
762        print
763        for (attribute, stats) in attributeStatsTable.items():
764                print "All values for attribute " + attribute + ":"
765                print stats
766
767
768print
769print "All done. Enjoy your data!"
770
771
772w.end("osm")
773
774
Note: See TracBrowser for help on using the repository browser.