Skip to content
Snippets Groups Projects
io_scene_open_street_map.py 12.3 KiB
Newer Older
  • Learn to ignore specific revisions
  • Campbell Barton's avatar
    Campbell Barton committed
    # ***** BEGIN GPL LICENSE BLOCK *****
    
    #
    # This program is free software; you can redistribute it and/or
    # modify it under the terms of the GNU General Public License
    # as published by the Free Software Foundation; either version 2
    # of the License, or (at your option) any later version.
    #
    # This program is distributed in the hope that it will be useful,
    # but WITHOUT ANY WARRANTY; without even the implied warranty of
    # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    # GNU General Public License for more details.
    #
    # You should have received a copy of the GNU General Public License
    # along with this program; if not, write to the Free Software Foundation,
    # Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    #
    # ***** END GPL LICENCE BLOCK *****
    
    # <pep8 compliant>
    
    bl_info = {
        "name": "Open Street Map (.osm)",
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        "author": "Michael Anthrax Schlachter, ideasman42, littleneo",
        "version": (0, 2),
    
        "blender": (2, 6, 3),
        "location": "File > Import",
        "description": "Load Open Street Map File",
        "wiki_url": "",
        "tracker_url": "",
        "category": "Import-Export"}
    
    # originally written for blender 2.4x by (manthrax _at_ hotmail.com),
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    # updated for blender 2.6x by ideasman42
    
    # If you use it for something cool, send me an email and let me know!
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    # littleneo (v0.2 - september 2012) :
    
    ## added a lat/lon to UTM conversion option geoToUTM() :
    # the generated mesh now almost match the exported osm image of the same area, just play with scale and rotate a bit.
    
    ## added ability to retrieve openstreet tags as vertex groups :
    # add every way id as a vgroup if the tag exists in osmkeys.
    # add nodes id as vgroup  if the tag value pair exists in osmvals.
    # any elligible is added in the corresponding osmkeys or value ('_V_...'  vgroup
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    ## added bpy access to addon : 
    # bpy.ops.import_mesh.osm(filepath="....", scale=100, utm=False, tag=False)
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    ## added custom properties to generated object so it can be used as source for other osm related scripts :
    # osmfile   path to osm xml file
    # tagged   bool. osm tags included
    # utm       bool. default True
    
    # modified the matrix approximation for utm = False
    # ! corrected lat/lon x/y mismatch
    
    import bpy
    from mathutils import Vector, Matrix
    import math
    from math import radians, sin, cos, tan, sqrt
    
    ## add more osm tags here.
    # http://wiki.openstreetmap.org/wiki/Map_Features
    osmkeys = [
        'highway',
        'barrier',
        'wall',
        'cycleway',
        'bicycle',
        'waterway',
        'railway',
        'aeroway',
        'aerialway',
        'power',
        'man_made',
        'building',
        'leisure',
        'amenity',
        'office',
        'shop',
        'craft',
        'emergency',
        'tourism',
        'historic',
        'landuse',
        'military',
        'natural',
    ]
    
    # just a try for nodes, for retrieving v tag values.
    # keyname must exists in osmkeys
    osmvals = {'highway':['traffic_signals']}
    
    # vertex group name -> vertex group index lookup
    grouplookup = {}
    
    def parseBranch(nodes, bm, nmap, obj, scale=100.0, tag=False, UTM=False):
        vgroups = bm.verts.layers.deform.verify()
    
        tidx = 0
        inNode = 0
        dlat = dlong = clat = clong = minlat = maxlat = minlong = maxlong = 0.0
        for node in nodes:
            if node.localName == "bounds":
                if node.hasAttributes():
                    for i in range(node.attributes.length):
                        at = node.attributes.item(i)
                        if at.name == "minlat":
                            minlat = float(at.nodeValue)
                        elif at.name == "minlon":
                            minlong = float(at.nodeValue)
                        elif at.name == "maxlat":
                            maxlat = float(at.nodeValue)
                        elif at.name == "maxlon":
                            maxlong = float(at.nodeValue)
                    dlat = maxlat - minlat
                    dlong = maxlong - minlong
                    clat = (maxlat + minlat) * 0.5
                    clong = (maxlong + minlong) * 0.5
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                    
                    if UTM : 
                        dlong, dlat = geoToUTM(dlong, dlat)
                        clong, clat = geoToUTM(clong, clat)
    
    
                    print(dlat, dlong, clat, clong)
    
            if node.localName == "way":
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                
                wayid = node.getAttribute('id')
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                if tag :
                    group = obj.vertex_groups.new('way_%s'%wayid)
                    gid = len(obj.vertex_groups) -1
    
                '''
                if node.hasAttributes():
                    for i in range(node.attributes.length):
                        at=node.attributes.item(i)
                        print(at.name)
                '''
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                if tag :
                    metagid = []
                    for ch in node.childNodes:
                        if ch.localName == "tag":
                            key = ch.getAttribute('k')
                            if key in osmkeys :
                                metagid.append(grouplookup[key])
    
    
                for ch in node.childNodes:
                    if ch.localName == "nd":
                        for i in range(ch.attributes.length):
                            at = ch.attributes.item(i)
                            if at.name == "ref":
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                                vid = int(at.nodeValue)
                                refs.append(vid)
                                if tag :
                                    vert = nmap[vid]
                                    weigths = vert[vgroups]
                                    weigths[gid] = 1.0
                                    for mid in metagid : weigths[mid] = 1.0
                    
    
                first = 1
                for r in refs:
                    if first == 0:
                        edge = bm.edges.get((nmap[pr], nmap[r]))
                        if edge is None:
                            edge = bm.edges.new((nmap[pr], nmap[r]))
                        del edge  # don't actually use it
                    else:
                        first = 0
                    pr = r
    
            if node.localName == "node":
                if node.hasAttributes():
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
                    nid = node.getAttribute('id')
                    nlong = node.getAttribute('lon')
                    nlat = node.getAttribute('lat')
    
                    # is this test necessary ? maybe for faulty .osm files
                    if (nid != '') and (nlat != '') and (nlong != '') :
    
                        if UTM : 
                            nlong, nlat = geoToUTM(float(nlong),float(nlat))
                        else :
                            nlat = float(nlat)
                            nlong = float(nlong)
                        
                        x = (nlong - clong) * scale / dlat
                        y = (nlat - clat) * scale / dlat
                        vert = bm.verts.new((x, y, 0.0))
                        nmap[int(nid)] = vert
                        if tag :
                            metagid = []
                            for ch in node.childNodes:
                                if ch.localName == "tag":
                                    key = ch.getAttribute('k')
                                    val = ch.getAttribute('v')
                                    if key in osmvals and val in osmvals[key] :
                                        metagid.append(grouplookup[key])
                                        metagid.append(grouplookup['_V_'+val])
                                        weigths = vert[vgroups]
                                        group = obj.vertex_groups.new('node_%s'%nid)
                                        gid = len(obj.vertex_groups) -1
                                        weigths[gid] = 1.0
                                        for mid in metagid : weigths[mid] = 1.0
    
                    else : print('node is missing some elements : %s %s %s'%(nid,nlat,nlong))
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
            tidx += parseBranch(node.childNodes, bm, nmap, obj,scale,tag,UTM)
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    def read(context, filepath, scale=100.0, tag=False, utm=False) :
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        # create mesh
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        name = bpy.path.display_name_from_filepath(filepath)
        me = bpy.data.meshes.new(name)    
        obj = bpy.data.objects.new(name, me)
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        # osm tags option
        if tag :
            tvid = 0
            for gid, grname in enumerate(osmkeys) :
                obj.vertex_groups.new('_'+grname)
                grouplookup[grname] = gid + tvid
                if grname in osmvals :
                    for val in osmvals[grname] :
                        tvid += 1
                        obj.vertex_groups.new('_V_'+val)
                        grouplookup['_V_'+val] = gid + tvid
    
        # get xml then feed bmesh
        print("Reading xml...")
        xmldoc = minidom.parse(filepath)
        
        print("Starting parse: %r..." % filepath)
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        tidx = parseBranch(xmldoc.childNodes, bm, nmap, obj,scale,tag,utm)
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        # fast approximation of utm for not too big area
    
    Campbell Barton's avatar
    Campbell Barton committed
        if utm is False :
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
            global_matrix = Matrix(((0.65, 0.0, 0.0, 0.0),
                                    (0.0, 1.0, 0.0, 0.0),
                                    (0.0, 0.0, 1.0, 0.0),
                                    (0.0, 0.0, 0.0, 1.0)))
            me.transform(global_matrix)
    
    
        # create the object in the scene
        scene = context.scene
        scene.objects.link(obj)
        scene.objects.active = obj
        obj.select = True
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        # entry points for other addons
        obj['osmfile'] = filepath
        obj['tagged'] = tag
        obj['utm'] = utm
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    # given lat and longitude in degrees, returns x and y in UTM kilometers.
    # accuracy : supposed to be centimeter :)
    # http://fr.wikipedia.org/wiki/Projection_UTM
    # http://fr.wikipedia.org/wiki/WGS_84
    # http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
    # http://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/alg0071.pdf
    def geoToUTM(lon, lat) :
    
        # if abs(lat) > 80 : lat = 80 #wrong coords.
        
        # UTM zone, longitude origin, then lat lon in radians
        z = int( (lon + 180)  / 6 ) + 1
        lon0 = radians(6*z - 183)
        lat = radians(lat)
        lon = radians(lon)
    
        # CONSTANTS (see refs.)
        # rayon de la terre à l'équateur
        a = 6378.137
        K0 = 0.9996
        # flattening consts
        f  = 0.0033528106647474805  # 1 / 298.257223563
        e2 = 0.0066943799901413165  # 2*f - f**2
        e4 = 4.481472345240445e-05  # e2**2
        e6 = 3.0000678794349315e-07 # e2**3
    
        # lat0. 10000 for South, 0 for North
        N0 = 10000 if lat < 0 else 0
    
        # wiki is your friend (don't ask me Im just a writing monkey.)
        A = (lon - lon0) * cos(lat)
        C = (e2 / (1 - e2)) * cos(lat)**2
        T = tan(lat)**2
        vlat = 1 / sqrt( 1 - e2 * sin(lat)**2 )
        slat = (1-(e2/4)-((3*e4)/64)-((5*e6)/256))*lat - (((3*e2)/8)+((3*e4)/32)+((45*e6)/1024))*sin(lat*2) + (((15*e4)/256) + ((45*e6)/1024) )*sin(lat*4) - ((35*e6)/3072)*sin(lat*6)
        E = 500 + (K0 * a * vlat) * (A + (1-T+C)*((A**3)/6) + (5 - 18 * T + T**2) * ((A**5)/120) )
        N = N0 + (K0 * a) * ( slat+vlat*tan(lat)* (A**2/2 + (5-T+9*C+4*C**2) * (A**4/24) + (61-58*T+T**2) * A**6/720) )
        return E,N
    
    
    ## for testing
    #if __name__ == "__main__":
    #    read("/data/downloads/osm_parser/map.osm", bpy.context)
    
    
    # ----------------------------------------------------------------------------
    # blender integration
    
    from bpy.types import Operator
    from bpy_extras.io_utils import ImportHelper
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
    from bpy.props import StringProperty, FloatProperty, BoolProperty
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        #bl_idname = "import.open_street_map"
        bl_idname = "import_mesh.osm"
    
        bl_label = "Import OpenStreetMap (.osm)"
    
        # ExportHelper mixin class uses this
        filename_ext = ".osm"
    
        filter_glob = StringProperty(
                default="*.osm",
                options={'HIDDEN'},
                )
    
        # List of operator properties, the attributes will be assigned
        # to the class instance from the operator settings before calling.
        scale = FloatProperty(
                name="Scale",
                default=100.0,
                )
    
    
    Jerome Mahieux's avatar
    Jerome Mahieux committed
        utm = BoolProperty(
                name="in UTM coordinates",
                default=True,
                )
    
        tag = BoolProperty(
                name="retrieve .osm tags as vertex groups",
                default=False,
                )
    
        def execute(self, context) :
            return read(context, self.filepath, self.scale, self.tag, self.utm)
    
    
    
    # Only needed if you want to add into a dynamic menu
    def menu_func_export(self, context):
        self.layout.operator(ImportOSM.bl_idname)
    
    
    def register():
        bpy.utils.register_class(ImportOSM)
        bpy.types.INFO_MT_file_import.append(menu_func_export)
    
    
    def unregister():
        bpy.utils.unregister_class(ImportOSM)
        bpy.types.INFO_MT_file_import.remove(menu_func_export)