io_scene_open_street_map.py

# ***** 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)",
    "author": "Michael Anthrax Schlachter, ideasman42, littleneo",
    "version": (0, 2),
    "blender": (2, 63, 0),
    "location": "File > Import",
    "description": "Load Open Street Map File",
    "doc_url": "",
    "tracker_url": "https://developer.blender.org/maniphest/task/edit/form/2/",
    "category": "Import-Export"}

# originally written for blender 2.4x by (manthrax _at_ hotmail.com),
# updated for blender 2.6x by ideasman42
# If you use it for something cool, send me an email and let me know!

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
    dlong = clat = clong = minlat = maxlat = minlong = maxlong = 0.0
    dlat = 1.0  # avoid divide by zero

    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

                if UTM:
                    dlong, dlat = geoToUTM(dlong, dlat)
                    clong, clat = geoToUTM(clong, clat)

                print(dlat, dlong, clat, clong)

        if node.localName == "way":
            wayid = node.getAttribute('id')
            nid = None
            refs = []
            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)
            '''

            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":
                            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 = True
            for r in refs:
                if first is False:
                    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 = False
                pr = r

        if node.localName == "node":
            if node.hasAttributes():
                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))

        tidx += 1
        # if tidx > 1000:
        #     break
        tidx += parseBranch(node.childNodes, bm, nmap, obj, scale, tag, UTM)

    return tidx


def read(context, filepath, scale=100.0, tag=False, utm=False):
    import bmesh
    from xml.dom import minidom

    # create mesh
    bm = bmesh.new()
    name = bpy.path.display_name_from_filepath(filepath)
    me = bpy.data.meshes.new(name)
    obj = bpy.data.objects.new(name, me)

    # 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)
    nmap = {}
    tidx = parseBranch(xmldoc.childNodes, bm, nmap, obj, scale, tag, utm)

    bm.to_mesh(me)

    # fast approximation of utm for not too big area
    if utm is False:
        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
    context.collection.objects.link(obj)
    context.view_layer.objects.active = obj
    obj.select_set(True)

    # entry points for other addons
    obj['osmfile'] = filepath
    obj['tagged'] = tag
    obj['utm'] = utm

    print("Parse done... %d" % tidx)

    return {'FINISHED'}


# 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.0
    vlat = 1.0 / sqrt(1.0 - e2 * sin(lat) ** 2.0)
    slat = ((1.0 - (e2    / 4.0) - ((3.0 * e4)    / 64)   - ((5.0  * e6) / 256.0))  * lat -
            (((3.0 * e2)  / 8.0) + ((3.0 * e4)    / 32.0) + ((45.0 * e6) / 1024.0)) * sin(lat * 2.0) +
            (((15.0 * e4) / 256.0) + ((45.0 * e6) / 1024.0)) *
            sin(lat * 4.0) - ((35.0 * e6) / 3072.0) * sin(lat * 6.0))

    E = (500.0 + (K0 * a * vlat) * (A + (1.0 - T + C) *
         ((A ** 3.0) / 6.0) + (5.0 - 18.0 * T + T**2) * ((A ** 5.0) / 120.0)))
    N = (N0 + (K0 * a) * (slat + vlat * tan(lat) *
         (A ** 2.0 / 2.0 + (5.0 - T + 9.0 * C + 4.0 * C ** 2.0) *
         (A ** 4.0 / 24.0) + (61.0 - 58.0 * T + T ** 2) * A ** 6.0 / 720.0)))
    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

from bpy.props import StringProperty, FloatProperty, BoolProperty


class ImportOSM(Operator, ImportHelper):
    """Import OSM"""
    #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,
            )
    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.TOPBAR_MT_file_import.append(menu_func_export)


def unregister():
    bpy.utils.unregister_class(ImportOSM)
    bpy.types.TOPBAR_MT_file_import.remove(menu_func_export)