io_scene_open_street_map.py

   1 # ***** BEGIN GPL LICENSE BLOCK *****
   2 #
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7 #
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software Foundation,
  15 # Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  16 #
  17 # ***** END GPL LICENCE BLOCK *****
  18
  19 bl_info = {
  20     "name": "Open Street Map (.osm)",
  21     "author": "Michael Anthrax Schlachter, ideasman42, littleneo",
  22     "version": (0, 2),
  23     "blender": (2, 63, 0),
  24     "location": "File > Import",
  25     "description": "Load Open Street Map File",
  26     "doc_url": "",
  27     "tracker_url": "https://developer.blender.org/maniphest/task/edit/form/2/",
  28     "category": "Import-Export"}
  29
  30 # originally written for blender 2.4x by (manthrax _at_ hotmail.com),
  31 # updated for blender 2.6x by ideasman42
  32 # If you use it for something cool, send me an email and let me know!
  33
  34 import bpy
  35 from mathutils import Vector, Matrix
  36 import math
  37 from math import radians, sin, cos, tan, sqrt
  38
  39 # add more osm tags here.
  40 # http://wiki.openstreetmap.org/wiki/Map_Features
  41 osmkeys = [
  42     'highway',
  43     'barrier',
  44     'wall',
  45     'cycleway',
  46     'bicycle',
  47     'waterway',
  48     'railway',
  49     'aeroway',
  50     'aerialway',
  51     'power',
  52     'man_made',
  53     'building',
  54     'leisure',
  55     'amenity',
  56     'office',
  57     'shop',
  58     'craft',
  59     'emergency',
  60     'tourism',
  61     'historic',
  62     'landuse',
  63     'military',
  64     'natural',
  65 ]
  66
  67 # just a try for nodes, for retrieving v tag values.
  68 # keyname must exists in osmkeys
  69 osmvals = {'highway': ['traffic_signals']}
  70
  71 # vertex group name -> vertex group index lookup
  72 grouplookup = {}
  73
  74
  75 def parseBranch(nodes, bm, nmap, obj, scale=100.0, tag=False, UTM=False):
  76     vgroups = bm.verts.layers.deform.verify()
  77     tidx = 0
  78     inNode = 0
  79     dlong = clat = clong = minlat = maxlat = minlong = maxlong = 0.0
  80     dlat = 1.0  # avoid divide by zero
  81
  82     for node in nodes:
  83         if node.localName == "bounds":
  84             if node.hasAttributes():
  85                 for i in range(node.attributes.length):
  86                     at = node.attributes.item(i)
  87                     if at.name == "minlat":
  88                         minlat = float(at.nodeValue)
  89                     elif at.name == "minlon":
  90                         minlong = float(at.nodeValue)
  91                     elif at.name == "maxlat":
  92                         maxlat = float(at.nodeValue)
  93                     elif at.name == "maxlon":
  94                         maxlong = float(at.nodeValue)
  95                 dlat = maxlat - minlat
  96                 dlong = maxlong - minlong
  97                 clat = (maxlat + minlat) * 0.5
  98                 clong = (maxlong + minlong) * 0.5
  99
 100                 if UTM:
 101                     dlong, dlat = geoToUTM(dlong, dlat)
 102                     clong, clat = geoToUTM(clong, clat)
 103
 104                 print(dlat, dlong, clat, clong)
 105
 106         if node.localName == "way":
 107             wayid = node.getAttribute('id')
 108             nid = None
 109             refs = []
 110             if tag:
 111                 group = obj.vertex_groups.new('way_%s' % wayid)
 112                 gid = len(obj.vertex_groups) - 1
 113             '''
 114             if node.hasAttributes():
 115                 for i in range(node.attributes.length):
 116                     at=node.attributes.item(i)
 117                     print(at.name)
 118             '''
 119
 120             if tag:
 121                 metagid = []
 122                 for ch in node.childNodes:
 123                     if ch.localName == "tag":
 124                         key = ch.getAttribute('k')
 125                         if key in osmkeys:
 126                             metagid.append(grouplookup[key])
 127
 128             for ch in node.childNodes:
 129                 if ch.localName == "nd":
 130                     for i in range(ch.attributes.length):
 131                         at = ch.attributes.item(i)
 132                         if at.name == "ref":
 133                             vid = int(at.nodeValue)
 134                             refs.append(vid)
 135                             if tag:
 136                                 vert = nmap[vid]
 137                                 weigths = vert[vgroups]
 138                                 weigths[gid] = 1.0
 139                                 for mid in metagid:
 140                                     weigths[mid] = 1.0
 141
 142             first = True
 143             for r in refs:
 144                 if first is False:
 145                     edge = bm.edges.get((nmap[pr], nmap[r]))
 146                     if edge is None:
 147                         edge = bm.edges.new((nmap[pr], nmap[r]))
 148                     del edge  # don't actually use it
 149                 else:
 150                     first = False
 151                 pr = r
 152
 153         if node.localName == "node":
 154             if node.hasAttributes():
 155                 nid = node.getAttribute('id')
 156                 nlong = node.getAttribute('lon')
 157                 nlat = node.getAttribute('lat')
 158
 159                 # is this test necessary ? maybe for faulty .osm files
 160                 if (nid != '') and (nlat != '') and (nlong != ''):
 161
 162                     if UTM:
 163                         nlong, nlat = geoToUTM(float(nlong), float(nlat))
 164                     else:
 165                         nlat = float(nlat)
 166                         nlong = float(nlong)
 167
 168                     x = (nlong - clong) * scale / dlat
 169                     y = (nlat - clat) * scale / dlat
 170                     vert = bm.verts.new((x, y, 0.0))
 171                     nmap[int(nid)] = vert
 172                     if tag:
 173                         metagid = []
 174                         for ch in node.childNodes:
 175                             if ch.localName == "tag":
 176                                 key = ch.getAttribute('k')
 177                                 val = ch.getAttribute('v')
 178                                 if key in osmvals and val in osmvals[key]:
 179                                     metagid.append(grouplookup[key])
 180                                     metagid.append(grouplookup['_V_' + val])
 181                                     weigths = vert[vgroups]
 182                                     group = obj.vertex_groups.new('node_%s' % nid)
 183                                     gid = len(obj.vertex_groups) - 1
 184                                     weigths[gid] = 1.0
 185                                     for mid in metagid:
 186                                         weigths[mid] = 1.0
 187                 else:
 188                     print('node is missing some elements : %s %s %s' % (nid, nlat, nlong))
 189
 190         tidx += 1
 191         # if tidx > 1000:
 192         #     break
 193         tidx += parseBranch(node.childNodes, bm, nmap, obj, scale, tag, UTM)
 194
 195     return tidx
 196
 197
 198 def read(context, filepath, scale=100.0, tag=False, utm=False):
 199     import bmesh
 200     from xml.dom import minidom
 201
 202     # create mesh
 203     bm = bmesh.new()
 204     name = bpy.path.display_name_from_filepath(filepath)
 205     me = bpy.data.meshes.new(name)
 206     obj = bpy.data.objects.new(name, me)
 207
 208     # osm tags option
 209     if tag:
 210         tvid = 0
 211         for gid, grname in enumerate(osmkeys):
 212             obj.vertex_groups.new('_' + grname)
 213             grouplookup[grname] = gid + tvid
 214             if grname in osmvals:
 215                 for val in osmvals[grname]:
 216                     tvid += 1
 217                     obj.vertex_groups.new('_V_' + val)
 218                     grouplookup['_V_' + val] = gid + tvid
 219
 220     # get xml then feed bmesh
 221     print("Reading xml...")
 222     xmldoc = minidom.parse(filepath)
 223
 224     print("Starting parse: %r..." % filepath)
 225     nmap = {}
 226     tidx = parseBranch(xmldoc.childNodes, bm, nmap, obj, scale, tag, utm)
 227
 228     bm.to_mesh(me)
 229
 230     # fast approximation of utm for not too big area
 231     if utm is False:
 232         global_matrix = Matrix(((0.65, 0.0, 0.0, 0.0),
 233                                 (0.0, 1.0, 0.0, 0.0),
 234                                 (0.0, 0.0, 1.0, 0.0),
 235                                 (0.0, 0.0, 0.0, 1.0)))
 236         me.transform(global_matrix)
 237
 238     # create the object in the scene
 239     context.collection.objects.link(obj)
 240     context.view_layer.objects.active = obj
 241     obj.select_set(True)
 242
 243     # entry points for other addons
 244     obj['osmfile'] = filepath
 245     obj['tagged'] = tag
 246     obj['utm'] = utm
 247
 248     print("Parse done... %d" % tidx)
 249
 250     return {'FINISHED'}
 251
 252
 253 # given lat and longitude in degrees, returns x and y in UTM kilometers.
 254 # accuracy : supposed to be centimeter :)
 255 # http://fr.wikipedia.org/wiki/Projection_UTM
 256 # http://fr.wikipedia.org/wiki/WGS_84
 257 # http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
 258 # http://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/alg0071.pdf
 259 def geoToUTM(lon, lat):
 260
 261     # if abs(lat) > 80 : lat = 80 #wrong coords.
 262
 263     # UTM zone, longitude origin, then lat lon in radians
 264     z = int((lon + 180) / 6) + 1
 265     lon0 = radians(6 * z - 183)
 266     lat = radians(lat)
 267     lon = radians(lon)
 268
 269     # CONSTANTS (see refs.)
 270     # rayon de la terre à l'équateur
 271     a = 6378.137
 272     K0 = 0.9996
 273     # flattening consts
 274     f  = 0.0033528106647474805   # 1 / 298.257223563
 275     e2 = 0.0066943799901413165   # 2*f - f**2
 276     e4 = 4.481472345240445e-05   # e2**2
 277     e6 = 3.0000678794349315e-07  # e2**3
 278
 279     # lat0. 10000 for South, 0 for North
 280     N0 = 10000 if lat < 0 else 0
 281
 282     # wiki is your friend (don't ask me Im just a writing monkey.)
 283     A = (lon - lon0) * cos(lat)
 284     C = (e2 / (1 - e2)) * cos(lat) ** 2
 285     T = tan(lat) ** 2.0
 286     vlat = 1.0 / sqrt(1.0 - e2 * sin(lat) ** 2.0)
 287     slat = ((1.0 - (e2    / 4.0) - ((3.0 * e4)    / 64)   - ((5.0  * e6) / 256.0))  * lat -
 288             (((3.0 * e2)  / 8.0) + ((3.0 * e4)    / 32.0) + ((45.0 * e6) / 1024.0)) * sin(lat * 2.0) +
 289             (((15.0 * e4) / 256.0) + ((45.0 * e6) / 1024.0)) *
 290             sin(lat * 4.0) - ((35.0 * e6) / 3072.0) * sin(lat * 6.0))
 291
 292     E = (500.0 + (K0 * a * vlat) * (A + (1.0 - T + C) *
 293          ((A ** 3.0) / 6.0) + (5.0 - 18.0 * T + T**2) * ((A ** 5.0) / 120.0)))
 294     N = (N0 + (K0 * a) * (slat + vlat * tan(lat) *
 295          (A ** 2.0 / 2.0 + (5.0 - T + 9.0 * C + 4.0 * C ** 2.0) *
 296          (A ** 4.0 / 24.0) + (61.0 - 58.0 * T + T ** 2) * A ** 6.0 / 720.0)))
 297     return E, N
 298
 299 ## for testing
 300 #if __name__ == "__main__":
 301 #    read("/data/downloads/osm_parser/map.osm", bpy.context)
 302
 303
 304 # ----------------------------------------------------------------------------
 305 # blender integration
 306
 307 from bpy.types import Operator
 308 from bpy_extras.io_utils import ImportHelper
 309
 310 from bpy.props import StringProperty, FloatProperty, BoolProperty
 311
 312
 313 class ImportOSM(Operator, ImportHelper):
 314     """Import OSM"""
 315     #bl_idname = "import.open_street_map"
 316     bl_idname = "import_mesh.osm"
 317     bl_label = "Import OpenStreetMap (.osm)"
 318
 319     # ExportHelper mixin class uses this
 320     filename_ext = ".osm"
 321     filter_glob: StringProperty(
 322             default="*.osm",
 323             options={'HIDDEN'},
 324             )
 325     # List of operator properties, the attributes will be assigned
 326     # to the class instance from the operator settings before calling.
 327     scale: FloatProperty(
 328             name="Scale",
 329             default=100.0,
 330             )
 331     utm: BoolProperty(
 332             name="in UTM coordinates",
 333             default=True,
 334             )
 335     tag: BoolProperty(
 336             name="retrieve .osm tags as vertex groups",
 337             default=False,
 338             )
 339
 340     def execute(self, context):
 341         return read(context, self.filepath, self.scale, self.tag, self.utm)
 342
 343
 344 # Only needed if you want to add into a dynamic menu
 345 def menu_func_export(self, context):
 346     self.layout.operator(ImportOSM.bl_idname)
 347
 348
 349 def register():
 350     bpy.utils.register_class(ImportOSM)
 351     bpy.types.TOPBAR_MT_file_import.append(menu_func_export)
 352
 353
 354 def unregister():
 355     bpy.utils.unregister_class(ImportOSM)
 356     bpy.types.TOPBAR_MT_file_import.remove(menu_func_export)