Code cleanup: style
[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": "",
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     dlat = dlong = clat = clong = minlat = maxlat = minlong = maxlong = 0.0
82
83     for node in nodes:
84         if node.localName == "bounds":
85             if node.hasAttributes():
86                 for i in range(node.attributes.length):
87                     at = node.attributes.item(i)
88                     if at.name == "minlat":
89                         minlat = float(at.nodeValue)
90                     elif at.name == "minlon":
91                         minlong = float(at.nodeValue)
92                     elif at.name == "maxlat":
93                         maxlat = float(at.nodeValue)
94                     elif at.name == "maxlon":
95                         maxlong = float(at.nodeValue)
96                 dlat = maxlat - minlat
97                 dlong = maxlong - minlong
98                 clat = (maxlat + minlat) * 0.5
99                 clong = (maxlong + minlong) * 0.5
100
101                 if UTM:
102                     dlong, dlat = geoToUTM(dlong, dlat)
103                     clong, clat = geoToUTM(clong, clat)
104
105                 print(dlat, dlong, clat, clong)
106
107         if node.localName == "way":
108             wayid = node.getAttribute('id')
109             nid = None
110             refs = []
111             if tag:
112                 group = obj.vertex_groups.new('way_%s' % wayid)
113                 gid = len(obj.vertex_groups) - 1
114             '''
115             if node.hasAttributes():
116                 for i in range(node.attributes.length):
117                     at=node.attributes.item(i)
118                     print(at.name)
119             '''
120
121             if tag:
122                 metagid = []
123                 for ch in node.childNodes:
124                     if ch.localName == "tag":
125                         key = ch.getAttribute('k')
126                         if key in osmkeys:
127                             metagid.append(grouplookup[key])
128
129             for ch in node.childNodes:
130                 if ch.localName == "nd":
131                     for i in range(ch.attributes.length):
132                         at = ch.attributes.item(i)
133                         if at.name == "ref":
134                             vid = int(at.nodeValue)
135                             refs.append(vid)
136                             if tag:
137                                 vert = nmap[vid]
138                                 weigths = vert[vgroups]
139                                 weigths[gid] = 1.0
140                                 for mid in metagid:
141                                     weigths[mid] = 1.0
142
143             first = True
144             for r in refs:
145                 if first is False:
146                     edge = bm.edges.get((nmap[pr], nmap[r]))
147                     if edge is None:
148                         edge = bm.edges.new((nmap[pr], nmap[r]))
149                     del edge  # don't actually use it
150                 else:
151                     first = False
152                 pr = r
153
154         if node.localName == "node":
155             if node.hasAttributes():
156                 nid = node.getAttribute('id')
157                 nlong = node.getAttribute('lon')
158                 nlat = node.getAttribute('lat')
159
160                 # is this test necessary ? maybe for faulty .osm files
161                 if (nid != '') and (nlat != '') and (nlong != ''):
162
163                     if UTM:
164                         nlong, nlat = geoToUTM(float(nlong), float(nlat))
165                     else:
166                         nlat = float(nlat)
167                         nlong = float(nlong)
168
169                     x = (nlong - clong) * scale / dlat
170                     y = (nlat - clat) * scale / dlat
171                     vert = bm.verts.new((x, y, 0.0))
172                     nmap[int(nid)] = vert
173                     if tag:
174                         metagid = []
175                         for ch in node.childNodes:
176                             if ch.localName == "tag":
177                                 key = ch.getAttribute('k')
178                                 val = ch.getAttribute('v')
179                                 if key in osmvals and val in osmvals[key]:
180                                     metagid.append(grouplookup[key])
181                                     metagid.append(grouplookup['_V_' + val])
182                                     weigths = vert[vgroups]
183                                     group = obj.vertex_groups.new('node_%s' % nid)
184                                     gid = len(obj.vertex_groups) - 1
185                                     weigths[gid] = 1.0
186                                     for mid in metagid:
187                                         weigths[mid] = 1.0
188                 else:
189                     print('node is missing some elements : %s %s %s' % (nid, nlat, nlong))
190
191         tidx += 1
192         # if tidx > 1000:
193         #     break
194         tidx += parseBranch(node.childNodes, bm, nmap, obj, scale, tag, UTM)
195
196     return tidx
197
198
199 def read(context, filepath, scale=100.0, tag=False, utm=False):
200     import bmesh
201     from xml.dom import minidom
202
203     # create mesh
204     bm = bmesh.new()
205     name = bpy.path.display_name_from_filepath(filepath)
206     me = bpy.data.meshes.new(name)
207     obj = bpy.data.objects.new(name, me)
208
209     # osm tags option
210     if tag:
211         tvid = 0
212         for gid, grname in enumerate(osmkeys):
213             obj.vertex_groups.new('_' + grname)
214             grouplookup[grname] = gid + tvid
215             if grname in osmvals:
216                 for val in osmvals[grname]:
217                     tvid += 1
218                     obj.vertex_groups.new('_V_' + val)
219                     grouplookup['_V_' + val] = gid + tvid
220
221     # get xml then feed bmesh
222     print("Reading xml...")
223     xmldoc = minidom.parse(filepath)
224
225     print("Starting parse: %r..." % filepath)
226     nmap = {}
227     tidx = parseBranch(xmldoc.childNodes, bm, nmap, obj, scale, tag, utm)
228
229     bm.to_mesh(me)
230
231     # fast approximation of utm for not too big area
232     if utm is False:
233         global_matrix = Matrix(((0.65, 0.0, 0.0, 0.0),
234                                 (0.0, 1.0, 0.0, 0.0),
235                                 (0.0, 0.0, 1.0, 0.0),
236                                 (0.0, 0.0, 0.0, 1.0)))
237         me.transform(global_matrix)
238
239     # create the object in the scene
240     scene = context.scene
241     scene.objects.link(obj)
242     scene.objects.active = obj
243     obj.select = 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.INFO_MT_file_import.append(menu_func_export)
354
355
356 def unregister():
357     bpy.utils.unregister_class(ImportOSM)
358     bpy.types.INFO_MT_file_import.remove(menu_func_export)