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