WikiProject Falkland Islands/Converting Government-supplied data

The Falkland Islands Government supplied data to OpenStreetMap in .dxf format, exported from AutoCAD. This page explains how this data was converted into .osm files for upload to OpenStreetMap.

First, the files were converted into shapefiles using the instructions at Convert dxf File to Shapefile Using Grass

The Shapefiles were then converted to .osm format using the following Python script, adapted from the one at http://boston.freemap.in/osm/files/mgis_to_osm.py

    #!/usr/bin/python
    # Hacked together from script at http://boston.freemap.in/osm/files/mgis_to_osm.py
    # On Fedora 7, needs gdal-python package & create symlink from /usr/lib/libproj.so.0.5.2 to /usr/lib/libproj.so

   VERSION="0.1"

   import ogr

   def parse_shp_for_massgis( filename ):
        #ogr.RegisterAll()

       dr = ogr.GetDriverByName("ESRI Shapefile")
        poDS = dr.Open( filename )

       if poDS == None:
            raise "Open failed."

       poLayer = poDS.GetLayer( 0 )

       poLayer.ResetReading()

       ret = []

       poFeature = poLayer.GetNextFeature()
        while poFeature:
            tags = {}

           # SEGMENT ID
            tags["massgis:seg_id"] = 2
            # STREET ID
            tags["massgis:way_id"] = 1

           # COPY DOWN THE GEOMETRY
            geom = []
                
            rawgeom = poFeature.GetGeometryRef()
            #print dir( rawgeom )
            for i in range( rawgeom.GetPointCount() ):
                geom.append( (rawgeom.GetX(i), rawgeom.GetY(i)) )
        
            ret.append( (geom, tags) )
            poFeature = poLayer.GetNextFeature()
            
        return ret

   import osr

   fk_wkt = \
    """PROJCS["UTM Zone 21, Southern Hemisphere",
        GEOGCS["International 1909 (Hayford)",
            DATUM["D_unknown",
                SPHEROID["intl",6378388,297]],
            PRIMEM["Greenwich",0],
            UNIT["Degree",0.017453292519943295]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",-57],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",10000000],
        UNIT["Meter",1]]"""

   from_proj = osr.SpatialReference()
    from_proj.ImportFromWkt( fk_wkt )

   to_proj = osr.SpatialReference()
    to_proj.SetWellKnownGeogCS( "EPSG:4326" )

   tr = osr.CoordinateTransformation( from_proj, to_proj )

   def unproject( point ):
        pt = tr.TransformPoint( point[0], point[1] )
        return (pt[1], pt[0])

   def round_point( point, accuracy=2 ):
        return tuple( round(x,accuracy) for x in point )

   def compile_nodelist( parsed_massgis, first_id=1 ):
        nodelist = {}
        
        i = first_id
        for geom, tags in parsed_massgis:
            if len( geom )==0:
                continue
            
            for point in geom:
                r_point = round_point( point )
                if r_point not in nodelist:
                    nodelist[ r_point ] = (i, unproject( point ))
                    i += 1
                
        return (i, nodelist)

   def adjacent( left, right ):
        left_left = round_point(left[0])
        left_right = round_point(left[-1])
        right_left = round_point(right[0])
        right_right = round_point(right[-1])
        
        return ( left_left == right_left or
                 left_left == right_right or
                 left_right == right_left or
                 left_right == right_right )
                 
    def glom( left, right ):
        
        left = list( left )
        right = list( right )
        
        left_left = round_point(left[0])
        left_right = round_point(left[-1])
        right_left = round_point(right[0])
        right_right = round_point(right[-1])
        
        if left_left == right_left:
            left.reverse()
            return left[:-1] + right
            
        if left_left == right_right:
            return right[:-1] + left
            
        if left_right == right_left:
            return left[:-1] + right
            
        if left_right == right_right:
            right.reverse()
            return left[:-1] + right
            
        raise 'segments are not adjacent'

   def glom_once( segments ):
        if len(segments)==0:
            return segments
        
        unsorted = list( segments )
        x = unsorted.pop(0)
        
        while unsorted:
            for i in range( len( unsorted ) ):
                y = unsorted[i]
                if adjacent( x, y ):
                    y = unsorted.pop(i)
                    x = glom( x, y )
                    break
            # Sorted and unsorted lists have no adjacent segments
            else:
                break
                
        return x, unsorted
        
    def glom_all( segments ):
        unsorted = segments
        chunks = []
        
        while unsorted:
            chunk, unsorted = glom_once( unsorted )
            chunks.append( chunk )
            
        return chunks
            
                    

   def compile_waylist( parsed_massgis, blank_way_id ):
        waylist = {}
        
        #Group by massgis:way_id
        for geom, tags in parsed_massgis:
            way_key = tags.copy()
            del way_key['massgis:seg_id']
            way_key = ( way_key['massgis:way_id'], tuple( way_key.iteritems() ) )
            
            if way_key not in waylist:
                waylist[way_key] = []
                
            waylist[way_key].append( geom )
        
        ret = {}
        for (way_id, way_key), segments in waylist.iteritems():
            
            if way_id != blank_way_id:
                ret[way_key] = glom_all( segments )
            else:
                ret[way_key] = segments
            
        return ret
                

   import time
    from xml.sax.saxutils import escape
    def massgis_to_osm( massgis_filename, osm_filename, blank_way_id ):
        
        import_guid = time.strftime( '%Y%m%d%H%M%S' )

       print "parsing shpfile"
        parsed_features = parse_shp_for_massgis( massgis_filename )
        
        print "compiling nodelist"
        i, nodelist = compile_nodelist( parsed_features )
        
        print "compiling waylist"
        waylist = compile_waylist( parsed_features, blank_way_id )
        
        print "constructing osm xml file"
        ret = []
        ret.append( "<?xml version='1.0' encoding='UTF-8'?>" )
        ret.append( "<osm version='0.5' generator='JOSM'>" )
        
        for id, (lat, lon) in nodelist.values():
            ret.append( "  <node id='-%d' action='create' visible='true' lat='%f' lon='%f' >" % (id, lat, lon) )
            #ret.append( "    <tag k=\"source\" v=\"massgis_import_v%s_%s\" />" % (VERSION, import_guid) )
            #ret.append( "    <tag k=\"attribution\" v=\"Office of Geographic and Environmental Information (MassGIS)\" />" )
            ret.append( "  </node>" )
            
        for waykey, segments in waylist.iteritems():
            for segment in segments:
                ret.append( "  <way id='-%d' action='modify' visible='true'>" % i )
                
                ids = [ nodelist[ round_point( point ) ][0] for point in segment ]
                for id in ids:
                    ret.append( "    <nd ref='-%d' />" % id )
                    
                # Following two lines add tags to the ways - comment them out if not required, or edit to suit
                ret.append( "    <tag k='highway' v='unsurfaced' />" )
                ret.append( "    <tag k='source' v='Falkland Islands Government (Public Works Department)' />" )
                
                ret.append( "  </way>" )
                
                i += 1
            
        ret.append( "</osm>" )
        
        print "writing to disk"
        fp = open( osm_filename, "w" )
        fp.write( "\n".join( ret ) )
        fp.close()
        
    if __name__ == '__main__':
        import sys, os.path
        if len(sys.argv) < 2:
            print "%s filename.shp [filename.osm]" % sys.argv[0]
            sys.exit()
        shape = sys.argv[1]
        if len(sys.argv) > 2:
            osm = sys.argv[2]
        else:
            osm = shape[:-4] + ".osm" 
        id = os.path.basename(shape).split("_")[-1]    
        massgis_to_osm( shape, osm, id )