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