File indexing completed on 2024-04-14 14:16:42
0001 #!/usr/bin/python3 0002 0003 """ 0004 This script is designed to act as assistance in converting shapefiles 0005 to OpenStreetMap data. This file is optimized and tested with MassGIS 0006 shapefiles, converted to EPSG:4326 before being passed to the script. 0007 You can perform this conversion with 0008 0009 ogr2ogr -t_srs EPSG:4326 new_file.shp old_file.shp 0010 0011 It is expected that you will modify the fixed_tags, tag_mapping, and 0012 boring_tags attributes of this script before running. You should read, 0013 or at least skim, the code up until it says: 0014 0015 DO NOT CHANGE AFTER THIS LINE. 0016 0017 to accommodate your own data. 0018 """ 0019 0020 __author__ = "Christopher Schmidt <crschmidt@crschmidt.net>" 0021 __version__ = "$Id$" 0022 0023 gdal_install = """ 0024 Installing GDAL depends on your platform. Information is available at: 0025 0026 http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries 0027 0028 For Debian-based systems: 0029 0030 apt-get install python-gdal 0031 0032 will usually suffice. 0033 """ 0034 0035 # These tags are attached to all exterior ways. You can put any key/value pairs 0036 # in this dictionary. 0037 0038 fixed_tags = {} 0039 feat_dict = {} 0040 node_dict = {} 0041 non_geom = 0 0042 eflag = False 0043 0044 nodes = [] #(id, lon, lat, tags) 0045 ways = [] #(id, node_refs, tags) 0046 relations = [] #(id, ways) 0047 0048 element_indices = [] 0049 0050 non_polygons = ['Admin-1 aggregation', 'Admin-1 minor island', 'Admin-1 scale rank'] 0051 0052 # Here are a number of functions: These functions define tag mappings. The API 0053 # For these functions is that they are passed the attributes from a feature, 0054 # and they return a list of two-tuples which match to key/value pairs. 0055 0056 def access(data): 0057 """Access restrictions.""" 0058 keys = { 0059 'Y': 'yes', 0060 'N': 'private', 0061 'L': 'restricted' 0062 } 0063 if 'pub_access' in data: 0064 if data['pub_access'] in keys: 0065 return [('access', keys[data['pub_access']])] 0066 return None 0067 0068 def protection(data): 0069 keys = { 0070 'P': 'perpetuity', 0071 'T': 'temporary', 0072 'L': 'limited', 0073 } 0074 if 'lev_prot' in data: 0075 if data['lev_prot'] in keys: 0076 return [('protected', keys[data['lev_prot']])] 0077 return None 0078 0079 def owner_type(data): 0080 """See wiki:Key:ownership""" 0081 keys = { 0082 'F': 'national', 0083 'S': 'state', 0084 'C': 'county', 0085 'M': 'municipal', 0086 'N': 'private_nonprofit', 0087 'P': 'private', 0088 'B': 'public_nonprofit', 0089 'L': 'land_trust', 0090 'G': 'conservation_rganization', 0091 'I': 'inholding', 0092 } 0093 if 'owner_type' in data: 0094 if data['owner_type'] in keys: 0095 return [['ownership', keys[data['owner_type']]]] 0096 0097 def purpose(data): 0098 """Based on a discussion on IRC""" 0099 keys = { 0100 'R': [('leisure', 'recreation_ground')], 0101 'C': [('leisure', 'nature_reserve'), ('landuse', 'conservation')], 0102 'B': [('landuse','conservation'), ('leisure','recreation_ground')], 0103 'H': [('historical', 'yes')], 0104 'A': [('agricultural', 'yes'), ('landuse','farm')], 0105 'W': [('landuse', 'resevoir')], 0106 'S': [('scenic','yes')], 0107 'F': [('landuse','land')], 0108 'Q': [('landuse','conservation')], 0109 'U': [('water','yes')] 0110 } 0111 if 'prim_purp' in data: 0112 if data['prim_purp'] in keys: 0113 return keys[data['prim_purp']] 0114 0115 def bathymetry_map(data): 0116 elevation = 0 0117 if 'depth' in data: 0118 elevation = data['depth'] 0119 tag = [('marble:feature', 'bathymetry'), ('ele', elevation)] 0120 return tag 0121 0122 0123 def road_map(data): 0124 keys = { 0125 'Ferry Route': [('route','ferry')], 0126 'Major Highway': [('highway','motorway')], 0127 'Beltway': [('highway','primary')], 0128 'Track': [('highway','tertiary')], 0129 'Unknown': [('highway','unclassified')], 0130 'Secondary Highway': [('highway','trunk')], 0131 'Bypass': [('highway','secondary')], 0132 'Road': [('highway','primary')] 0133 } 0134 if 'type' in data: 0135 if data['type'] in keys: 0136 return keys[data['type']] 0137 0138 def city_map(data): 0139 population = 0 0140 capital = 'no' 0141 country = '' 0142 if data['featurecla'] == 'Admin-0 capital' or data['featurecla'] == 'Admin-1 capital' or data['featurecla'] == 'Admin-0 region capital' or data['featurecla'] == 'Admin-1 region capital': 0143 capital = 'yes' 0144 if 'pop_max' in data: 0145 population = data['pop_max'] 0146 if 'adm0name' in data: 0147 country = data['adm0name'] 0148 0149 # National capitals in OSM are recongised by tag "admin_level"="2" ( http://wiki.openstreetmap.org/wiki/Key:capital#Using_relations_for_capitals ). 0150 # In Natural Earth .shp files these capitals are tagged as "Admin-0 capital". 0151 if data['featurecla'] == 'Admin-0 capital': 0152 admin_level = "2" 0153 temp = [('is_in:country', country), ('capital', capital), ('admin_level', admin_level), ('population', population), ('place', 'city')] 0154 else: 0155 temp = [('is_in:country', country), ('capital', capital), ('population', population), ('place', 'city')] 0156 0157 return temp 0158 0159 def mountain_map(data): 0160 elevation = 0 0161 if 'elevation' in data: 0162 elevation = data['elevation'] 0163 temp = [('natural', 'peak'), ('ele', elevation)] 0164 return temp 0165 0166 def feature_class(data): 0167 global non_fcla_dict 0168 keys = { 0169 'Lake': [('natural', 'water')], 0170 'Alkaline Lake': [('natural', 'water')], 0171 'Reservoir': [('natural', 'water')], 0172 'Road': [(road_map,None)], 0173 'Ferry': [('route','ferry')], 0174 'River': [('waterway', 'river')], 0175 'Coastline': [('natural', 'coastline')], 0176 'Minor coastline': [('natural', 'coastline')], 0177 'Ocean': [('natural', 'water')], 0178 'Land': [('marble_land', 'landmass')], 0179 'Minor island': [('marble_land', 'landmass'), ('place', 'island')], 0180 'Reefs': [('natural', 'reef')], 0181 'Admin-0 country': [('marble_land', 'landmass')], 0182 'Admin-0 sovereignty': [('marble_land', 'landmass')], 0183 'Admin-0 map unit': [('marble_land', 'landmass')], 0184 'Adm-0 scale ranks': [('marble_land', 'landmass')], 0185 'International boundary (verify)': [('boundary', 'administrative'), ('admin_level', '2')], 0186 'Overlay limit': [('boundary', 'administrative'), ('admin_level', '2')], 0187 'Disputed (please verify)': [('boundary', 'administrative'), ('admin_level', '2')], 0188 'Line of control (please verify)': [('boundary', 'administrative'), ('admin_level', '2')], 0189 'Indefinite (please verify)': [('boundary', 'administrative'), ('admin_level', '2')], 0190 'Lease limit': [('boundary', 'administrative'), ('admin_level', '2')], 0191 'Indeterminant frontier': [('boundary', 'administrative'), ('admin_level', '2')], 0192 'Admin-0 lease': [('marble_land', 'landmass')], 0193 'Admin-0 claim area': [('marble_land', 'landmass')], 0194 'Admin-0 breakaway and disputed': [('marble_land', 'landmass')], 0195 'Admin-0 overlay': [('marble_land', 'landmass')], 0196 'Admin-0 indeterminant': [('marble_land', 'landmass')], 0197 'Admin-1 aggregation': [('boundary', 'administrative'), ('admin_level', '4')], 0198 'Admin-1 minor island': [('boundary', 'administrative'), ('admin_level', '4')], 0199 'Admin-1 scale rank': [('boundary', 'administrative'), ('admin_level', '4')], 0200 'Admin-1 region boundary': [('boundary', 'administrative'), ('admin_level', '4')], 0201 'Railroad': [('railway', 'rail')], 0202 'Railroad ferry': [('route', 'ferry')], 0203 'Urban area': [('settlement', 'yes')], 0204 'Timezone': [('marble_land', 'landmass')], 0205 'Historic place': [(city_map,None)], 0206 'Populated place': [(city_map,None)], 0207 'Scientific station': [(city_map,None)], 0208 'Meteorological Station': [(city_map,None)], 0209 'Admin-0 capital': [(city_map,None)], 0210 'Admin-1 capital': [(city_map,None)], 0211 'Admin-0 region capital': [(city_map,None)], 0212 'Admin-1 region capital': [(city_map,None)], 0213 'Admin-0 capital alt': [(city_map,None)], 0214 'Lake Centerline': [('waterway', 'river')], 0215 'Port': [('harbour', 'yes')], 0216 'Island': [('marble_land', 'landmass'), ('place', 'island')], 0217 'Island group': [('marble_land', 'landmass'), ('place', 'island')], 0218 'Wetlands': [('natural', 'wetland')], 0219 'Basin': [('landuse', 'basin')], 0220 'Desert': [('natural', 'desert')], #no tag support in marble, no clue as to how to render 0221 'Depression': [('natural', 'sinkhole')], #no tag support in marble, no clue as to how to render 0222 'Range/mtn': [('marble_land', 'landmass')], # no clue as to how to render 0223 'Geoarea': [('marble_land', 'landmass')], # no clue as to how to render 0224 'Plain' : [('marble_land', 'landmass')], # no clue as to how to render 0225 'Tundra': [('natural', 'tundra')], #no tag support in marble, no clue as to how to render 0226 'Foothills' : [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0227 'Lowland' : [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0228 'Dragons-be-here': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0229 'Delta': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0230 'Peninsula': [('marble_land', 'landmass')], 0231 'Plateau' : [('marble_land', 'landmass')], 0232 'Pen/cape' : [('natural', 'cape')], #osm tag still in proposal stage, no clue as to how to render 0233 'Gorge': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0234 'Coast': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0235 'Continent': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, redundant data 0236 'Isthmus': [('marble_land', 'landmass')], #no tag support in OSM, nothing in Marble, no clue as to how to render 0237 'Valley': [('natural', 'valley')], #no tag support in marble 0238 'waterfall' : [('waterway', 'waterfall')], #no tag support in marble 0239 'cape': [('natural', 'cape')], #osm tag still in proposal stage, no clue as to how to render 'NONE' 0240 'pole': [], 0241 'plain': [], 0242 'island group': [], 0243 'island': [], 0244 'mountain': [(mountain_map,None)], #this tag has to be mapped to Mountain feature of Marble 0245 'spot elevation': [], 0246 'plateau': [], 0247 'depression': [], 0248 'pass': [], 0249 'sound': [('natural', 'water')], 0250 'river': [('waterway', 'river')], 0251 'generic': [('natural', 'water')], #still have to confirm what generic means exactly in marine_polygons? 0252 'lagoon': [('natural', 'water'), ('water', 'lagoon')], 0253 'reef': [('natural', 'reef')], 0254 'gulf': [('natural', 'water')], 0255 'inlet': [('natural', 'water')], 0256 'strait': [('natural', 'water')], 0257 'bay': [('natural', 'water')], 0258 'fjord': [('natural', 'water')], 0259 'sea': [('natural', 'water')], 0260 'ocean': [('natural', 'water')], 0261 'channel': [('natural', 'water')], 0262 'Playa': [('natural', 'water'), ('water', 'lake'), ('salt', 'yes')], 0263 #glacier:type is a proposed OSM tag - http://wiki.openstreetmap.org/wiki/Proposed_features/Glaciers_tags 0264 'Antarctic Ice Shelf': [('natural', 'glacier'), ('glacier:type','shelf')], 0265 'Antarctic Ice Shelf Edge': [('glacier:edge', 'calving_line')], #marble does not support this osm tag 0266 'Glaciated areas': [('natural', 'glacier')], 0267 'Admin-0 Tiny Countries': [], 0268 'Admin-0 Tiny GeoUnit': [], 0269 'Admin-0 Tiny GeoSubunit': [], 0270 'Admin-0 Tiny Countries Pacific': [], 0271 'Pacific Groupings': [], 0272 'Admin-1 boundary': [('boundary', 'administrative'), ('admin_level', '4')], 0273 'Map unit boundary':[], 0274 'Marine Indicator Treaty':[('boundary', 'administrative'), ('admin_level', '2'), ('maritime', 'yes'), ('border_type', 'territorial')], 0275 'Marine Indicator Median':[('boundary', 'administrative'), ('admin_level', '2'), ('maritime', 'yes'), ('border_type', 'territorial')], 0276 'Boundary Treaty':[('boundary', 'administrative'), ('admin_level', '2'), ('maritime', 'yes'), ('border_type', 'territorial')], 0277 'Marine Indicator 200 mi nl':[('boundary', 'maritime'), ('border_type', 'eez')], 0278 'Marine Indicator Disputed':[('boundary', 'administrative'), ('admin_level', '2'), ('maritime', 'yes'), ('border_type', 'territorial'), ('marble:disputed', 'yes')], 0279 'Claim boundary': [('boundary', 'administrative'), ('admin_level', '2')], 0280 'Reference line': [('boundary', 'administrative'), ('admin_level', '2')], 0281 'Breakaway': [('boundary', 'administrative'), ('admin_level', '2')], 0282 'Elusive frontier': [('boundary', 'administrative'), ('admin_level', '2')], 0283 'Country': [('marble_land', 'landmass')], 0284 '1st Order Admin Lines': [('boundary', 'administrative'), ('admin_level', '4')], 0285 'Claim': [('boundary', 'administrative'), ('admin_level', '4')], 0286 'Airport': [('aeroway', 'aerodrome')], 0287 'Date line': [('marble_line', 'date')], 0288 'Bathymetry': [(bathymetry_map,None)] 0289 } 0290 0291 if 'featurecla' in data: 0292 if data['featurecla'] in feat_dict: 0293 feat_dict[data['featurecla']] += 1 0294 else: 0295 feat_dict[data['featurecla']] = 1 0296 if data['featurecla'] in keys: 0297 if len(keys[data['featurecla']]) == 0: 0298 return keys[data['featurecla']] 0299 if hasattr(keys[data['featurecla']][0][0], '__call__'): 0300 return keys[data['featurecla']][0][0](data) 0301 else: 0302 return keys[data['featurecla']] 0303 else: 0304 if data['featurecla'] in non_fcla_dict: 0305 non_fcla_dict[data['featurecla']] += 1 0306 else: 0307 non_fcla_dict[data['featurecla']] = 1 0308 0309 0310 0311 def name_map(data): 0312 if 'name' in data: 0313 return [('name', data['name'])] 0314 0315 0316 def name_tags(data): 0317 """This function returns two things: a 'pretty' name to use, and 0318 may return a landuse of either 'cemetery' or 'forest' if the name 0319 contains those words; based on evaluation the dataset in question.""" 0320 tags = [] 0321 name = data.get('site_name', None) 0322 if not name: 0323 return 0324 name = name.title() 0325 0326 if "cemetery" in name.lower(): 0327 tags.append(['landuse', 'cemetery']) 0328 elif "forest" in name.lower(): 0329 tags.append(['landuse', 'forest']) 0330 0331 tags.append(['name', name]) 0332 return tags 0333 0334 def cal_date(data): 0335 """Return YYYY-MM-DD or YYYY formatted dates, based on 0336 (m)m/(d)d/yyyy dates""" 0337 date = data.get('cal_date_r', None) 0338 if not date: return 0339 try: 0340 m, d, y = map(int, date.split("/")) 0341 if m == 1 and d == 1: 0342 return [['start_date', '%4i' % y]] 0343 return [['start_date', '%04i-%02i-%02i' % (y, m, d)]] 0344 except: 0345 print("Invalid date: %s" % date) 0346 eflag = True 0347 return None 0348 0349 # The most important part of the code: define a set of key/value pairs 0350 # to iterate over to generate keys. This is a list of two-tuples: first 0351 # is a 'key', which is only used if the second value is a string. In 0352 # that case, it is a map of lowercased fielnames to OSM tag names: so 0353 # fee_owner maps to 'owner' in the OSM output. 0354 0355 # if the latter is callable (has a __call__; is a function), then that 0356 # method is called, passing in a dict of feature attributes with 0357 # lowercased key names. Those functions can then return a list of 0358 # two-tuples to be used as tags, or nothin' to skip the tags. 0359 0360 0361 tag_mapping = [ 0362 ('fee_owner', 'owner'), 0363 ('cal_date', cal_date), 0364 ('pub_access', access), 0365 ('lev_prot', protection), 0366 ('owner_type', owner_type), 0367 ('prim_purp', purpose), 0368 ('site_name', name_tags), 0369 ('featurecla', feature_class), 0370 ('name', name_map) 0371 ] 0372 0373 # These tags are not exported, even with the source data; this should be 0374 # used for tags which are usually calculated in a GIS. AREA and LEN are 0375 # common. 0376 0377 boring_tags = [ 'AREA', 'LEN', 'GIS_ACRES'] 0378 0379 # Namespace is used to prefix existing data attributes. If 'None', or 0380 # '--no-source' is set, then source attributes are not exported, only 0381 # attributes in tag_mapping. 0382 0383 namespace = "natural_earth" 0384 #namespace = None 0385 0386 # Uncomment the "DONT_RUN = False" line to get started. 0387 0388 #DONT_RUN = True 0389 DONT_RUN = False 0390 0391 # =========== DO NOT CHANGE AFTER THIS LINE. =========================== 0392 # Below here is regular code, part of the file. This is not designed to 0393 # be modified by users. 0394 # ====================================================================== 0395 0396 import sys 0397 0398 try: 0399 try: 0400 from osgeo import ogr 0401 except ImportError: 0402 import ogr 0403 eflag = True 0404 except ImportError: 0405 __doc__ += gdal_install 0406 if DONT_RUN: 0407 print(__doc__) 0408 sys.exit(2) 0409 print("OGR Python Bindings not installed.\n%s" % gdal_install) 0410 sys.exit(1) 0411 eflag = True 0412 0413 def close_file(): 0414 """ Internal. Close an open file.""" 0415 global open_file 0416 if not open_file.closed: 0417 open_file.write("</osm>") 0418 open_file.close() 0419 0420 def start_new_file(): 0421 """ Internal. Open a new file, closing existing file if necessary.""" 0422 global open_file, file_counter, node_dict, file_name 0423 file_counter += 1 0424 if open_file: 0425 close_file() 0426 open_file = open("%s.%s.osm" % (file_name, file_counter), "w") 0427 print("<?xml version='1.0' encoding='UTF-8'?>" , end = '\n', file = open_file) 0428 print("<osm version='0.5'>" , end = '\n', file = open_file) 0429 node_dict = {} 0430 0431 def write_osm_files(): 0432 global nodes, ways, relations, element_indices, open_file, file_counter, file_name 0433 current = [0, 0, 0] 0434 previous = [0, 0, 0] 0435 for indices in element_indices: 0436 start_new_file() 0437 current[0] = indices[0] - 1 0438 current[1] = indices[1] - 1 0439 current[2] = indices[2] - 1 0440 while current[0] >= previous[0]: 0441 write_node(nodes[current[0]]) 0442 current[0] -= 1 0443 while current[1] >= previous[1]: 0444 write_way(ways[current[1]]) 0445 current[1] -= 1 0446 while current[2] >= previous[2]: 0447 write_relation_multipolygon(relations[current[2]]) 0448 current[2] -= 1 0449 previous = indices[:] 0450 close_file() 0451 element_indices = [] 0452 0453 def clean_attr(val): 0454 """Internal. Hacky way to make attribute XML safe.""" 0455 val = str(val) 0456 val = val.replace("&", "&").replace("'", """).replace("<", "<").replace(">", ">").strip() 0457 return val 0458 0459 def check_featurecla(f): 0460 """ 0461 Checks if featurecla field is present in the feature f. 0462 If present it implies that shp data is from Natural Earth dataset 0463 """ 0464 if 'featurecla' in f.keys(): 0465 return True 0466 else: 0467 return False 0468 0469 def add_point(f): 0470 """Adds a point geometry to the OSM file""" 0471 global id_counter 0472 airport_metadata = None 0473 pt = f.GetGeometryRef() 0474 if check_featurecla(f): 0475 if f['featurecla'] == 'Airport': 0476 airport_metadata = f 0477 f = None 0478 node_id = add_node(id_counter, pt.GetX(0), pt.GetY(0), 'POINT', f) 0479 if node_id == id_counter: 0480 id_counter += 1 0481 if airport_metadata != None: 0482 add_way_around_node(airport_metadata) 0483 0484 0485 0486 def add_relation_multipolygon(geom, f): 0487 """ Writes the multipolygon relation to the OSM file, returns 0 if no relation is formed""" 0488 global id_counter, file_counter, counter, file_name, open_file, namespace 0489 rel_ways = [] 0490 rel_id = 0 0491 way_id = add_way(geom.GetGeometryRef(0), f, True) 0492 if way_id == None: 0493 print('Error in writing relation') 0494 return None 0495 rel_ways.append(way_id) 0496 0497 if geom.GetGeometryCount() > 1: 0498 for i in range(1, geom.GetGeometryCount()): 0499 way_id = add_way(geom.GetGeometryRef(i), f, False) 0500 if way_id == None: 0501 print('Error in writing relation') 0502 return None 0503 rel_ways.append(way_id) 0504 rel_id = id_counter 0505 if check_featurecla(f): 0506 if f['featurecla'] in non_polygons: 0507 return 0 #means no relation is there 0508 relations.append((rel_id, rel_ways)) 0509 id_counter += 1 0510 0511 return rel_id #if rel_id return 0, means no relations is there 0512 0513 def write_relation_multipolygon(relation): 0514 global open_file 0515 print("<relation id='-%s'><tag k='type' v='multipolygon' />" % relation[0] , end = '\n', file = open_file) 0516 print('<member type="way" ref="-%s" role="outer" />' % relation[1][0] , end = '\n', file = open_file) 0517 for way in relation[1][1:]: 0518 print('<member type="way" ref="-%s" role="inner" />' % way , end = '\n', file = open_file) 0519 print("</relation>" , end = '\n', file = open_file) 0520 0521 def write_tags(f): 0522 """Writes the tags associated with a way or a relation""" 0523 global id_counter, file_counter, counter, file_name, open_file, namespace 0524 field_count = f.GetFieldCount() 0525 fields = {} 0526 for field in range(field_count): 0527 value = f.GetFieldAsString(field) 0528 name = f.GetFieldDefnRef(field).GetName() 0529 if namespace and name and value and name not in boring_tags: 0530 print(" <tag k='%s:%s' v='%s' />" % (namespace, name, clean_attr(value)) , end = '\n', file = open_file) 0531 fields[name.lower()] = value 0532 tags={} 0533 for tag_name, map_value in tag_mapping: 0534 if hasattr(map_value, '__call__'): 0535 tag_values = map_value(fields) 0536 if tag_values: 0537 for tag in tag_values: 0538 tags[tag[0]] = tag[1] 0539 else: 0540 if tag_name in fields: 0541 tags[map_value] = fields[tag_name].title() 0542 for key, value in tags.items(): 0543 if key and value: 0544 print(" <tag k='%s' v='%s' />" % (key, clean_attr(value)) , end = '\n', file = open_file) 0545 for name, value in fixed_tags.items(): 0546 print(" <tag k='%s' v='%s' />" % (name, clean_attr(value)) , end = '\n', file = open_file) 0547 if f.GetGeometryRef().GetGeometryName() == 'POLYGON' or f.GetGeometryRef().GetGeometryName() == 'MULTIPOLYGON': 0548 if check_featurecla(f): 0549 if f['featurecla'] not in non_polygons: 0550 print(" <tag k='area' v='yes' />" , end = '\n', file = open_file) 0551 0552 def add_way(geom, f, tag_flag): 0553 """ Writes the way of a particular geometry to the OSM file""" 0554 global open_file, id_counter, ways 0555 ids = add_way_nodes(geom, f) 0556 if len(ids) == 0: 0557 print('Error in writing way') 0558 return None 0559 way_id = id_counter 0560 id_counter += 1 0561 node_refs = ids 0562 if tag_flag: 0563 tags = f 0564 else: 0565 tags = None 0566 ways.append((way_id, node_refs, tags)) 0567 return way_id 0568 0569 def write_way(way): 0570 global open_file 0571 print("<way id='-%s'>" % way[0] , end = '\n', file = open_file) 0572 for i in way[1]: 0573 print(" <nd ref='-%s' />" % i , end = '\n', file = open_file) 0574 if way[2]: 0575 write_tags(way[2]) 0576 print("</way>", end = '\n', file = open_file) 0577 0578 def add_way_nodes(geom, f): 0579 """Writes the nodes of a particular way""" 0580 global open_file, id_counter 0581 ids = [] 0582 geom_name = geom.GetGeometryName() 0583 pt_count = geom.GetPointCount() 0584 if geom_name == 'LINESTRING': 0585 range_count = range(geom.GetPointCount()) 0586 else: 0587 range_count = range(geom.GetPointCount() - 1) 0588 if range_count == 0 or pt_count == 0: 0589 print( "Degenerate ", geom_name , end = '\n', file = sys.stderr) 0590 return 0591 #if geom_name != 'LINESTRING': 0592 # pt_count -= 1 0593 for count in range(pt_count): 0594 node_id = add_node(id_counter, geom.GetX(count), geom.GetY(count), geom_name, f) 0595 if node_id == id_counter: #means a new node is created, if not means node already exists 0596 id_counter += 1 0597 ids.append(node_id) 0598 return ids 0599 0600 0601 def add_node(num_id, lon, lat, geom_name, f): 0602 """ Writes the node to the OSM file""" 0603 global open_file, node_dict 0604 key = (lon, lat) 0605 if geom_name == 'POINT': 0606 nodes.append((num_id, lon, lat, f)) 0607 node_dict[key] = num_id 0608 else: 0609 if key in node_dict: 0610 num_id = node_dict[key] 0611 else: 0612 nodes.append((num_id, lon, lat, None)) 0613 node_dict[key] = num_id 0614 return num_id 0615 0616 def write_node(node): 0617 global open_file 0618 if node[3] == None: 0619 print("<node id='-{}' visible='true' lon='{:.10f}' lat='{:.10f}' />".format(node[0], node[1], node[2]), end = '\n', file = open_file) 0620 else: 0621 print("<node id='-{}' visible='true' lon='{:.10f}' lat='{:.10f}' >".format(node[0], node[1], node[2]), end = '\n', file = open_file) 0622 write_tags(node[3]) 0623 print("</node>", end = '\n', file = open_file) 0624 0625 def add_way_around_node(f): 0626 """ Writes a way around a single point""" 0627 global id_counter, ways 0628 nid = id_counter - 1 0629 ways.append((id_counter, [nid], f)) 0630 id_counter += 1 0631 0632 0633 open_file = None 0634 0635 file_name = None 0636 0637 id_counter = 1 0638 0639 file_counter = 0 0640 counter = 0 0641 0642 geom_counter = {} 0643 0644 class AppError(Exception): pass 0645 0646 def run(filenames, slice_count=1, obj_count=5000000, output_location=None, no_source=False): 0647 """Run the converter. Requires open_file, file_name, id_counter, 0648 file_counter, counter to be defined in global space; not really a very good 0649 singleton.""" 0650 global id_counter, file_counter, counter, file_name, open_file, namespace, non_geom, non_fcla_dict, nodes, ways, relations, geom_counter 0651 open_file = None 0652 0653 file_name = None 0654 0655 id_counter = 1 0656 0657 file_counter = 0 0658 counter = 0 0659 0660 geom_counter = {} 0661 node_dict = {} 0662 if output_location: 0663 file_name = output_location 0664 # start_new_file() 0665 for filename in filenames: 0666 non_geom = 0 0667 non_fcla_dict = {} 0668 0669 if no_source: 0670 namespace=None 0671 0672 ds = ogr.Open(filename) 0673 if not ds: 0674 raise AppError("OGR Could not open the file %s" % filename) 0675 eflag = True 0676 l = ds.GetLayer(0) 0677 0678 max_objs_per_file = obj_count 0679 0680 extent = l.GetExtent() 0681 #if extent[0] < -180 or extent[0] > 180 or extent[2] < -90 or extent[2] > 90: 0682 # raise AppError("Extent does not look like degrees; are you sure it is? \n(%s, %s, %s, %s)" % (extent[0], extent[2], extent[1], extent[3])) 0683 slice_width = (extent[1] - extent[0]) / slice_count 0684 0685 seen = {} 0686 0687 print("Running %s slices with %s base filename against shapefile %s" % ( 0688 slice_count, file_name, filename)) 0689 0690 for i in range(slice_count): 0691 0692 l.ResetReading() 0693 l.SetSpatialFilterRect(extent[0] + slice_width * i, extent[2], extent[0] + (slice_width * (i + 1)), extent[3]) 0694 0695 #start_new_file() 0696 f = l.GetNextFeature() 0697 0698 obj_counter = 0 0699 last_obj_split = 0 0700 0701 while f: 0702 start_id_counter = id_counter 0703 if f.GetFID() in seen: 0704 f = l.GetNextFeature() 0705 continue 0706 seen[f.GetFID()] = True 0707 0708 if (obj_counter - last_obj_split) > max_objs_per_file: 0709 print("Splitting file with %s objs" % (obj_counter - last_obj_split)) 0710 #start_new_file() 0711 last_obj_split = obj_counter 0712 element_indices.append((len(nodes), len(ways), len(relations))) 0713 node_dict = {} 0714 feat_dict = f.items() 0715 geom = f.GetGeometryRef() 0716 geom_name = geom.GetGeometryName() 0717 if geom_name in geom_counter: 0718 geom_counter[geom_name] += 1 0719 else: 0720 geom_counter[geom_name] = 1 0721 if geom_name == 'POLYGON': 0722 rel_id = add_relation_multipolygon(geom, f) 0723 if rel_id == None: 0724 f = l.GetNextFeature() 0725 continue 0726 elif geom_name == 'LINESTRING': 0727 way_id = add_way(geom, f, True) 0728 if way_id == None: 0729 f = l.GetNextFeature() 0730 continue 0731 elif geom_name == 'MULTILINESTRING': 0732 for i in range(geom.GetGeometryCount()): 0733 way_id = add_way(geom.GetGeometryRef(i), f, True) 0734 if way_id == None: 0735 f = l.GetNextFeature() 0736 continue 0737 elif geom_name == 'MULTIPOLYGON': 0738 for i in range(geom.GetGeometryCount()): 0739 rel_id = add_relation_multipolygon(geom.GetGeometryRef(i), f) 0740 if rel_id == None: 0741 f = l.GetNextFeature() 0742 continue 0743 elif geom_name == 'POINT': 0744 add_point(f) 0745 else: 0746 ids = [] 0747 non_geom += 1 0748 0749 counter += 1 0750 f = l.GetNextFeature() 0751 obj_counter += (id_counter - start_id_counter) 0752 0753 # for node in nodes: 0754 # write_node(node) 0755 # for way in ways: 0756 # write_way(way) 0757 # for relation in relations: 0758 # write_relation_multipolygon(relation) 0759 0760 element_indices.append((len(nodes), len(ways), len(relations))) 0761 write_osm_files() 0762 0763 # close_file() 0764 nodes = [] #(id, lon, lat, tags) 0765 ways = [] #(id, node_refs, tags) 0766 relations = [] #(id, ways) 0767 0768 if __name__ == "__main__": 0769 if DONT_RUN: 0770 print(__doc__) 0771 sys.exit(2) 0772 0773 from optparse import OptionParser 0774 0775 parse = OptionParser(usage="%prog [args] filename.shp", version=__version__) 0776 parse.add_option("-s", "--slice-count", dest="slice_count", 0777 help="Number of horizontal slices of data", default=1, 0778 action="store", type="int") 0779 parse.add_option("-o", "--obj-count", 0780 dest="obj_count", 0781 help="Target Maximum number of objects in a single .osm file", 0782 default=5000000, type="int") 0783 parse.add_option("-n", "--no-source", dest="no_source", 0784 help="Do not store source attributes as tags.", 0785 action="store_true", default=False) 0786 parse.add_option("-l", "--output-location", 0787 dest="output_location", help="base filepath for output files.", 0788 default="poly_output") 0789 (options, args) = parse.parse_args() 0790 0791 if not len(args): 0792 print("No shapefile name given!") 0793 parse.print_help() 0794 sys.exit(3) 0795 0796 kw = {} 0797 for key in ('slice_count', 'obj_count', 'output_location', 'no_source'): 0798 kw[key] = getattr(options, key) 0799 0800 try: 0801 run(args, **kw) 0802 except AppError as E: 0803 print("An error occurred: \n%s" % E) 0804 eflag = True 0805 0806 print() 0807 print('Geometry types present: ') 0808 for key in geom_counter: 0809 print(key, geom_counter[key]) 0810 print() 0811 print('Feature type present: ') 0812 for key in feat_dict: 0813 print(key, feat_dict[key]) 0814 print() 0815 0816 0817 if eflag: 0818 print('Conversion not Successful :') 0819 else: 0820 if len(non_fcla_dict) == 0 and non_geom == 0: 0821 print('Conversion Successful') 0822 else: 0823 print('Conversion not Successful :') 0824 if len(non_fcla_dict) != 0: 0825 print('Unknown features present in SHP file: ', len(non_fcla_dict)) 0826 print() 0827 for key in non_fcla_dict: 0828 print(key, non_fcla_dict[key]) 0829 if non_geom != 0: 0830 print('Unknown geometry present in SHP file: ', non_geom) 0831 0832