addons-contrib: objects.link/unlink syntax update
[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     context.collection.objects.link(obj)
242     context.view_layer.objects.active = obj
243     obj.select_set(True)
244
245     # entry points for other addons
246     obj['osmfile'] = filepath
247     obj['tagged'] = tag
248     obj['utm'] = utm
249
250     print("Parse done... %d" % tidx)
251
252     return {'FINISHED'}
253
254
255 # given lat and longitude in degrees, returns x and y in UTM kilometers.
256 # accuracy : supposed to be centimeter :)
257 # http://fr.wikipedia.org/wiki/Projection_UTM
258 # http://fr.wikipedia.org/wiki/WGS_84
259 # http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
260 # http://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/alg0071.pdf
261 def geoToUTM(lon, lat):
262
263     # if abs(lat) > 80 : lat = 80 #wrong coords.
264
265     # UTM zone, longitude origin, then lat lon in radians
266     z = int((lon + 180) / 6) + 1
267     lon0 = radians(6 * z - 183)
268     lat = radians(lat)
269     lon = radians(lon)
270
271     # CONSTANTS (see refs.)
272     # rayon de la terre √† l'√©quateur
273     a = 6378.137
274     K0 = 0.9996
275     # flattening consts
276     f  = 0.0033528106647474805   # 1 / 298.257223563
277     e2 = 0.0066943799901413165   # 2*f - f**2
278     e4 = 4.481472345240445e-05   # e2**2
279     e6 = 3.0000678794349315e-07  # e2**3
280
281     # lat0. 10000 for South, 0 for North
282     N0 = 10000 if lat < 0 else 0
283
284     # wiki is your friend (don't ask me Im just a writing monkey.)
285     A = (lon - lon0) * cos(lat)
286     C = (e2 / (1 - e2)) * cos(lat) ** 2
287     T = tan(lat) ** 2.0
288     vlat = 1.0 / sqrt(1.0 - e2 * sin(lat) ** 2.0)
289     slat = ((1.0 - (e2    / 4.0) - ((3.0 * e4)    / 64)   - ((5.0  * e6) / 256.0))  * lat -
290             (((3.0 * e2)  / 8.0) + ((3.0 * e4)    / 32.0) + ((45.0 * e6) / 1024.0)) * sin(lat * 2.0) +
291             (((15.0 * e4) / 256.0) + ((45.0 * e6) / 1024.0)) *
292             sin(lat * 4.0) - ((35.0 * e6) / 3072.0) * sin(lat * 6.0))
293
294     E = (500.0 + (K0 * a * vlat) * (A + (1.0 - T + C) *
295          ((A ** 3.0) / 6.0) + (5.0 - 18.0 * T + T**2) * ((A ** 5.0) / 120.0)))
296     N = (N0 + (K0 * a) * (slat + vlat * tan(lat) *
297          (A ** 2.0 / 2.0 + (5.0 - T + 9.0 * C + 4.0 * C ** 2.0) *
298          (A ** 4.0 / 24.0) + (61.0 - 58.0 * T + T ** 2) * A ** 6.0 / 720.0)))
299     return E, N
300
301 ## for testing
302 #if __name__ == "__main__":
303 #    read("/data/downloads/osm_parser/map.osm", bpy.context)
304
305
306 # ----------------------------------------------------------------------------
307 # blender integration
308
309 from bpy.types import Operator
310 from bpy_extras.io_utils import ImportHelper
311
312 from bpy.props import StringProperty, FloatProperty, BoolProperty
313
314
315 class ImportOSM(Operator, ImportHelper):
316     """Import OSM"""
317     #bl_idname = "import.open_street_map"
318     bl_idname = "import_mesh.osm"
319     bl_label = "Import OpenStreetMap (.osm)"
320
321     # ExportHelper mixin class uses this
322     filename_ext = ".osm"
323     filter_glob: StringProperty(
324             default="*.osm",
325             options={'HIDDEN'},
326             )
327     # List of operator properties, the attributes will be assigned
328     # to the class instance from the operator settings before calling.
329     scale: FloatProperty(
330             name="Scale",
331             default=100.0,
332             )
333     utm: BoolProperty(
334             name="in UTM coordinates",
335             default=True,
336             )
337     tag: BoolProperty(
338             name="retrieve .osm tags as vertex groups",
339             default=False,
340             )
341
342     def execute(self, context):
343         return read(context, self.filepath, self.scale, self.tag, self.utm)
344
345
346 # Only needed if you want to add into a dynamic menu
347 def menu_func_export(self, context):
348     self.layout.operator(ImportOSM.bl_idname)
349
350
351 def register():
352     bpy.utils.register_class(ImportOSM)
353     bpy.types.TOPBAR_MT_file_import.append(menu_func_export)
354
355
356 def unregister():
357     bpy.utils.unregister_class(ImportOSM)
358     bpy.types.TOPBAR_MT_file_import.remove(menu_func_export)