WikiProject Falkland Islands/Converting Government-supplied data
Jump to navigation
Jump to search
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 )