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