From f92baca50ee48053adacdd5b6fe693af67679f8a Mon Sep 17 00:00:00 2001 From: Harald Husum Date: Fri, 16 Sep 2022 19:30:39 +0200 Subject: [PATCH 1/2] Format with black --- n50merge.py | 855 ++++++------ n50osm.py | 3631 ++++++++++++++++++++++++++++----------------------- utm.py | 463 +++---- 3 files changed, 2705 insertions(+), 2244 deletions(-) diff --git a/n50merge.py b/n50merge.py index 428f28d..2bfa3c3 100644 --- a/n50merge.py +++ b/n50merge.py @@ -19,7 +19,7 @@ import_folder = "~/Jottacloud/osm/n50/" # Folder containing import highway files (default folder tried first) -n50_parts = ['coastline', 'water', 'wood', 'landuse'] +n50_parts = ["coastline", "water", "wood", "landuse"] bbox_margin = 1 # meters @@ -29,353 +29,402 @@ # Output message to console -def message (text): - sys.stderr.write(text) - sys.stderr.flush() +def message(text): + + sys.stderr.write(text) + sys.stderr.flush() # Format time -def timeformat (sec): - if sec > 3600: - return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) - elif sec > 60: - return "%i:%02i minutes" % (sec / 60, sec % 60) - else: - return "%i seconds" % sec +def timeformat(sec): + + if sec > 3600: + return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) + elif sec > 60: + return "%i:%02i minutes" % (sec / 60, sec % 60) + else: + return "%i seconds" % sec # Get name or id of municipality from GeoNorge api -def get_municipality_name (query): - - if query.isdigit(): - url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query - else: - url = "https://ws.geonorge.no/kommuneinfo/v1/sok?knavn=" + urllib.parse.quote(query) - - request = urllib.request.Request(url, headers=header) - - try: - file = urllib.request.urlopen(request) - except urllib.error.HTTPError as e: - if e.code == 404: # Not found - sys.exit("\tMunicipality '%s' not found\n\n" % query) - else: - raise - - if query.isdigit(): - result = json.load(file) - file.close() - municipality_name = result['kommunenavnNorsk'] - return (query, municipality_name) - - else: - result = json.load(file) - file.close() - if result['antallTreff'] == 1: - municipality_id = result['kommuner'][0]['kommunenummer'] - municipality_name = result['kommuner'][0]['kommunenavnNorsk'] - return (municipality_id, municipality_name) - else: - municipalities = [] - for municipality in result['kommuner']: - municipalities.append(municipality['kommunenummer'] + " " + municipalities['kommunenavnNorsk']) - sys.exit("\tMore than one municipality found: %s\n\n" % ", ".join(municipalities)) + +def get_municipality_name(query): + + if query.isdigit(): + url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query + else: + url = "https://ws.geonorge.no/kommuneinfo/v1/sok?knavn=" + urllib.parse.quote( + query + ) + + request = urllib.request.Request(url, headers=header) + + try: + file = urllib.request.urlopen(request) + except urllib.error.HTTPError as e: + if e.code == 404: # Not found + sys.exit("\tMunicipality '%s' not found\n\n" % query) + else: + raise + + if query.isdigit(): + result = json.load(file) + file.close() + municipality_name = result["kommunenavnNorsk"] + return (query, municipality_name) + + else: + result = json.load(file) + file.close() + if result["antallTreff"] == 1: + municipality_id = result["kommuner"][0]["kommunenummer"] + municipality_name = result["kommuner"][0]["kommunenavnNorsk"] + return (municipality_id, municipality_name) + else: + municipalities = [] + for municipality in result["kommuner"]: + municipalities.append( + municipality["kommunenummer"] + + " " + + municipalities["kommunenavnNorsk"] + ) + sys.exit( + "\tMore than one municipality found: %s\n\n" % ", ".join(municipalities) + ) # Build dict data structure from XML. # Works for both N50 and OSM. -def prepare_data (root, tree, nodes, ways, relations): - - # Create dict of nodes - # Each node is a tuple of (lon, lat), corresponding to GeoJSON format x,y - - for node in root.iter("node"): - nodes[ node.attrib['id'] ] = { - 'xml': node, - 'coord':( float(node.attrib['lon']), float(node.attrib['lat']) ) - } - - # Create dict of ways - - for way in root.iter("way"): - way_id = way.attrib['id'] - - way_nodes = [] - way_coordinates = [] - incomplete = False - - # Get way nodes + determine if way is complete - - for node in way.iter("nd"): - node_id = node.attrib['ref'] - if node_id in nodes: - way_nodes.append(node_id) - way_coordinates.append(nodes[node_id]['coord']) - else: - incomplete = True - way_nodes = [] - way_coordinates = [] - break - - ways[ way_id ] = { - 'xml': way, - 'incomplete': incomplete, - 'nodes': way_nodes, - 'coordinates': way_coordinates, - 'parents': set() # Built in next loop - } - - # Create dict of relations - - for relation in root.iter("relation"): - relation_id = relation.attrib['id'] - - incomplete = False - members = [] - for member in relation.iter("member"): - way_id = member.attrib['ref'] - members.append(way_id) - if way_id in ways: - ways[ way_id ]['parents'].add(relation_id) # Incomplete members will be missing - else: - incomplete = True - - # Identify relations which are islands - - island = False - for tag in relation.iter("tag"): - if tag.attrib['k'] == "place" and tag.attrib['v'] in ["islet", "island"]: - island = True - break - - if not incomplete: # Do not store incomplete relations - relations[ relation_id ] = { - 'xml': relation, - 'members': members, - 'island': island - } + +def prepare_data(root, tree, nodes, ways, relations): + + # Create dict of nodes + # Each node is a tuple of (lon, lat), corresponding to GeoJSON format x,y + + for node in root.iter("node"): + nodes[node.attrib["id"]] = { + "xml": node, + "coord": (float(node.attrib["lon"]), float(node.attrib["lat"])), + } + + # Create dict of ways + + for way in root.iter("way"): + way_id = way.attrib["id"] + + way_nodes = [] + way_coordinates = [] + incomplete = False + + # Get way nodes + determine if way is complete + + for node in way.iter("nd"): + node_id = node.attrib["ref"] + if node_id in nodes: + way_nodes.append(node_id) + way_coordinates.append(nodes[node_id]["coord"]) + else: + incomplete = True + way_nodes = [] + way_coordinates = [] + break + + ways[way_id] = { + "xml": way, + "incomplete": incomplete, + "nodes": way_nodes, + "coordinates": way_coordinates, + "parents": set(), # Built in next loop + } + + # Create dict of relations + + for relation in root.iter("relation"): + relation_id = relation.attrib["id"] + + incomplete = False + members = [] + for member in relation.iter("member"): + way_id = member.attrib["ref"] + members.append(way_id) + if way_id in ways: + ways[way_id]["parents"].add( + relation_id + ) # Incomplete members will be missing + else: + incomplete = True + + # Identify relations which are islands + + island = False + for tag in relation.iter("tag"): + if tag.attrib["k"] == "place" and tag.attrib["v"] in ["islet", "island"]: + island = True + break + + if not incomplete: # Do not store incomplete relations + relations[relation_id] = { + "xml": relation, + "members": members, + "island": island, + } # Load relevant OSM elements for chosen municipality + def load_osm(): - global osm_root, osm_tree + global osm_root, osm_tree - message ("Load existing OSM elements from Overpass ...\n") + message("Load existing OSM elements from Overpass ...\n") - query = '[timeout:90];(area[ref=%s][admin_level=7][place=municipality];)->.a;' % municipality_id + \ - '( nwr["natural"](area.a); nwr["waterway"](area.a); nwr["landuse"](area.a); nwr["leisure"](area.a);' + \ - 'nwr["aeroway"](area.a); nwr["seamark:type"="rock"](area.a); );' + \ - '(._;>;<;);out meta;' + query = ( + "[timeout:90];(area[ref=%s][admin_level=7][place=municipality];)->.a;" + % municipality_id + + '( nwr["natural"](area.a); nwr["waterway"](area.a); nwr["landuse"](area.a); nwr["leisure"](area.a);' + + 'nwr["aeroway"](area.a); nwr["seamark:type"="rock"](area.a); );' + + "(._;>;<;);out meta;" + ) - request = urllib.request.Request(overpass_api + "?data=" + urllib.parse.quote(query), headers=header) - file = urllib.request.urlopen(request) - data = file.read() - file.close() + request = urllib.request.Request( + overpass_api + "?data=" + urllib.parse.quote(query), headers=header + ) + file = urllib.request.urlopen(request) + data = file.read() + file.close() - osm_root = ET.fromstring(data) - osm_tree = ET.ElementTree(osm_root) + osm_root = ET.fromstring(data) + osm_tree = ET.ElementTree(osm_root) - prepare_data (osm_root, osm_tree, osm_nodes, osm_ways, osm_relations) + prepare_data(osm_root, osm_tree, osm_nodes, osm_ways, osm_relations) - message ("\tLoaded %i ways, %i relations\n" % (len(osm_ways), len(osm_relations))) + message("\tLoaded %i ways, %i relations\n" % (len(osm_ways), len(osm_relations))) - # Get and display top contributors + # Get and display top contributors - users = {} - for element in osm_root: - if "user" in element.attrib: - user = element.attrib['user'] - if element.tag in ["way", "relation"]: - if user not in users: - users[ user ] = 0 - users[ user ] += 1 + users = {} + for element in osm_root: + if "user" in element.attrib: + user = element.attrib["user"] + if element.tag in ["way", "relation"]: + if user not in users: + users[user] = 0 + users[user] += 1 - sorted_users = sorted(users.items(), key=lambda x: x[1], reverse=True) + sorted_users = sorted(users.items(), key=lambda x: x[1], reverse=True) - message ("\tTop contributors:\n") - for i, user in enumerate(sorted_users): - if user[1] > 10 and i < 10 or user[1] >= 100: - message ("\t\t%s (%i)\n" % (user[0], user[1])) + message("\tTop contributors:\n") + for i, user in enumerate(sorted_users): + if user[1] > 10 and i < 10 or user[1] >= 100: + message("\t\t%s (%i)\n" % (user[0], user[1])) # Load N50 import file. # The file should only contain new elements (negative id's). + def load_n50(): - global n50_root, n50_tree + global n50_root, n50_tree - message ("Load N50 import file elements ...\n") + message("Load N50 import file elements ...\n") - if os.path.isfile(filename): - file = open(filename) - else: - file = open(os.path.expanduser(import_folder + filename)) + if os.path.isfile(filename): + file = open(filename) + else: + file = open(os.path.expanduser(import_folder + filename)) - data = file.read() - file.close() + data = file.read() + file.close() - n50_root = ET.fromstring(data) - n50_tree = ET.ElementTree(n50_root) + n50_root = ET.fromstring(data) + n50_tree = ET.ElementTree(n50_root) - prepare_data (n50_root, n50_tree, n50_nodes, n50_ways, n50_relations) + prepare_data(n50_root, n50_tree, n50_nodes, n50_ways, n50_relations) - for element in n50_root: - if int(element.attrib['id']) > 0: - sys.exit("\t*** Please do not import existing OSM elements\n") + for element in n50_root: + if int(element.attrib["id"]) > 0: + sys.exit("\t*** Please do not import existing OSM elements\n") - message ("\tLoaded %i ways, %i relations\n" % (len(n50_ways), len(n50_relations))) + message("\tLoaded %i ways, %i relations\n" % (len(n50_ways), len(n50_relations))) # Filter elements according to given part (coastline, water, wood, landuse/other). # Found elements are returned in 'found_elements' parameter (set). -def filter_parts (part, n50_elements, found_elements): - - for element_id, element in iter(n50_elements.items()): - all_tags = element['xml'].findall('tag') - if all_tags != None: - for tag in all_tags: - key = tag.attrib['k'] - value = tag.attrib['v'] - if part == "coastline" and (key == "natural" and value == "coastline" or key == "seamark:type") or \ - part == "water" and (key == "natural" and value in ["water", "wetland", "glacier"] or key == "waterway") or \ - part == "wood" and key == "natural" and value == "wood" or \ - part == "landuse" and key in ["landuse", "leisure", "aeroway"]: - found_elements.add(element_id) - break +def filter_parts(part, n50_elements, found_elements): + + for element_id, element in iter(n50_elements.items()): + + all_tags = element["xml"].findall("tag") + if all_tags != None: + for tag in all_tags: + key = tag.attrib["k"] + value = tag.attrib["v"] + if ( + part == "coastline" + and ( + key == "natural" + and value == "coastline" + or key == "seamark:type" + ) + or part == "water" + and ( + key == "natural" + and value in ["water", "wetland", "glacier"] + or key == "waterway" + ) + or part == "wood" + and key == "natural" + and value == "wood" + or part == "landuse" + and key in ["landuse", "leisure", "aeroway"] + ): + found_elements.add(element_id) + break # Split N50 import file into parts (coastline, water, wood, landuse/other). + def split_n50(): - message ("Splitting N50 import file ...\n") + message("Splitting N50 import file ...\n") - for part in n50_parts: - message ("\t%-9s: " % part.title()) + for part in n50_parts: + message("\t%-9s: " % part.title()) - relations = set() - ways = set() - nodes = set() + relations = set() + ways = set() + nodes = set() - root = ET.Element("osm", version="0.6") - tree = ET.ElementTree(root) + root = ET.Element("osm", version="0.6") + tree = ET.ElementTree(root) - # Identify elements to include + # Identify elements to include - filter_parts(part, n50_relations, relations) - filter_parts(part, n50_ways, ways) - filter_parts(part, n50_nodes, nodes) + filter_parts(part, n50_relations, relations) + filter_parts(part, n50_ways, ways) + filter_parts(part, n50_nodes, nodes) - count_tagged_nodes = len(nodes) # Count before additional nodes are added + count_tagged_nodes = len(nodes) # Count before additional nodes are added - for relation in relations: - for member in n50_relations[ relation ]['members']: - ways.add(member) + for relation in relations: + for member in n50_relations[relation]["members"]: + ways.add(member) - # Collect additional elements for islands + # Collect additional elements for islands - if part in ["coastline", "water"]: - for relation_id, relation in iter(n50_relations.items()): - if relation['island']: + if part in ["coastline", "water"]: + for relation_id, relation in iter(n50_relations.items()): + if relation["island"]: - for member in relation['members']: - if member in ways: - relations.add(relation_id) - for island_member in relation['members']: - ways.add(island_member) - break + for member in relation["members"]: + if member in ways: + relations.add(relation_id) + for island_member in relation["members"]: + ways.add(island_member) + break - for way in ways: - for node in n50_ways[ way ]['nodes']: - nodes.add(node) + for way in ways: + for node in n50_ways[way]["nodes"]: + nodes.add(node) - # Build output tree + # Build output tree - for node in n50_root.iter("node"): - if node.attrib['id'] in nodes: - root.append(node) + for node in n50_root.iter("node"): + if node.attrib["id"] in nodes: + root.append(node) - for way in n50_root.iter("way"): - if way.attrib['id'] in ways: - root.append(way) + for way in n50_root.iter("way"): + if way.attrib["id"] in ways: + root.append(way) - for relation in n50_root.iter("relation"): - if relation.attrib['id'] in relations: - root.append(relation) + for relation in n50_root.iter("relation"): + if relation.attrib["id"] in relations: + root.append(relation) - # Generate file + # Generate file - root.set("generator", "n50merge v"+version) - root.set("upload", "false") - indent_tree(root) + root.set("generator", "n50merge v" + version) + root.set("upload", "false") + indent_tree(root) - part_filename = output_filename.replace("merged.osm", "") + part + ".osm" + part_filename = output_filename.replace("merged.osm", "") + part + ".osm" - tree.write(part_filename, encoding='utf-8', method='xml', xml_declaration=True) + tree.write(part_filename, encoding="utf-8", method="xml", xml_declaration=True) - message ("Saved %i elements to file '%s'\n" % (len(relations) + len(ways) + count_tagged_nodes, part_filename)) + message( + "Saved %i elements to file '%s'\n" + % (len(relations) + len(ways) + count_tagged_nodes, part_filename) + ) # Identify duplicate ways in existing OSM. # Note: WIP. It identifies the ways but it does not merge. + def merge_osm(): - message ("Merge OSM ways ...\n") + message("Merge OSM ways ...\n") - count = 0 - count_down = len(osm_ways) + count = 0 + count_down = len(osm_ways) - ways_list = list(osm_ways.keys()) + ways_list = list(osm_ways.keys()) - for i, way_id1 in enumerate(ways_list): - message ("\r\t%i " % count_down) - way1 = osm_ways[ way_id1 ] - count_down -= 1 + for i, way_id1 in enumerate(ways_list): + message("\r\t%i " % count_down) + way1 = osm_ways[way_id1] + count_down -= 1 - for way_id2 in ways_list[i+1:]: - way2 = osm_ways[ way_id2 ] - if (way1['coordinates'] == way2['coordinates'] or way1['coordinates'] == way2['coordinates'][::-1]) \ - and not way1['incomplete'] and not way2['incomplete']: - count += 1 - way1['xml'].append(ET.Element("tag", k="MATCH", v=way_id2)) - way1['xml'].set("action", "modify") + for way_id2 in ways_list[i + 1 :]: + way2 = osm_ways[way_id2] + if ( + ( + way1["coordinates"] == way2["coordinates"] + or way1["coordinates"] == way2["coordinates"][::-1] + ) + and not way1["incomplete"] + and not way2["incomplete"] + ): + count += 1 + way1["xml"].append(ET.Element("tag", k="MATCH", v=way_id2)) + way1["xml"].set("action", "modify") - message ("\r\tFound %i identical ways\n" % count) + message("\r\tFound %i identical ways\n" % count) # Merge tags from N50 element into OSM element, if there are no conflicts (e.g. natural=coastline + natural=wood). -def merge_tags (n50_xml, osm_xml): - n50_tags = {} - for tag in n50_xml.findall("tag"): - n50_tags[ tag.attrib['k'] ] = tag.attrib['v'] +def merge_tags(n50_xml, osm_xml): + + n50_tags = {} + for tag in n50_xml.findall("tag"): + n50_tags[tag.attrib["k"]] = tag.attrib["v"] - osm_tags = {} - for tag in osm_xml.findall("tag"): - key = tag.attrib['k'] - value = tag.attrib['v'] - osm_tags[ key ] = value - if key in n50_tags and value != n50_tags[key]: - return False + osm_tags = {} + for tag in osm_xml.findall("tag"): + key = tag.attrib["k"] + value = tag.attrib["v"] + osm_tags[key] = value + if key in n50_tags and value != n50_tags[key]: + return False - for key, value in iter(n50_tags.items()): - if key not in osm_tags or value != osm_tags[ key ]: - osm_xml.append(ET.Element("tag", k=key, v=value)) - osm_xml.set("action", "modify") + for key, value in iter(n50_tags.items()): + if key not in osm_tags or value != osm_tags[key]: + osm_xml.append(ET.Element("tag", k=key, v=value)) + osm_xml.set("action", "modify") - return True + return True # Merge N50 elements with existing OSM elements. @@ -383,66 +432,73 @@ def merge_tags (n50_xml, osm_xml): # Relations will be merged if all members have been merged. # Note that ways which differs only slightly (e.g. one more node, or a few centimeters off) will not be merged. + def merge_n50(): - message ("Merge N50 with OSM ...\n") + message("Merge N50 with OSM ...\n") - lap_time = time.time() - count_ways = 0 - count_down = len(n50_ways) - swap_nodes = {} + lap_time = time.time() + count_ways = 0 + count_down = len(n50_ways) + swap_nodes = {} - # Loop all N50 ways and match against all OSM ways + # Loop all N50 ways and match against all OSM ways - for n50_way_id in list(n50_ways.keys()): - n50_way = n50_ways[ n50_way_id ] - message ("\r\t%i " % count_down) - count_down -= 1 + for n50_way_id in list(n50_ways.keys()): + n50_way = n50_ways[n50_way_id] + message("\r\t%i " % count_down) + count_down -= 1 - for osm_way_id, osm_way in iter(osm_ways.items()): - if n50_way['coordinates'] == osm_way['coordinates'] and not n50_way['incomplete'] and not osm_way['incomplete']: + for osm_way_id, osm_way in iter(osm_ways.items()): + if ( + n50_way["coordinates"] == osm_way["coordinates"] + and not n50_way["incomplete"] + and not osm_way["incomplete"] + ): - # Check if conflicting tags, for example natural=water + natural=wood. Also update OSM tags. + # Check if conflicting tags, for example natural=water + natural=wood. Also update OSM tags. - if merge_tags(n50_way['xml'], osm_way['xml']): + if merge_tags(n50_way["xml"], osm_way["xml"]): - # Swap ref if member in relations + # Swap ref if member in relations - for parent_id in n50_way['parents']: - if parent_id in n50_relations: - parent = n50_relations[ parent_id ] - for i, member in enumerate(parent['members'][:]): - if member == n50_way_id: - parent['members'][i] = osm_way_id - member_xml = parent['xml'].find("member[@ref='%s']" % n50_way_id) - member_xml.set("ref", osm_way_id) + for parent_id in n50_way["parents"]: + if parent_id in n50_relations: + parent = n50_relations[parent_id] + for i, member in enumerate(parent["members"][:]): + if member == n50_way_id: + parent["members"][i] = osm_way_id + member_xml = parent["xml"].find( + "member[@ref='%s']" % n50_way_id + ) + member_xml.set("ref", osm_way_id) - # Mark nodes for swapping, in case nodes are used by other ways + # Mark nodes for swapping, in case nodes are used by other ways - for i in range(len(n50_way['nodes'])): - swap_nodes[ n50_way['nodes'][i] ] = osm_way['nodes'][i] + for i in range(len(n50_way["nodes"])): + swap_nodes[n50_way["nodes"][i]] = osm_way["nodes"][i] - # Remove merged way from N50 xml + # Remove merged way from N50 xml - n50_root.remove(n50_way['xml']) - del n50_ways[ n50_way_id ] - count_ways += 1 + n50_root.remove(n50_way["xml"]) + del n50_ways[n50_way_id] + count_ways += 1 - break + break - message ("\r\t \r") + message("\r\t \r") - # Swap affected nodes + # Swap affected nodes - for way_id, way in iter(n50_ways.items()): - for i, node in enumerate(way['nodes'][:]): - if node in swap_nodes: - way['nodes'][i] = swap_nodes[ node ] - node_xml = way['xml'].find("nd[@ref='%s']" % node) - node_xml.set("ref", swap_nodes[ node]) + for way_id, way in iter(n50_ways.items()): + for i, node in enumerate(way["nodes"][:]): + if node in swap_nodes: + way["nodes"][i] = swap_nodes[node] + node_xml = way["xml"].find("nd[@ref='%s']" % node) + node_xml.set("ref", swap_nodes[node]) - # Delete N50 nodes which have been replaced, unless there is a tag conflict - ''' + # Delete N50 nodes which have been replaced, unless there is a tag conflict + """ for n50_node_xml in n50_root.findall("node"): # List n50_node_id = n50_node_xml.attrib['id'] if n50_node_id in swap_nodes: @@ -453,149 +509,170 @@ def merge_n50(): del n50_nodes[ n50_node_id ] n50_root.remove(n50_node_xml) - ''' - for node in swap_nodes: - if n50_nodes[node]['xml'].find("tag") == None or not merge_tags(n50_node_xml, osm_node_xml): - n50_root.remove( n50_nodes[node]['xml'] ) - del n50_nodes[ node ] + """ + for node in swap_nodes: + if n50_nodes[node]["xml"].find("tag") == None or not merge_tags( + n50_node_xml, osm_node_xml + ): + n50_root.remove(n50_nodes[node]["xml"]) + del n50_nodes[node] - # Delete N50 relations which have been replaced by OSM relations, unless there is a tag conflict + # Delete N50 relations which have been replaced by OSM relations, unless there is a tag conflict - count_relations = 0 + count_relations = 0 - for n50_relation_id in list(n50_relations.keys()): - n50_relation = n50_relations[ n50_relation_id ] + for n50_relation_id in list(n50_relations.keys()): + n50_relation = n50_relations[n50_relation_id] - # Check if all members have been replaced, and compile set of possible OSM relations to check + # Check if all members have been replaced, and compile set of possible OSM relations to check - all_parents = set() - for member in n50_relation['members']: - if member in osm_ways: - all_parents.update(osm_ways[ member ]['parents']) - else: - all_parents = set() - break + all_parents = set() + for member in n50_relation["members"]: + if member in osm_ways: + all_parents.update(osm_ways[member]["parents"]) + else: + all_parents = set() + break - # Delete N50 relation if matching OSM relation is found, unless tag conflict + # Delete N50 relation if matching OSM relation is found, unless tag conflict - if all_parents: - for parent_id in all_parents: - if parent_id in osm_relations and set(osm_relations[ parent_id ]['members']) == set(n50_relation['members']): - if merge_tags(n50_relation['xml'], osm_relations[ parent_id ]['xml']): - n50_root.remove(n50_relation['xml']) - del n50_relations[ n50_relation_id ] - count_relations += 1 - break + if all_parents: + for parent_id in all_parents: + if parent_id in osm_relations and set( + osm_relations[parent_id]["members"] + ) == set(n50_relation["members"]): + if merge_tags(n50_relation["xml"], osm_relations[parent_id]["xml"]): + n50_root.remove(n50_relation["xml"]) + del n50_relations[n50_relation_id] + count_relations += 1 + break - message ("\r\tSwapped %i ways, %i relations\n" % (count_ways, count_relations)) - message ("\tRun time %s\n" % (timeformat(time.time() - lap_time))) + message("\r\tSwapped %i ways, %i relations\n" % (count_ways, count_relations)) + message("\tRun time %s\n" % (timeformat(time.time() - lap_time))) # Indent XML output + def indent_tree(elem, level=0): - i = "\n" + level*" " + i = "\n" + level * " " if len(elem): if not elem.text or not elem.text.strip(): elem.text = i + " " if not elem.tail or not elem.tail.strip(): elem.tail = i for elem in elem: - indent_tree(elem, level+1) + indent_tree(elem, level + 1) if not elem.tail or not elem.tail.strip(): elem.tail = i else: if level and (not elem.tail or not elem.tail.strip()): elem.tail = i + # Output merged N50/OSM tree to file. + def save_osm(): - message ("Saving file ...\n") + message("Saving file ...\n") - # Merge remaining N50 tree into OSM tree - if n50_root: - for element in n50_root: - osm_root.append(element) + # Merge remaining N50 tree into OSM tree + if n50_root: + for element in n50_root: + osm_root.append(element) - osm_root.set("generator", "n50merge v"+version) - osm_root.set("upload", "false") - indent_tree(osm_root) + osm_root.set("generator", "n50merge v" + version) + osm_root.set("upload", "false") + indent_tree(osm_root) - osm_tree.write(output_filename, encoding='utf-8', method='xml', xml_declaration=True) + osm_tree.write( + output_filename, encoding="utf-8", method="xml", xml_declaration=True + ) - message ("\tSaved to file '%s'\n" % output_filename) + message("\tSaved to file '%s'\n" % output_filename) # Main program -if __name__ == '__main__': - - start_time = time.time() - message ("\nn50merge v%s\n\n" % version) +if __name__ == "__main__": - osm_nodes = {} - osm_ways = {} - osm_relations = {} - osm_root = None - osm_tree = None + start_time = time.time() + message("\nn50merge v%s\n\n" % version) - n50_nodes = {} - n50_ways = {} - n50_relations = {} - n50_root = None - n50_tree = None + osm_nodes = {} + osm_ways = {} + osm_relations = {} + osm_root = None + osm_tree = None - debug = False # Include debug tags and unused segments - osm_merge = False # Also merge overlapping lines in OSM only + n50_nodes = {} + n50_ways = {} + n50_relations = {} + n50_root = None + n50_tree = None - # Parse parameters + debug = False # Include debug tags and unused segments + osm_merge = False # Also merge overlapping lines in OSM only - if len(sys.argv) < 2: - message ("Please provide 1) municipality, and 2) N50 filename.\n") - message ("Option: -split\n\n") - sys.exit() + # Parse parameters - # Get municipality + if len(sys.argv) < 2: + message("Please provide 1) municipality, and 2) N50 filename.\n") + message("Option: -split\n\n") + sys.exit() - municipality_query = sys.argv[1] - [municipality_id, municipality_name] = get_municipality_name(municipality_query) - if municipality_id is None: - sys.exit("Municipality '%s' not found\n" % municipality_query) - else: - message ("Municipality: %s %s\n" % (municipality_id, municipality_name)) + # Get municipality - # Determine filename and check if file exists - - if len(sys.argv) > 2 and ".osm" in sys.argv[2]: - filename = sys.argv[2] - elif len(sys.argv) > 2 and sys.argv[2] in n50_parts: - filename = "n50_%s_%s_Arealdekke_%s.osm" % (municipality_id, municipality_name, sys.argv[2]) - else: - filename = "n50_%s_%s_Arealdekke.osm" % (municipality_id, municipality_name) - - if os.path.isfile(filename) or os.path.isfile(os.path.expanduser(import_folder + filename)): - message ("N50 filename: %s\n" % filename) - output_filename = filename.replace(".osm", "") + "_merged.osm" - elif "-osm" in sys.argv: - filename = "" - output_filename = "n50_%s_%s_merged.osm" % (municipality_id, municipality_name.replace(" ", "_")) - else: - sys.exit("\t*** File '%s' not found\n\n" % filename) - - message ("\n") + municipality_query = sys.argv[1] + [municipality_id, municipality_name] = get_municipality_name(municipality_query) + if municipality_id is None: + sys.exit("Municipality '%s' not found\n" % municipality_query) + else: + message("Municipality: %s %s\n" % (municipality_id, municipality_name)) + + # Determine filename and check if file exists + + if len(sys.argv) > 2 and ".osm" in sys.argv[2]: + filename = sys.argv[2] + elif len(sys.argv) > 2 and sys.argv[2] in n50_parts: + filename = "n50_%s_%s_Arealdekke_%s.osm" % ( + municipality_id, + municipality_name, + sys.argv[2], + ) + else: + filename = "n50_%s_%s_Arealdekke.osm" % (municipality_id, municipality_name) + + if os.path.isfile(filename) or os.path.isfile( + os.path.expanduser(import_folder + filename) + ): + message("N50 filename: %s\n" % filename) + output_filename = filename.replace(".osm", "") + "_merged.osm" + elif "-osm" in sys.argv: + filename = "" + output_filename = "n50_%s_%s_merged.osm" % ( + municipality_id, + municipality_name.replace(" ", "_"), + ) + else: + sys.exit("\t*** File '%s' not found\n\n" % filename) - # Process data + message("\n") - load_n50() + # Process data - if "-split" in sys.argv: - split_n50() - else: - load_osm() - merge_n50() - save_osm() + load_n50() - duration = time.time() - start_time - message ("\tTotal run time %s (%i ways per second)\n\n" % (timeformat(duration), int(len(n50_ways) / duration))) + if "-split" in sys.argv: + split_n50() + else: + load_osm() + merge_n50() + save_osm() + + duration = time.time() - start_time + message( + "\tTotal run time %s (%i ways per second)\n\n" + % (timeformat(duration), int(len(n50_ways) / duration)) + ) diff --git a/n50osm.py b/n50osm.py index 19e349f..6aa4969 100644 --- a/n50osm.py +++ b/n50osm.py @@ -25,238 +25,273 @@ lake_ele_size = 2000 # Minimum square meters for fetching elevation -data_categories = ["AdministrativeOmrader", "Arealdekke", "BygningerOgAnlegg", "Hoyde", "Restriksjonsomrader", "Samferdsel", "Stedsnavn"] +data_categories = [ + "AdministrativeOmrader", + "Arealdekke", + "BygningerOgAnlegg", + "Hoyde", + "Restriksjonsomrader", + "Samferdsel", + "Stedsnavn", +] avoid_objects = [ # Object types to exclude from output - 'ÅpentOmråde', 'Tregruppe', # Arealdekke - 'GangSykkelveg', 'VegSenterlinje', 'Vegsperring', # Samferdsel - 'Forsenkningskurve', 'Hjelpekurve', 'Høydekurve', # Hoyde - 'PresentasjonTekst' # Stedsnavn + "ÅpentOmråde", + "Tregruppe", # Arealdekke + "GangSykkelveg", + "VegSenterlinje", + "Vegsperring", # Samferdsel + "Forsenkningskurve", + "Hjelpekurve", + "Høydekurve", # Hoyde + "PresentasjonTekst", # Stedsnavn ] -auxiliary_objects = ['Arealbrukgrense', 'Dataavgrensning', 'FiktivDelelinje', \ - 'InnsjøElvSperre', 'InnsjøInnsjøSperre', 'ElvBekkKant', 'Havflate', 'Innsjøkant', 'InnsjøkantRegulert', 'FerskvannTørrfallkant'] +auxiliary_objects = [ + "Arealbrukgrense", + "Dataavgrensning", + "FiktivDelelinje", + "InnsjøElvSperre", + "InnsjøInnsjøSperre", + "ElvBekkKant", + "Havflate", + "Innsjøkant", + "InnsjøkantRegulert", + "FerskvannTørrfallkant", +] avoid_tags = [ # N50 properties to exclude from output (unless debug) - 'oppdateringsdato', 'datafangstdato', - 'målemetode', 'nøyaktighet' + "oppdateringsdato", + "datafangstdato", + "målemetode", + "nøyaktighet", ] osm_tags = { - # Arealdekke - - 'Alpinbakke': { 'landuse': 'winter_sports', 'piste:type': 'downhill', 'area': 'yes' }, - 'BymessigBebyggelse': { 'landuse': 'retail' }, #'landuse': 'residential', 'residential': 'urban' }, - 'DyrketMark': { 'landuse': 'farmland' }, - 'ElvBekk': { 'waterway': 'stream' }, - 'FerskvannTørrfall': { 'waterway': 'riverbank', 'intermittent': 'yes' }, - 'Foss': { 'waterway': 'waterfall' }, - 'Golfbane': { 'leisure': 'golf_course' }, - 'Gravplass': { 'landuse': 'cemetery' }, - 'HavElvSperre': { 'natural': 'coastline' }, - 'HavInnsjøSperre': { 'natural': 'coastline' }, - 'Hyttefelt': { 'landuse': 'residential', 'residential': 'cabin' }, - 'Industriområde': { 'landuse': 'industrial' }, - 'Innsjø': { 'natural': 'water' }, - 'InnsjøRegulert': { 'natural': 'water', 'water': 'reservoir' }, - 'Kystkontur': { 'natural': 'coastline' }, - 'Lufthavn': { 'aeroway': 'aerodrome' }, - 'Myr': { 'natural': 'wetland', 'wetland': 'bog' }, - 'Park': { 'leisure': 'park' }, - 'Rullebane': { 'aeroway': 'runway' }, - 'Skjær': { 'seamark:type': 'rock' }, - 'Skog': { 'natural': 'wood' }, - 'Skytefelt': { 'leisure': 'pitch', 'sport': 'shooting' }, - 'SnøIsbre': { 'natural': 'glacier' }, - 'SportIdrettPlass': { 'leisure': 'pitch' }, - 'Steinbrudd': { 'landuse': 'quarry'}, - 'Steintipp': { 'landuse': 'landfill' }, - 'Tettbebyggelse': { 'landuse': 'residential' }, - - # Samferdsel - - 'Barmarksløype': { 'highway': 'track' }, - 'Traktorveg': { 'highway': 'track' }, - 'Sti': { 'highway': 'path' }, - - # Hoyde - - 'Terrengpunkt': { 'natural': 'hill' }, - 'TrigonometriskPunkt': { 'natural': 'hill' }, - - # Restriksjonsomrader - - 'Naturvernområde': { 'boundary': 'protected_area' }, - 'Allmenning': { 'boundary': 'protected_area', 'protect_class': '27'}, # Public land - - # BygningerOgAnlegg - - 'Bygning': { 'building': 'yes' }, - 'Campingplass': { 'tourism': 'camp_site' }, - 'Dam': { 'waterway': 'dam' }, - 'Flytebrygge': { 'man_made': 'pier', 'floating': 'yes' }, - 'Gruve': { 'man_made': 'adit' }, # Could be shaft - 'Hoppbakke': { 'piste:type': 'ski_jump' }, - 'KaiBrygge': { 'man_made': 'quay' }, - 'Ledning': { 'power': 'line' }, - 'LuftledningLH': { 'power': 'line' }, - 'Lysløype': { 'highway': 'track', 'lit': 'yes', 'trailblazed': 'yes' }, - 'MastTele': { 'man_made': 'mast', 'tower:type': 'communication' }, - 'Molo': { 'man_made': 'breakwater' }, - 'Navigasjonsinstallasjon': { 'man_made': 'lighthouse' }, # Only lighthouses, it seems - 'Parkeringsområde': { 'amenity': 'parking' }, - 'Pir': { 'man_made': 'pier' }, - 'Reingjerde': { 'barrier': 'fence' }, - 'Rørgate': { 'man_made': 'pipeline' }, # Also "tømmerrenne" - 'Skitrekk': { 'aerialway': 'drag_lift' }, # Could be other aerialway values - 'Skytebaneinnretning': { 'leisure': 'pitch', 'sport': 'shooting' }, - 'Tank': { 'man_made': 'tank' }, - 'Taubane': { 'aerialway': 'cable_car' }, # Could be other aerial values, e.g. gondola, goods - 'Tårn': { 'man_made': 'tower' }, # Any massive or substantial tower - 'Vindkraftverk': { 'power': 'generator', 'generator:source': 'wind', 'generator:type': 'horizontal_axis' } - + # Arealdekke + "Alpinbakke": {"landuse": "winter_sports", "piste:type": "downhill", "area": "yes"}, + "BymessigBebyggelse": { + "landuse": "retail" + }, #'landuse': 'residential', 'residential': 'urban' }, + "DyrketMark": {"landuse": "farmland"}, + "ElvBekk": {"waterway": "stream"}, + "FerskvannTørrfall": {"waterway": "riverbank", "intermittent": "yes"}, + "Foss": {"waterway": "waterfall"}, + "Golfbane": {"leisure": "golf_course"}, + "Gravplass": {"landuse": "cemetery"}, + "HavElvSperre": {"natural": "coastline"}, + "HavInnsjøSperre": {"natural": "coastline"}, + "Hyttefelt": {"landuse": "residential", "residential": "cabin"}, + "Industriområde": {"landuse": "industrial"}, + "Innsjø": {"natural": "water"}, + "InnsjøRegulert": {"natural": "water", "water": "reservoir"}, + "Kystkontur": {"natural": "coastline"}, + "Lufthavn": {"aeroway": "aerodrome"}, + "Myr": {"natural": "wetland", "wetland": "bog"}, + "Park": {"leisure": "park"}, + "Rullebane": {"aeroway": "runway"}, + "Skjær": {"seamark:type": "rock"}, + "Skog": {"natural": "wood"}, + "Skytefelt": {"leisure": "pitch", "sport": "shooting"}, + "SnøIsbre": {"natural": "glacier"}, + "SportIdrettPlass": {"leisure": "pitch"}, + "Steinbrudd": {"landuse": "quarry"}, + "Steintipp": {"landuse": "landfill"}, + "Tettbebyggelse": {"landuse": "residential"}, + # Samferdsel + "Barmarksløype": {"highway": "track"}, + "Traktorveg": {"highway": "track"}, + "Sti": {"highway": "path"}, + # Hoyde + "Terrengpunkt": {"natural": "hill"}, + "TrigonometriskPunkt": {"natural": "hill"}, + # Restriksjonsomrader + "Naturvernområde": {"boundary": "protected_area"}, + "Allmenning": {"boundary": "protected_area", "protect_class": "27"}, # Public land + # BygningerOgAnlegg + "Bygning": {"building": "yes"}, + "Campingplass": {"tourism": "camp_site"}, + "Dam": {"waterway": "dam"}, + "Flytebrygge": {"man_made": "pier", "floating": "yes"}, + "Gruve": {"man_made": "adit"}, # Could be shaft + "Hoppbakke": {"piste:type": "ski_jump"}, + "KaiBrygge": {"man_made": "quay"}, + "Ledning": {"power": "line"}, + "LuftledningLH": {"power": "line"}, + "Lysløype": {"highway": "track", "lit": "yes", "trailblazed": "yes"}, + "MastTele": {"man_made": "mast", "tower:type": "communication"}, + "Molo": {"man_made": "breakwater"}, + "Navigasjonsinstallasjon": {"man_made": "lighthouse"}, # Only lighthouses, it seems + "Parkeringsområde": {"amenity": "parking"}, + "Pir": {"man_made": "pier"}, + "Reingjerde": {"barrier": "fence"}, + "Rørgate": {"man_made": "pipeline"}, # Also "tømmerrenne" + "Skitrekk": {"aerialway": "drag_lift"}, # Could be other aerialway values + "Skytebaneinnretning": {"leisure": "pitch", "sport": "shooting"}, + "Tank": {"man_made": "tank"}, + "Taubane": { + "aerialway": "cable_car" + }, # Could be other aerial values, e.g. gondola, goods + "Tårn": {"man_made": "tower"}, # Any massive or substantial tower + "Vindkraftverk": { + "power": "generator", + "generator:source": "wind", + "generator:type": "horizontal_axis", + }, } - # OSM tagging; first special cases + def tag_object(feature_type, geometry_type, properties, feature): - tags = {} - missing_tags = set() - - # First special object cases - - if feature_type == "ElvBekk": - if geometry_type == "område": - tags['waterway'] = "riverbank" - elif "vannBredde" in properties and properties['vannBredde'] > "2": # >3 meter - tags['waterway'] = "river" - else: - tags['waterway'] = "stream" - - elif feature_type == "Skytefelt" and data_category == "Restriksjonsomrader": # Eception to Arealdekke - tags['landuse'] = "military" - - elif feature_type == "Bygning": - if "bygningstype" in properties: - if properties['bygningstype'] == "956": # Turisthytte - if "betjeningsgrad" in properties: - if properties['betjeningsgrad'] == "B": # Betjent - tags['tourism'] = "alpine_hut" - elif properties['betjeningsgrad'] == "S": # Selvbetjent - tags['tourism'] = "wilderness_hut" - elif properties['betjeningsgrad'] in ["U", "D", "R"]: # Ubetjent, dagstur, rastebu - tags['amenity'] = "shelter" - tags['shelter_type'] = "basic_hut" - else: - tags['amenity'] = "shelter" - tags['shelter_type'] = "lean_to" - if "hytteeier" in properties: - if properties['hytteeier'] == "1": - tags['operator'] = "DNT" - elif properties['hytteeier'] == "3": - tags['operator'] = "Fjellstyre" - elif properties['hytteeier'] == "4": - tags['operator'] = "Statskog" - - elif properties['bygningstype'] in building_tags: - for key, value in iter(building_tags[ properties['bygningstype'] ].items()): - if geometry_type == "område" or key != "building": # or len(building_tags[ properties['bygningstype'] ]) > 1: - tags[ key ] = value - - if geometry_type != "posisjon" and "building" not in tags: - tags['building'] = "yes" # No building tag for single nodes - - elif feature_type == "Lufthavn": - if "lufthavntype" in properties and properties['lufthavntype'] == "H": - tags['aeroway'] = "heliport" - else: - tags['aeroway'] = "aerodrome" - if "trafikktype" in properties: - if properties['trafikktype'] == "I": - tags['aeroway:type'] = "international" - elif properties['trafikktype'] == "N": - tags['aeroway:type'] = "regional" - elif properties['trafikktype'] == "A": - tags['aeroway:type'] = "airfield" - if "iataKode" in properties and properties['iataKode'] != "XXX": - tags['iata'] = properties['iataKode'] - if "icaoKode" in properties and properties['icaoKode'] != "XXXX": - tags['icao'] = properties['icaoKode'] - - elif geometry_type == "område" and feature_type == "SportIdrettPlass": - if len(feature['coordinates']) > 1: - tags['leisure'] = "track" - tags['area'] = "yes" - else: - tags['leisure'] = "pitch" - - # Then conversion dict - - elif feature_type in osm_tags: - tags.update( osm_tags[feature_type] ) - - # Collect set of remaining object types not handled - - elif feature_type not in auxiliary_objects: - missing_tags.add(feature_type) - - # Additional tagging based on object properties from GML - - if "høyde" in properties: - tags['ele'] = properties['høyde'] - if "lavesteRegulerteVannstand" in properties: - tags['ele:min'] = properties['lavesteRegulerteVannstand'] - if "vatnLøpenummer" in properties: - tags['ref:nve:vann'] = properties['vatnLøpenummer'] - - if "navn" in properties: - tags['name'] = properties['navn'] - if "fulltekst" in properties: - tags['name'] = properties['fulltekst'] - if "stedsnummer" in properties: - tags['ssr:stedsnr'] = properties['stedsnummer'] -# if "skriftkode" in properties: -# tags['SKRIFTKODE'] = properties['skriftkode'] - - if "merking" in properties and properties['merking'] == "JA": - tags['trailblazed'] = "yes" - - if "verneform" in properties: - if properties['verneform'] in ["NP", "NPS"]: - tags['boundary'] = "national_park" - elif properties['verneform'] in ["LVO", "NM"]: - tags['boundary'] = "protected_area" - else: - tags['leisure'] = "nature_reserve" - - if "lengde" in properties and feature_type == "Hoppbakke": - tags['ref'] = "K" + properties['lengde'] - - return (tags, missing_tags) + tags = {} + missing_tags = set() + + # First special object cases + + if feature_type == "ElvBekk": + if geometry_type == "område": + tags["waterway"] = "riverbank" + elif "vannBredde" in properties and properties["vannBredde"] > "2": # >3 meter + tags["waterway"] = "river" + else: + tags["waterway"] = "stream" + + elif ( + feature_type == "Skytefelt" and data_category == "Restriksjonsomrader" + ): # Eception to Arealdekke + tags["landuse"] = "military" + + elif feature_type == "Bygning": + if "bygningstype" in properties: + if properties["bygningstype"] == "956": # Turisthytte + if "betjeningsgrad" in properties: + if properties["betjeningsgrad"] == "B": # Betjent + tags["tourism"] = "alpine_hut" + elif properties["betjeningsgrad"] == "S": # Selvbetjent + tags["tourism"] = "wilderness_hut" + elif properties["betjeningsgrad"] in [ + "U", + "D", + "R", + ]: # Ubetjent, dagstur, rastebu + tags["amenity"] = "shelter" + tags["shelter_type"] = "basic_hut" + else: + tags["amenity"] = "shelter" + tags["shelter_type"] = "lean_to" + if "hytteeier" in properties: + if properties["hytteeier"] == "1": + tags["operator"] = "DNT" + elif properties["hytteeier"] == "3": + tags["operator"] = "Fjellstyre" + elif properties["hytteeier"] == "4": + tags["operator"] = "Statskog" + + elif properties["bygningstype"] in building_tags: + for key, value in iter( + building_tags[properties["bygningstype"]].items() + ): + if ( + geometry_type == "område" or key != "building" + ): # or len(building_tags[ properties['bygningstype'] ]) > 1: + tags[key] = value + + if geometry_type != "posisjon" and "building" not in tags: + tags["building"] = "yes" # No building tag for single nodes + + elif feature_type == "Lufthavn": + if "lufthavntype" in properties and properties["lufthavntype"] == "H": + tags["aeroway"] = "heliport" + else: + tags["aeroway"] = "aerodrome" + if "trafikktype" in properties: + if properties["trafikktype"] == "I": + tags["aeroway:type"] = "international" + elif properties["trafikktype"] == "N": + tags["aeroway:type"] = "regional" + elif properties["trafikktype"] == "A": + tags["aeroway:type"] = "airfield" + if "iataKode" in properties and properties["iataKode"] != "XXX": + tags["iata"] = properties["iataKode"] + if "icaoKode" in properties and properties["icaoKode"] != "XXXX": + tags["icao"] = properties["icaoKode"] + + elif geometry_type == "område" and feature_type == "SportIdrettPlass": + if len(feature["coordinates"]) > 1: + tags["leisure"] = "track" + tags["area"] = "yes" + else: + tags["leisure"] = "pitch" + + # Then conversion dict + + elif feature_type in osm_tags: + tags.update(osm_tags[feature_type]) + + # Collect set of remaining object types not handled + + elif feature_type not in auxiliary_objects: + missing_tags.add(feature_type) + + # Additional tagging based on object properties from GML + + if "høyde" in properties: + tags["ele"] = properties["høyde"] + if "lavesteRegulerteVannstand" in properties: + tags["ele:min"] = properties["lavesteRegulerteVannstand"] + if "vatnLøpenummer" in properties: + tags["ref:nve:vann"] = properties["vatnLøpenummer"] + + if "navn" in properties: + tags["name"] = properties["navn"] + if "fulltekst" in properties: + tags["name"] = properties["fulltekst"] + if "stedsnummer" in properties: + tags["ssr:stedsnr"] = properties["stedsnummer"] + # if "skriftkode" in properties: + # tags['SKRIFTKODE'] = properties['skriftkode'] + + if "merking" in properties and properties["merking"] == "JA": + tags["trailblazed"] = "yes" + + if "verneform" in properties: + if properties["verneform"] in ["NP", "NPS"]: + tags["boundary"] = "national_park" + elif properties["verneform"] in ["LVO", "NM"]: + tags["boundary"] = "protected_area" + else: + tags["leisure"] = "nature_reserve" + + if "lengde" in properties and feature_type == "Hoppbakke": + tags["ref"] = "K" + properties["lengde"] + + return (tags, missing_tags) # Output message -def message (output_text): - sys.stdout.write (output_text) - sys.stdout.flush() +def message(output_text): + + sys.stdout.write(output_text) + sys.stdout.flush() # Format time -def timeformat (sec): - if sec > 3600: - return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) - elif sec > 60: - return "%i:%02i minutes" % (sec / 60, sec % 60) - else: - return "%i seconds" % sec +def timeformat(sec): + + if sec > 3600: + return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) + elif sec > 60: + return "%i:%02i minutes" % (sec / 60, sec % 60) + else: + return "%i seconds" % sec # Calculate coordinate area of polygon in square meters @@ -265,609 +300,693 @@ def timeformat (sec): # > 0: Counter-clockwise # = 0: Polygon not closed -def polygon_area (polygon): - if polygon[0] == polygon[-1]: - lat_dist = math.pi * 6371009.0 / 180.0 +def polygon_area(polygon): + + if polygon[0] == polygon[-1]: + lat_dist = math.pi * 6371009.0 / 180.0 - coord = [] - for node in polygon: - y = node[1] * lat_dist - x = node[0] * lat_dist * math.cos(math.radians(node[1])) - coord.append((x,y)) + coord = [] + for node in polygon: + y = node[1] * lat_dist + x = node[0] * lat_dist * math.cos(math.radians(node[1])) + coord.append((x, y)) - area = 0.0 - for i in range(len(coord) - 1): - area += (coord[i+1][1] - coord[i][1]) * (coord[i+1][0] + coord[i][0]) # (x2-x1)(y2+y1) + area = 0.0 + for i in range(len(coord) - 1): + area += (coord[i + 1][1] - coord[i][1]) * ( + coord[i + 1][0] + coord[i][0] + ) # (x2-x1)(y2+y1) - return int(area / 2.0) - else: - return 0 + return int(area / 2.0) + else: + return 0 # Calculate coordinate area of multipolygon, i.e. excluding inner polygons -def multipolygon_area (multipolygon): - if type(multipolygon) is list and len(multipolygon) > 0 and type(multipolygon[0]) is list and \ - multipolygon[0][0] == multipolygon[0][-1]: +def multipolygon_area(multipolygon): - area = polygon_area(multipolygon[0]) - for patch in multipolygon[1:]: - inner_area = polygon_area(patch) - if inner_area: - area -= inner_area - else: - return None - return area + if ( + type(multipolygon) is list + and len(multipolygon) > 0 + and type(multipolygon[0]) is list + and multipolygon[0][0] == multipolygon[0][-1] + ): - else: - return None + area = polygon_area(multipolygon[0]) + for patch in multipolygon[1:]: + inner_area = polygon_area(patch) + if inner_area: + area -= inner_area + else: + return None + return area + + else: + return None # Calculate centroid of polygon # Source: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon -def polygon_centroid (polygon): - if polygon[0] == polygon[-1]: - x = 0 - y = 0 - det = 0 +def polygon_centroid(polygon): - for i in range(len(polygon) - 1): - d = polygon[i][0] * polygon[i+1][1] - polygon[i+1][0] * polygon[i][1] - det += d - x += (polygon[i][0] + polygon[i+1][0]) * d # (x1 + x2) (x1*y2 - x2*y1) - y += (polygon[i][1] + polygon[i+1][1]) * d # (y1 + y2) (x1*y2 - x2*y1) + if polygon[0] == polygon[-1]: + x = 0 + y = 0 + det = 0 - return (x / (3.0 * det), y / (3.0 * det) ) + for i in range(len(polygon) - 1): + d = polygon[i][0] * polygon[i + 1][1] - polygon[i + 1][0] * polygon[i][1] + det += d + x += (polygon[i][0] + polygon[i + 1][0]) * d # (x1 + x2) (x1*y2 - x2*y1) + y += (polygon[i][1] + polygon[i + 1][1]) * d # (y1 + y2) (x1*y2 - x2*y1) - else: - return None + return (x / (3.0 * det), y / (3.0 * det)) + + else: + return None # Tests whether point (x,y) is inside a polygon # Ray tracing method -def inside_polygon (point, polygon): - if polygon[0] == polygon[-1]: - x, y = point - n = len(polygon) - inside = False +def inside_polygon(point, polygon): + + if polygon[0] == polygon[-1]: + x, y = point + n = len(polygon) + inside = False - p1x, p1y = polygon[0] - for i in range(n): - p2x, p2y = polygon[i] - if y > min(p1y, p2y): - if y <= max(p1y, p2y): - if x <= max(p1x, p2x): - if p1y != p2y: - xints = (y-p1y) * (p2x-p1x) / (p2y-p1y) + p1x - if p1x == p2x or x <= xints: - inside = not inside - p1x, p1y = p2x, p2y + p1x, p1y = polygon[0] + for i in range(n): + p2x, p2y = polygon[i] + if y > min(p1y, p2y): + if y <= max(p1y, p2y): + if x <= max(p1x, p2x): + if p1y != p2y: + xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x + if p1x == p2x or x <= xints: + inside = not inside + p1x, p1y = p2x, p2y - return inside + return inside - else: - return None + else: + return None # Test whether point (x,y) is inside a multipolygon, i.e. not inside inner polygons -def inside_multipolygon (point, multipolygon): - if type(multipolygon) is list and len(multipolygon) > 0 and type(multipolygon[0]) is list and \ - multipolygon[0][0] == multipolygon[0][-1]: +def inside_multipolygon(point, multipolygon): - inside = inside_polygon(point, multipolygon[0]) - if inside: - for patch in multipolygon[1:]: - inside = (inside and not inside_polygon(point, patch)) - if not inside: - break + if ( + type(multipolygon) is list + and len(multipolygon) > 0 + and type(multipolygon[0]) is list + and multipolygon[0][0] == multipolygon[0][-1] + ): + + inside = inside_polygon(point, multipolygon[0]) + if inside: + for patch in multipolygon[1:]: + inside = inside and not inside_polygon(point, patch) + if not inside: + break - return inside + return inside - else: - return None + else: + return None # Calculate new node with given distance offset in meters # Works over short distances -def coordinate_offset (node, distance): - m = (1 / ((math.pi / 180.0) * 6378137.0)) # Degrees per meter +def coordinate_offset(node, distance): - latitude = node[1] + (distance * m) - longitude = node[0] + (distance * m) / math.cos( math.radians(node[1]) ) + m = 1 / ((math.pi / 180.0) * 6378137.0) # Degrees per meter - return (longitude, latitude) + latitude = node[1] + (distance * m) + longitude = node[0] + (distance * m) / math.cos(math.radians(node[1])) + + return (longitude, latitude) # Identify bounds of line coordinates # Returns lower left and upper right corners of square bounds + extra perimeter (in meters) + def get_bbox(coordinates, perimeter): - if type(coordinates) is tuple: - patch = [ coordinates ] - elif type(coordinates[0]) is tuple: - patch = coordinates - else: - patch = coordinates[0] - - min_node = list(patch[0]) - max_node = copy.deepcopy(min_node) + if type(coordinates) is tuple: + patch = [coordinates] + elif type(coordinates[0]) is tuple: + patch = coordinates + else: + patch = coordinates[0] + + min_node = list(patch[0]) + max_node = copy.deepcopy(min_node) - for node in patch: - for i in [0,1]: - min_node[i] = min(min_node[i], node[i]) - max_node[i] = max(max_node[i], node[i]) + for node in patch: + for i in [0, 1]: + min_node[i] = min(min_node[i], node[i]) + max_node[i] = max(max_node[i], node[i]) - if perimeter > 0: - min_node = coordinate_offset(min_node, - perimeter) - max_node = coordinate_offset(max_node, + perimeter) + if perimeter > 0: + min_node = coordinate_offset(min_node, -perimeter) + max_node = coordinate_offset(max_node, +perimeter) - return [min_node, max_node] + return [min_node, max_node] # Create feature with one point -def create_point (node, gml_id, note): - if debug: - entry = { - 'object': 'Debug', - 'type': 'Point', - 'coordinates': node, - 'members': [], - 'tags': {}, - 'extras': { - 'objekttype': 'Debug', - 'note': note - } - } - if gml_id: - entry['gml_id'] = gml_id +def create_point(node, gml_id, note): + + if debug: + entry = { + "object": "Debug", + "type": "Point", + "coordinates": node, + "members": [], + "tags": {}, + "extras": {"objekttype": "Debug", "note": note}, + } + if gml_id: + entry["gml_id"] = gml_id - features.append(entry) + features.append(entry) # Get list of coordinates from GML # Each point is a tuple of (lon, lat), corresponding to GeoJSON format x,y -def parse_coordinates (coord_text): - - global gml_id - - parse_count = 0 - split_coord = coord_text.split(" ") - coordinates = [] - for i in range(0, len(split_coord) - 1, 2): - x = float(split_coord[i]) - y = float(split_coord[i+1]) - [lat, lon] = utm.UtmToLatLon (x, y, 33, "N") - node = ( round(lon, coordinate_decimals), round(lat, coordinate_decimals) ) - parse_count += 1 - if not coordinates or node != coordinates[-1] or json_output: - coordinates.append(node) - else: -# message ("\t*** DELETED DUPLICATE NODE: %s %s\n" % (node, gml_id)) - create_point(node, gml_id, "deleted duplicate") - - # Remove single outlayer node - - if not json_output: - i = 0 - while i < len(coordinates): - if i > 1 and coordinates[i] == coordinates[i - 2]: -# message ("\t*** DELETED ARTEFACT NODE: %s %s\n" % (coordinates[ i - 1 ], gml_id)) - create_point(copy.deepcopy(coordinates[ i - 1 ]), gml_id, "deleted artefact") - coordinates.pop(i) - coordinates.pop(i - 1) - i -= 1 - else: - i += 1 - - if len(coordinates) > 2 and coordinates[0] == coordinates[-1] and coordinates[1] == coordinates[-2] : -# message ("\t*** DELETED ARTEFACT NODE: %s %s\n" % (coordinates[ 0 ], gml_id)) - coordinates = coordinates[ 1: -1 ] - - if len(coordinates) == 1 and parse_count > 1: - message ("\t*** SHORT WAY %s\n" % gml_id) - - return coordinates + +def parse_coordinates(coord_text): + + global gml_id + + parse_count = 0 + split_coord = coord_text.split(" ") + coordinates = [] + for i in range(0, len(split_coord) - 1, 2): + x = float(split_coord[i]) + y = float(split_coord[i + 1]) + [lat, lon] = utm.UtmToLatLon(x, y, 33, "N") + node = (round(lon, coordinate_decimals), round(lat, coordinate_decimals)) + parse_count += 1 + if not coordinates or node != coordinates[-1] or json_output: + coordinates.append(node) + else: + # message ("\t*** DELETED DUPLICATE NODE: %s %s\n" % (node, gml_id)) + create_point(node, gml_id, "deleted duplicate") + + # Remove single outlayer node + + if not json_output: + i = 0 + while i < len(coordinates): + if i > 1 and coordinates[i] == coordinates[i - 2]: + # message ("\t*** DELETED ARTEFACT NODE: %s %s\n" % (coordinates[ i - 1 ], gml_id)) + create_point( + copy.deepcopy(coordinates[i - 1]), gml_id, "deleted artefact" + ) + coordinates.pop(i) + coordinates.pop(i - 1) + i -= 1 + else: + i += 1 + + if ( + len(coordinates) > 2 + and coordinates[0] == coordinates[-1] + and coordinates[1] == coordinates[-2] + ): + # message ("\t*** DELETED ARTEFACT NODE: %s %s\n" % (coordinates[ 0 ], gml_id)) + coordinates = coordinates[1:-1] + + if len(coordinates) == 1 and parse_count > 1: + message("\t*** SHORT WAY %s\n" % gml_id) + + return coordinates # Get name or id of municipality from GeoNorge api -def get_municipality_name (query): - - if query.isdigit(): - url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query - else: - url = "https://ws.geonorge.no/kommuneinfo/v1/sok?knavn=" + urllib.parse.quote(query) - - request = urllib.request.Request(url, headers=header) - - try: - file = urllib.request.urlopen(request) - except urllib.error.HTTPError as e: - if e.code == 404: # Not found - sys.exit("\tMunicipality '%s' not found\n\n" % query) - else: - raise - - if query.isdigit(): - result = json.load(file) - file.close() - municipality_name = result['kommunenavnNorsk'] - return (query, municipality_name) - - else: - result = json.load(file) - file.close() - if result['antallTreff'] == 1: - municipality_id = result['kommuner'][0]['kommunenummer'] - municipality_name = result['kommuner'][0]['kommunenavnNorsk'] - return (municipality_id, municipality_name) - else: - municipalities = [] - for municipality in result['kommuner']: - municipalities.append(municipality['kommunenummer'] + " " + municipalities['kommunenavnNorsk']) - sys.exit("\tMore than one municipality found: %s\n\n" % ", ".join(municipalities)) + +def get_municipality_name(query): + + if query.isdigit(): + url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query + else: + url = "https://ws.geonorge.no/kommuneinfo/v1/sok?knavn=" + urllib.parse.quote( + query + ) + + request = urllib.request.Request(url, headers=header) + + try: + file = urllib.request.urlopen(request) + except urllib.error.HTTPError as e: + if e.code == 404: # Not found + sys.exit("\tMunicipality '%s' not found\n\n" % query) + else: + raise + + if query.isdigit(): + result = json.load(file) + file.close() + municipality_name = result["kommunenavnNorsk"] + return (query, municipality_name) + + else: + result = json.load(file) + file.close() + if result["antallTreff"] == 1: + municipality_id = result["kommuner"][0]["kommunenummer"] + municipality_name = result["kommuner"][0]["kommunenavnNorsk"] + return (municipality_id, municipality_name) + else: + municipalities = [] + for municipality in result["kommuner"]: + municipalities.append( + municipality["kommunenummer"] + + " " + + municipalities["kommunenavnNorsk"] + ) + sys.exit( + "\tMore than one municipality found: %s\n\n" % ", ".join(municipalities) + ) # Load conversion CSV table for tagging building types # Format in CSV: "key=value + key=value + ..." + def load_building_types(): -# file = open("building_types.csv") - url = "https://raw.githubusercontent.com/NKAmapper/building2osm/main/building_types.csv" - request = urllib.request.Request(url, headers=header) - file = urllib.request.urlopen(request) + # file = open("building_types.csv") + url = "https://raw.githubusercontent.com/NKAmapper/building2osm/main/building_types.csv" + request = urllib.request.Request(url, headers=header) + file = urllib.request.urlopen(request) - building_csv = csv.DictReader(TextIOWrapper(file, "utf-8"), fieldnames=["id", "name", "building_tag", "extra_tag"], delimiter=";") - next(building_csv) + building_csv = csv.DictReader( + TextIOWrapper(file, "utf-8"), + fieldnames=["id", "name", "building_tag", "extra_tag"], + delimiter=";", + ) + next(building_csv) - for row in building_csv: - tag_string = (row['building_tag'] + "+" + row['extra_tag']).strip().strip("+") + for row in building_csv: + tag_string = (row["building_tag"] + "+" + row["extra_tag"]).strip().strip("+") - if tag_string: - osm_tag = {} - tag_list = tag_string.replace(" ","").split("+") + if tag_string: + osm_tag = {} + tag_list = tag_string.replace(" ", "").split("+") - for tag_part in tag_list: - tag_split = tag_part.split("=") - osm_tag[ tag_split[0] ] = tag_split[1] + for tag_part in tag_list: + tag_split = tag_part.split("=") + osm_tag[tag_split[0]] = tag_split[1] - building_tags[ row['id'] ] = osm_tag + building_tags[row["id"]] = osm_tag - file.close() + file.close() # Compute length based on coordinates (not in meters) -def simple_length (coord): - length = 0 - for i in range(len(coord) - 2): - length += (coord[i+1][0] - coord[i][0])**2 + ((coord[i+1][1] - coord[i][1])**2) * 0.5 +def simple_length(coord): + + length = 0 + for i in range(len(coord) - 2): + length += (coord[i + 1][0] - coord[i][0]) ** 2 + ( + (coord[i + 1][1] - coord[i][1]) ** 2 + ) * 0.5 - return length + return length # Split patch if self-intersecting or touching polygon -def split_patch (coordinates): - for i in range(1, len(coordinates) - 1): - first = coordinates.index(coordinates[i]) - if first < i: - result1 = split_patch( coordinates[ : first ] + coordinates[ i: ] ) - result2 = split_patch( coordinates[ first : i + 1 ]) -# message ("\t*** SPLIT SELF-INTERSECTING/TOUCHING POLYGON: %s\n" % str(coordinates[i])) - - if simple_length(result1[0]) > simple_length(result2[0]): - result1.extend(result2) - return result1 - else: - result2.extend(result1) - return result2 +def split_patch(coordinates): - return [ coordinates ] + for i in range(1, len(coordinates) - 1): + first = coordinates.index(coordinates[i]) + if first < i: + result1 = split_patch(coordinates[:first] + coordinates[i:]) + result2 = split_patch(coordinates[first : i + 1]) + # message ("\t*** SPLIT SELF-INTERSECTING/TOUCHING POLYGON: %s\n" % str(coordinates[i])) + + if simple_length(result1[0]) > simple_length(result2[0]): + result1.extend(result2) + return result1 + else: + result2.extend(result1) + return result2 + + return [coordinates] # Get all app properties from nested XML; recursive search -def get_property(top, ns_app): - properties = {} - if ns_app in top.tag: - tag = top.tag[ len(ns_app)+2 : ] - value = top.text - value = value.strip() - if value: - properties[tag] = value +def get_property(top, ns_app): + + properties = {} + if ns_app in top.tag: + tag = top.tag[len(ns_app) + 2 :] + value = top.text + value = value.strip() + if value: + properties[tag] = value - for child in top: - properties.update(get_property(child, ns_app)) + for child in top: + properties.update(get_property(child, ns_app)) - return properties + return properties # Load N50 topo data from Kartverket -def load_n50_data (municipality_id, municipality_name, data_category): - global gml_id +def load_n50_data(municipality_id, municipality_name, data_category): - lap = time.time() + global gml_id - message ("\nLoad N50 data from Kartverket...\n") + lap = time.time() - source_date = ["9", "0"] # First and last source date ("datafangstdato") - update_date = ["9", "0"] # First and last update date ("oppdateringsdato") - stream_count = 0 - missing_tags = set() - - # Load latest N50 file for municipality from Kartverket - - filename = "Basisdata_%s_%s_25833_N50Kartdata_GML" % (municipality_id, municipality_name) - filename = filename.replace("Æ","E").replace("Ø","O").replace("Å","A")\ - .replace("æ","e").replace("ø","o").replace("å","a").replace(" ", "_") - message ("\tLoading file '%s'\n" % filename) - - request = urllib.request.Request("https://nedlasting.geonorge.no/geonorge/Basisdata/N50Kartdata/GML/" + filename + ".zip", headers=header) - file_in = urllib.request.urlopen(request) - zip_file = zipfile.ZipFile(BytesIO(file_in.read())) - -# for file_entry in zip_file.namelist(): -# message ("\t%s\n" % file_entry) - - filename2 = filename.replace("Kartdata", data_category) # For example "Arealdekke" - file = zip_file.open(filename2 + ".gml") - - tree = ET.parse(file) - file.close() - file_in.close() - root = tree.getroot() - - ns_gml = 'http://www.opengis.net/gml/3.2' - ns_app = 'http://skjema.geonorge.no/SOSI/produktspesifikasjon/N50/20170401' - - ns = { - 'gml': ns_gml, - 'app': ns_app - } - - # Loop features, parse and load into data structure and tag - - message ("\tParsing...\n") - - for feature in root: - if "featureMember" in feature.tag: - feature_type = feature[0].tag[ len(ns_app)+2 : ] - geometry_type = None - gml_id = feature[0].attrib["{%s}id" % ns_gml] - - if feature_type not in object_count: - object_count[feature_type] = 0 - object_count[feature_type] += 1 - - if feature_type in avoid_objects and not json_output: - continue - - entry = { - 'object': feature_type, - 'type': None, - 'gml_id': gml_id, - 'coordinates': [], - 'members': [], - 'tags': {}, - 'extras': { - 'objekttype': feature_type - } - } - - properties = {} # Attributes provided from GML - - for app in feature[0]: - - # Get app properties/attributes - - tag = app.tag[ len(ns_app)+2 : ] - if tag in ['posisjon', 'grense', 'område', 'senterlinje', 'geometri']: - geometry_type = tag - entry['extras']['geometri'] = geometry_type - else: - properties.update(get_property(app, ns_app)) - - # Get geometry/coordinates - - for geo in app: - # point - if geo.tag == "{%s}Point" % ns_gml: - entry['type'] = "Point" - entry['coordinates'] = parse_coordinates(geo[0].text)[0] - - # LineString - elif geo.tag == "{%s}LineString" % ns_gml: - if geometry_type != "geometri": - entry['type'] = "LineString" - entry['coordinates'] = parse_coordinates(geo[0].text) - else: - entry['type'] = "Point" - entry['coordinates'] = parse_coordinates(geo[0].text)[0] - - # Curve, stored as LineString - elif geo.tag == "{%s}Curve" % ns_gml: - entry['type'] = "LineString" - entry['extras']['type2'] = "curve" - entry['coordinates'] = [] - for patch in geo[0]: - coordinates = parse_coordinates(patch[0].text) - if entry['coordinates']: - entry['coordinates'] + coordinates[1:] - else: - entry['coordinates'] = coordinates # First patch - - # (Multi)Polygon - elif geo.tag == "{%s}Surface" % ns_gml: - entry['type'] = "Polygon" - entry['coordinates'] = [] # List of patches - - for patch in geo[0][0]: - role = patch.tag[ len(ns_gml)+2 : ] - - if patch[0].tag[ len(ns_gml)+2 : ] == "Ring": - coordinates = parse_coordinates(patch[0][0][0][0][0][0].text) # Ring->Curve->LineStringSegment - else: - coordinates = parse_coordinates(patch[0][0].text) # LinearRing - - if len(coordinates) > 1: - if len(coordinates) < 3: - message("\t*** POLYGON PATCH TOO SHORT: %s %s\n" % (role, gml_id)) - elif coordinates[0] != coordinates[-1]: - message("\t*** POLYGON PATCH NOT CLOSED: %s %s\n" % (role, gml_id)) - - if json_output: - entry['coordinates'].append( coordinates ) - else: - # Check for intersecting/touching polygon - entry['coordinates'].extend( split_patch( coordinates ) ) - else: - message("\t*** EMPTY POLYGON PATCH: %s %s\n" % (role, gml_id)) - - elif ns_gml in geo.tag: - message ("\t*** UNKNOWN GEOMETRY: %s\n" % geo.tag) - - # Store tags - - [ tags, new_missing_tags ] = tag_object(feature_type, geometry_type, properties, entry) - entry['tags'].update(tags) - entry['extras'].update(properties) - missing_tags.update(new_missing_tags) - - if n50_tags and not debug: - for key, value in iter(properties.items()): - if key not in avoid_tags: - entry['tags'][ "N50_" + key ] = value - - # Add to relevant list - - if not (entry['type'] == "LineString" and len(entry['coordinates']) <= 1): - if geometry_type == "grense": - if feature_type in ["Kystkontur", "HavElvSperre", "HavInnsjøSperre"]: - entry['used'] = 1 - else: - entry['used'] = 0 - - segments.append(entry) - elif not (entry['type'] == "Point" and not entry['tags']) or debug: # Omit untagged single points - features.append(entry) - else: - message ("\t*** SEGMENT TOO SHORT: %s\n" % gml_id) - - # Update min/max dates for information - - if "datafangstdato" in properties: - if properties['datafangstdato'] < source_date[0] and properties['datafangstdato'] > "1801": - source_date[0] = properties['datafangstdato'] - if properties['datafangstdato'] > source_date[1]: - source_date[1] = properties['datafangstdato'] - - if "oppdateringsdato" in properties: - if properties['oppdateringsdato'] < update_date[0] and properties ['oppdateringsdato'] > "1801": - update_date[0] = properties['oppdateringsdato'] - if properties['oppdateringsdato'] > update_date[1]: - update_date[1] = properties['oppdateringsdato'] - - if feature_type == "ElvBekk": - stream_count += 1 - - message("\tObjects loaded:\n") - for object_type in sorted(object_count): - if object_type not in auxiliary_objects: - message("\t\t%i\t%s\n" % (object_count[object_type], object_type)) - - if missing_tags: - message ("\tNot tagged: %s\n" % (", ".join(missing_tags.difference("Havflate")))) - if stream_count > 0: - message ("\t%i streams\n" % stream_count) - - message ("\tSource dates: %s - %s\n" % (source_date[0], source_date[1])) - message ("\tUpdate dates: %s - %s\n" % (update_date[0], update_date[1])) - message ("\t%i feature objects, %i segments\n" % (len(features), len(segments))) - message ("\tRun time %s\n" % (timeformat(time.time() - lap))) + message("\nLoad N50 data from Kartverket...\n") + + source_date = ["9", "0"] # First and last source date ("datafangstdato") + update_date = ["9", "0"] # First and last update date ("oppdateringsdato") + stream_count = 0 + missing_tags = set() + + # Load latest N50 file for municipality from Kartverket + + filename = "Basisdata_%s_%s_25833_N50Kartdata_GML" % ( + municipality_id, + municipality_name, + ) + filename = ( + filename.replace("Æ", "E") + .replace("Ø", "O") + .replace("Å", "A") + .replace("æ", "e") + .replace("ø", "o") + .replace("å", "a") + .replace(" ", "_") + ) + message("\tLoading file '%s'\n" % filename) + + request = urllib.request.Request( + "https://nedlasting.geonorge.no/geonorge/Basisdata/N50Kartdata/GML/" + + filename + + ".zip", + headers=header, + ) + file_in = urllib.request.urlopen(request) + zip_file = zipfile.ZipFile(BytesIO(file_in.read())) + + # for file_entry in zip_file.namelist(): + # message ("\t%s\n" % file_entry) + + filename2 = filename.replace("Kartdata", data_category) # For example "Arealdekke" + file = zip_file.open(filename2 + ".gml") + + tree = ET.parse(file) + file.close() + file_in.close() + root = tree.getroot() + + ns_gml = "http://www.opengis.net/gml/3.2" + ns_app = "http://skjema.geonorge.no/SOSI/produktspesifikasjon/N50/20170401" + + ns = {"gml": ns_gml, "app": ns_app} + + # Loop features, parse and load into data structure and tag + + message("\tParsing...\n") + + for feature in root: + if "featureMember" in feature.tag: + feature_type = feature[0].tag[len(ns_app) + 2 :] + geometry_type = None + gml_id = feature[0].attrib["{%s}id" % ns_gml] + + if feature_type not in object_count: + object_count[feature_type] = 0 + object_count[feature_type] += 1 + + if feature_type in avoid_objects and not json_output: + continue + + entry = { + "object": feature_type, + "type": None, + "gml_id": gml_id, + "coordinates": [], + "members": [], + "tags": {}, + "extras": {"objekttype": feature_type}, + } + + properties = {} # Attributes provided from GML + + for app in feature[0]: + + # Get app properties/attributes + + tag = app.tag[len(ns_app) + 2 :] + if tag in ["posisjon", "grense", "område", "senterlinje", "geometri"]: + geometry_type = tag + entry["extras"]["geometri"] = geometry_type + else: + properties.update(get_property(app, ns_app)) + + # Get geometry/coordinates + + for geo in app: + # point + if geo.tag == "{%s}Point" % ns_gml: + entry["type"] = "Point" + entry["coordinates"] = parse_coordinates(geo[0].text)[0] + + # LineString + elif geo.tag == "{%s}LineString" % ns_gml: + if geometry_type != "geometri": + entry["type"] = "LineString" + entry["coordinates"] = parse_coordinates(geo[0].text) + else: + entry["type"] = "Point" + entry["coordinates"] = parse_coordinates(geo[0].text)[0] + + # Curve, stored as LineString + elif geo.tag == "{%s}Curve" % ns_gml: + entry["type"] = "LineString" + entry["extras"]["type2"] = "curve" + entry["coordinates"] = [] + for patch in geo[0]: + coordinates = parse_coordinates(patch[0].text) + if entry["coordinates"]: + entry["coordinates"] + coordinates[1:] + else: + entry["coordinates"] = coordinates # First patch + + # (Multi)Polygon + elif geo.tag == "{%s}Surface" % ns_gml: + entry["type"] = "Polygon" + entry["coordinates"] = [] # List of patches + + for patch in geo[0][0]: + role = patch.tag[len(ns_gml) + 2 :] + + if patch[0].tag[len(ns_gml) + 2 :] == "Ring": + coordinates = parse_coordinates( + patch[0][0][0][0][0][0].text + ) # Ring->Curve->LineStringSegment + else: + coordinates = parse_coordinates( + patch[0][0].text + ) # LinearRing + + if len(coordinates) > 1: + if len(coordinates) < 3: + message( + "\t*** POLYGON PATCH TOO SHORT: %s %s\n" + % (role, gml_id) + ) + elif coordinates[0] != coordinates[-1]: + message( + "\t*** POLYGON PATCH NOT CLOSED: %s %s\n" + % (role, gml_id) + ) + + if json_output: + entry["coordinates"].append(coordinates) + else: + # Check for intersecting/touching polygon + entry["coordinates"].extend( + split_patch(coordinates) + ) + else: + message( + "\t*** EMPTY POLYGON PATCH: %s %s\n" + % (role, gml_id) + ) + + elif ns_gml in geo.tag: + message("\t*** UNKNOWN GEOMETRY: %s\n" % geo.tag) + + # Store tags + + [tags, new_missing_tags] = tag_object( + feature_type, geometry_type, properties, entry + ) + entry["tags"].update(tags) + entry["extras"].update(properties) + missing_tags.update(new_missing_tags) + + if n50_tags and not debug: + for key, value in iter(properties.items()): + if key not in avoid_tags: + entry["tags"]["N50_" + key] = value + + # Add to relevant list + + if not (entry["type"] == "LineString" and len(entry["coordinates"]) <= 1): + if geometry_type == "grense": + if feature_type in [ + "Kystkontur", + "HavElvSperre", + "HavInnsjøSperre", + ]: + entry["used"] = 1 + else: + entry["used"] = 0 + + segments.append(entry) + elif ( + not (entry["type"] == "Point" and not entry["tags"]) or debug + ): # Omit untagged single points + features.append(entry) + else: + message("\t*** SEGMENT TOO SHORT: %s\n" % gml_id) + + # Update min/max dates for information + + if "datafangstdato" in properties: + if ( + properties["datafangstdato"] < source_date[0] + and properties["datafangstdato"] > "1801" + ): + source_date[0] = properties["datafangstdato"] + if properties["datafangstdato"] > source_date[1]: + source_date[1] = properties["datafangstdato"] + + if "oppdateringsdato" in properties: + if ( + properties["oppdateringsdato"] < update_date[0] + and properties["oppdateringsdato"] > "1801" + ): + update_date[0] = properties["oppdateringsdato"] + if properties["oppdateringsdato"] > update_date[1]: + update_date[1] = properties["oppdateringsdato"] + + if feature_type == "ElvBekk": + stream_count += 1 + + message("\tObjects loaded:\n") + for object_type in sorted(object_count): + if object_type not in auxiliary_objects: + message("\t\t%i\t%s\n" % (object_count[object_type], object_type)) + + if missing_tags: + message("\tNot tagged: %s\n" % (", ".join(missing_tags.difference("Havflate")))) + if stream_count > 0: + message("\t%i streams\n" % stream_count) + + message("\tSource dates: %s - %s\n" % (source_date[0], source_date[1])) + message("\tUpdate dates: %s - %s\n" % (update_date[0], update_date[1])) + message("\t%i feature objects, %i segments\n" % (len(features), len(segments))) + message("\tRun time %s\n" % (timeformat(time.time() - lap))) # Create missing KantUtsnitt segments to get complete polygons + def create_border_segments(patch, members, gml_id, match): - # First create list of existing conncetions between coordinates i and i+1 of patch - - connection = [] - for i in range(len(patch)-1): - connection.append(False) - - for member in members: - segment = segments[member] - n0 = patch.index(segment['coordinates'][0]) - - for node in segment['coordinates'][1:]: - n1 = patch.index(node) - if abs(n1 - n0) == 1: - connection[ min(n0, n1) ] = True - elif abs(n1 - n0) == len(patch) - 2: - connection[ len(patch) - 2 ] = True - else: - message ("\t*** MISSING BORDER SEGMENT: Node %i - %i of %i/%i nodes %s\n" % (n0, n1, len(patch), match, gml_id)) - - n0 = n1 - - # Then create new segments for the parts of patch which have no segments - - start_index = 0 - while start_index < len(patch) - 2: - - while start_index < len(patch) - 1 and connection[start_index]: - start_index += 1 - - end_index = start_index - while end_index < len(patch) - 1 and not connection[end_index]: - end_index += 1 - - if end_index > start_index: - entry = { - 'object': 'KantUtsnitt', - 'type': 'LineString', - 'coordinates': patch[ start_index : end_index + 1], - 'members': [], - 'tags': { }, - 'extras': { - 'objekttype': 'KantUtsnitt' - }, - 'used': 1 - } - [ entry['min_bbox'], entry['max_bbox'] ] = get_bbox(entry['coordinates'], 0) - segments.append(entry) - members.append( len(segments) - 1 ) - start_index = end_index + # First create list of existing conncetions between coordinates i and i+1 of patch + + connection = [] + for i in range(len(patch) - 1): + connection.append(False) + + for member in members: + segment = segments[member] + n0 = patch.index(segment["coordinates"][0]) + + for node in segment["coordinates"][1:]: + n1 = patch.index(node) + if abs(n1 - n0) == 1: + connection[min(n0, n1)] = True + elif abs(n1 - n0) == len(patch) - 2: + connection[len(patch) - 2] = True + else: + message( + "\t*** MISSING BORDER SEGMENT: Node %i - %i of %i/%i nodes %s\n" + % (n0, n1, len(patch), match, gml_id) + ) + + n0 = n1 + + # Then create new segments for the parts of patch which have no segments + + start_index = 0 + while start_index < len(patch) - 2: + + while start_index < len(patch) - 1 and connection[start_index]: + start_index += 1 + + end_index = start_index + while end_index < len(patch) - 1 and not connection[end_index]: + end_index += 1 + + if end_index > start_index: + entry = { + "object": "KantUtsnitt", + "type": "LineString", + "coordinates": patch[start_index : end_index + 1], + "members": [], + "tags": {}, + "extras": {"objekttype": "KantUtsnitt"}, + "used": 1, + } + [entry["min_bbox"], entry["max_bbox"]] = get_bbox(entry["coordinates"], 0) + segments.append(entry) + members.append(len(segments) - 1) + start_index = end_index # Function for sorting member segments of polygon relation # Index 1 used to avoid equal 0/-1 positions + def segment_position(segment_index, patch): - return patch.index(segments[segment_index]['coordinates'][1]) + return patch.index(segments[segment_index]["coordinates"][1]) # Fix data structure: @@ -876,225 +995,296 @@ def segment_position(segment_index, patch): # - Order members of multipolygons # - Crate missing segments along municipality border -def split_polygons(): - - message ("Decompose polygons into segments...\n") - - # Create bbox for segments and line features - - for segment in segments: - [ segment['min_bbox'], segment['max_bbox'] ] = get_bbox(segment['coordinates'], 0) - - # Loop all polygons and patches - - lap = time.time() - split_count = 0 - - for feature in features: - - if feature['type'] == "Polygon": - matching_polygon = [] - - for patch in feature['coordinates']: - matching_segments = [] - matched_nodes = 0 - patch_set = set(patch) - - # Try matching with segments within the polygon's bbox - - [ patch_min_bbox, patch_max_bbox ] = get_bbox(patch, 0) - for i, segment in enumerate(segments): - - if patch_min_bbox[0] <= segment['max_bbox'][0] and patch_max_bbox[0] >= segment['min_bbox'][0] and \ - patch_min_bbox[1] <= segment['max_bbox'][1] and patch_max_bbox[1] >= segment['min_bbox'][1] and \ - set(segment['coordinates']) <= patch_set: - - # Note: If patch is a closed way, segment may wrap start/end of patch - - if len(segment['coordinates']) == 2: - node1 = patch.index(segment['coordinates'][0]) - node2 = patch.index(segment['coordinates'][1]) - if not(abs(node1-node2) == 1 or patch[0] == patch[-1] and abs(node1-node2) == len(patch) - 2): - continue - - matching_segments.append(i) - matched_nodes += len(segment['coordinates']) - 1 - - # Correct direction of coastline, lakes and riverways - - if feature['object'] == "Havflate" and \ - segment['object'] in ["Kystkontur", "HavElvSperre", "HavInnsjøSperre"] or \ - feature['object'] in ["Innsjø", "InnsjøRegulert", "ElvBekk", 'FerskvannTørrfall'] and \ - segment['object'] in ["Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant", "KantUtsnitt"]: - - segment['used'] += 1 - node1 = patch.index(segment['coordinates'][0]) - node2 = patch.index(segment['coordinates'][1]) - - if not (node1 + 1 == node2 or patch[0] == patch[-1] and node1 == len(patch) - 2 and node2 == 0): - segment['coordinates'].reverse() - segment['extras']['reversert'] = "yes" - - elif feature['object'] != "Havflate": - segment['used'] += 1 - - if matched_nodes == len(patch) - 1: - break - - if matching_segments: - # Use leftover nodes to create missing border segments - if matched_nodes < len(patch) - 1 and feature['object'] != "Havflate": - create_border_segments(patch, matching_segments, feature['gml_id'], matched_nodes) - - # Sort relation members for better presentation - matching_segments.sort(key=lambda segment_index: segment_position(segment_index, patch)) - matching_polygon.append(matching_segments) - split_count += len(matching_segments) - 1 - else: - message ("\t*** NO MATCH: %s\n" % (feature['gml_id'])) - feature['extras']['segmentering'] = "no" - - feature['members'] = matching_polygon +def split_polygons(): - message ("\t%i splits\n" % split_count) - message ("\tRun time %s\n" % (timeformat(time.time() - lap))) + message("Decompose polygons into segments...\n") + + # Create bbox for segments and line features + + for segment in segments: + [segment["min_bbox"], segment["max_bbox"]] = get_bbox(segment["coordinates"], 0) + + # Loop all polygons and patches + + lap = time.time() + split_count = 0 + + for feature in features: + + if feature["type"] == "Polygon": + matching_polygon = [] + + for patch in feature["coordinates"]: + matching_segments = [] + matched_nodes = 0 + patch_set = set(patch) + + # Try matching with segments within the polygon's bbox + + [patch_min_bbox, patch_max_bbox] = get_bbox(patch, 0) + + for i, segment in enumerate(segments): + + if ( + patch_min_bbox[0] <= segment["max_bbox"][0] + and patch_max_bbox[0] >= segment["min_bbox"][0] + and patch_min_bbox[1] <= segment["max_bbox"][1] + and patch_max_bbox[1] >= segment["min_bbox"][1] + and set(segment["coordinates"]) <= patch_set + ): + + # Note: If patch is a closed way, segment may wrap start/end of patch + + if len(segment["coordinates"]) == 2: + node1 = patch.index(segment["coordinates"][0]) + node2 = patch.index(segment["coordinates"][1]) + if not ( + abs(node1 - node2) == 1 + or patch[0] == patch[-1] + and abs(node1 - node2) == len(patch) - 2 + ): + continue + + matching_segments.append(i) + matched_nodes += len(segment["coordinates"]) - 1 + + # Correct direction of coastline, lakes and riverways + + if ( + feature["object"] == "Havflate" + and segment["object"] + in ["Kystkontur", "HavElvSperre", "HavInnsjøSperre"] + or feature["object"] + in [ + "Innsjø", + "InnsjøRegulert", + "ElvBekk", + "FerskvannTørrfall", + ] + and segment["object"] + in [ + "Innsjøkant", + "InnsjøkantRegulert", + "ElvBekkKant", + "KantUtsnitt", + ] + ): + + segment["used"] += 1 + node1 = patch.index(segment["coordinates"][0]) + node2 = patch.index(segment["coordinates"][1]) + + if not ( + node1 + 1 == node2 + or patch[0] == patch[-1] + and node1 == len(patch) - 2 + and node2 == 0 + ): + segment["coordinates"].reverse() + segment["extras"]["reversert"] = "yes" + + elif feature["object"] != "Havflate": + segment["used"] += 1 + + if matched_nodes == len(patch) - 1: + break + + if matching_segments: + # Use leftover nodes to create missing border segments + if ( + matched_nodes < len(patch) - 1 + and feature["object"] != "Havflate" + ): + create_border_segments( + patch, matching_segments, feature["gml_id"], matched_nodes + ) + + # Sort relation members for better presentation + matching_segments.sort( + key=lambda segment_index: segment_position(segment_index, patch) + ) + matching_polygon.append(matching_segments) + split_count += len(matching_segments) - 1 + else: + message("\t*** NO MATCH: %s\n" % (feature["gml_id"])) + feature["extras"]["segmentering"] = "no" + + feature["members"] = matching_polygon + + message("\t%i splits\n" % split_count) + message("\tRun time %s\n" % (timeformat(time.time() - lap))) # Fix coastline -def find_islands(): - - message ("Find islands...\n") - - lap = time.time() - island_count = 0 - - # Part 1: Identify islands described by inner parts of lakes and sea - - # First build list of other candidate relations - - candidates = [] - for feature in features: - if len(feature['members']) == 1 and len(feature['members'][0]) > 1 and \ - feature['object'] not in ["Innsjø", "InnsjøRegulert", "ElvBekk", "Havflate", "FerskvannTørrfall"]: - found = True - for member in feature['members'][0]: - if segments[ member ]['object'] not in ['Kystkontur', "Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant"]: - found = False - break - if found: - candidates.append(feature) - - # Loop all inner objects of multipolygon lakes and sea - - for feature in features: - if feature['object'] in ["Innsjø", "InnsjøRegulert", "ElvBekk", "Havflate", "FerskvannTørrfall"]: - for i in range(1, len(feature['members'])): - - # Do not use patch with intermittent edge - - found = True - for member in feature['members'][i]: - if segments[ member ]['object'] == "FerskvannTørrfallkant": - found = False - break - if not found: - continue - # Determine island type based on area - - area = polygon_area(feature['coordinates'][i]) - - if abs(area) > island_size: - island_type = "island" - else: - island_type = "islet" - -# if area < 0: -# message ("\t*** ISLAND WITH REVERSE COASTLINE: %s\n" % feature['gml_id']) - - # Tag closed way if possible - - if len(feature['members'][i]) == 1: - segment = segments[ feature['members'][i][0] ] - segment['tags']['place'] = island_type - segment['extras']['areal'] = str(int(abs(area))) - - else: - # Serach for already existing relation - - found = False - for feature2 in candidates: - if set(feature['members'][i]) == set(feature2['members'][0]): - feature2['tags']['place'] = island_type - feature2['extras']['areal'] = str(int(abs(area))) - break - - # Else create new polygon - - if not found: - entry = { - 'object': "Øy", - 'type': 'Polygon', - 'coordinates': copy.deepcopy(feature['coordinates'][i]), - 'members': [ copy.deepcopy(feature['members'][i]) ], - 'tags': { 'place': island_type }, - 'extras': { 'areal': str(int(abs(area))) } - } - - features.append(entry) - - island_count += 1 - - - # Part 2: Identify remaining islands - # Note: This section will also find all lakes, so check coastline direction - - # First build unordered list of segment coastline - - coastlines = [] -# for i, segment in enumerate(segments): -# if segment['object'] in ["Kystkontur", "Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant"] and i not in used_segments: -# coastlines.append(segment) - - for feature in features: - if feature['object'] in ["Havflate", "Innsjø", "InnsjøRegulert", "ElvBekk", "FerskvannTørrfall"] and \ - len(feature['members']) > 0: - for member in feature['members'][0]: # Only outer patch - segment = segments[ member ] - if segment['object'] in ["HavElvSperre", "HavInnsjøSperre", "InnsjøInnsjøSperre", \ - "InnsjøElvSperre", "FerskvannTørrfallkant", "FiktivDelelinje"]: - for member2 in feature['members'][0]: - segment2 = segments[ member2 ] - if segment2['object'] in ["Kystkontur", "Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant"]: - coastlines.append(segment2) - break - - # Merge coastline segments until exhausted - - while coastlines: - segment = coastlines[0] - island = [ coastlines[0] ] - coastlines.pop(0) - first_node = segment['coordinates'][0] - last_node = segment['coordinates'][-1] - - # Build coastline/island forward - - found = True - while found and first_node != last_node: - found = False - for segment in coastlines[:]: - if segment['coordinates'][0] == last_node: - last_node = segment['coordinates'][-1] - island.append(segment) - coastlines.remove(segment) - found = True - break +def find_islands(): - # Build coastline/island backward - ''' + message("Find islands...\n") + + lap = time.time() + island_count = 0 + + # Part 1: Identify islands described by inner parts of lakes and sea + + # First build list of other candidate relations + + candidates = [] + for feature in features: + if ( + len(feature["members"]) == 1 + and len(feature["members"][0]) > 1 + and feature["object"] + not in [ + "Innsjø", + "InnsjøRegulert", + "ElvBekk", + "Havflate", + "FerskvannTørrfall", + ] + ): + found = True + for member in feature["members"][0]: + if segments[member]["object"] not in [ + "Kystkontur", + "Innsjøkant", + "InnsjøkantRegulert", + "ElvBekkKant", + ]: + found = False + break + if found: + candidates.append(feature) + + # Loop all inner objects of multipolygon lakes and sea + + for feature in features: + if feature["object"] in [ + "Innsjø", + "InnsjøRegulert", + "ElvBekk", + "Havflate", + "FerskvannTørrfall", + ]: + for i in range(1, len(feature["members"])): + + # Do not use patch with intermittent edge + + found = True + for member in feature["members"][i]: + if segments[member]["object"] == "FerskvannTørrfallkant": + found = False + break + if not found: + continue + + # Determine island type based on area + + area = polygon_area(feature["coordinates"][i]) + + if abs(area) > island_size: + island_type = "island" + else: + island_type = "islet" + + # if area < 0: + # message ("\t*** ISLAND WITH REVERSE COASTLINE: %s\n" % feature['gml_id']) + + # Tag closed way if possible + + if len(feature["members"][i]) == 1: + segment = segments[feature["members"][i][0]] + segment["tags"]["place"] = island_type + segment["extras"]["areal"] = str(int(abs(area))) + + else: + # Serach for already existing relation + + found = False + for feature2 in candidates: + if set(feature["members"][i]) == set(feature2["members"][0]): + feature2["tags"]["place"] = island_type + feature2["extras"]["areal"] = str(int(abs(area))) + break + + # Else create new polygon + + if not found: + entry = { + "object": "Øy", + "type": "Polygon", + "coordinates": copy.deepcopy(feature["coordinates"][i]), + "members": [copy.deepcopy(feature["members"][i])], + "tags": {"place": island_type}, + "extras": {"areal": str(int(abs(area)))}, + } + + features.append(entry) + + island_count += 1 + + # Part 2: Identify remaining islands + # Note: This section will also find all lakes, so check coastline direction + + # First build unordered list of segment coastline + + coastlines = [] + # for i, segment in enumerate(segments): + # if segment['object'] in ["Kystkontur", "Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant"] and i not in used_segments: + # coastlines.append(segment) + + for feature in features: + if ( + feature["object"] + in ["Havflate", "Innsjø", "InnsjøRegulert", "ElvBekk", "FerskvannTørrfall"] + and len(feature["members"]) > 0 + ): + for member in feature["members"][0]: # Only outer patch + segment = segments[member] + if segment["object"] in [ + "HavElvSperre", + "HavInnsjøSperre", + "InnsjøInnsjøSperre", + "InnsjøElvSperre", + "FerskvannTørrfallkant", + "FiktivDelelinje", + ]: + for member2 in feature["members"][0]: + segment2 = segments[member2] + if segment2["object"] in [ + "Kystkontur", + "Innsjøkant", + "InnsjøkantRegulert", + "ElvBekkKant", + ]: + coastlines.append(segment2) + break + + # Merge coastline segments until exhausted + + while coastlines: + segment = coastlines[0] + island = [coastlines[0]] + coastlines.pop(0) + first_node = segment["coordinates"][0] + last_node = segment["coordinates"][-1] + + # Build coastline/island forward + + found = True + while found and first_node != last_node: + found = False + for segment in coastlines[:]: + if segment["coordinates"][0] == last_node: + last_node = segment["coordinates"][-1] + island.append(segment) + coastlines.remove(segment) + found = True + break + + # Build coastline/island backward + """ found = True while found and first_node != last_node: found = False @@ -1105,151 +1295,189 @@ def find_islands(): coastlines.remove(segment) found = True break - ''' + """ - # Add island to features list if closed chain of ways + # Add island to features list if closed chain of ways - if first_node == last_node: + if first_node == last_node: - members = [] - coordinates = [ first_node ] - for segment in island: - members.append(segments.index(segment)) - coordinates += segment['coordinates'][1:] + members = [] + coordinates = [first_node] + for segment in island: + members.append(segments.index(segment)) + coordinates += segment["coordinates"][1:] - area = polygon_area(coordinates) - if area < 0: - continue + area = polygon_area(coordinates) + if area < 0: + continue - if abs(area) > island_size: - island_type = "island" - else: - island_type = "islet" + if abs(area) > island_size: + island_type = "island" + else: + island_type = "islet" - # Reuse existing relation if possible + # Reuse existing relation if possible - found = False - for feature in candidates: - if set(members) == set(feature['members'][0]): - feature['tags']['place'] = island_type - feature['extras']['areal'] = str(int(abs(area))) - island_count += 1 - found = True - break - - # Else create new relation for island + found = False + for feature in candidates: + if set(members) == set(feature["members"][0]): + feature["tags"]["place"] = island_type + feature["extras"]["areal"] = str(int(abs(area))) + island_count += 1 + found = True + break - if not found: - entry = { - 'object': "Øy", - 'type': 'Polygon', - 'coordinates': [ coordinates ], - 'members': [ members ], - 'tags': copy.deepcopy(island[0]['tags']), - 'extras': copy.deepcopy(island[0]['extras']) - } + # Else create new relation for island - entry['tags']['place'] = island_type - entry['tags'].pop("natural", None) # Remove natural=coastline (already on segments) - entry['extras']['areal'] = str(int(abs(area))) + if not found: + entry = { + "object": "Øy", + "type": "Polygon", + "coordinates": [coordinates], + "members": [members], + "tags": copy.deepcopy(island[0]["tags"]), + "extras": copy.deepcopy(island[0]["extras"]), + } - features.append(entry) - island_count += 1 + entry["tags"]["place"] = island_type + entry["tags"].pop( + "natural", None + ) # Remove natural=coastline (already on segments) + entry["extras"]["areal"] = str(int(abs(area))) - # Remove Havflate objects, which are not used going forward + features.append(entry) + island_count += 1 - for feature in features[:]: - if feature['object'] == "Havflate": - features.remove(feature) + # Remove Havflate objects, which are not used going forward - message ("\t%i islands\n" % island_count) - message ("\tRun time %s\n" % (timeformat(time.time() - lap))) + for feature in features[:]: + if feature["object"] == "Havflate": + features.remove(feature) + message("\t%i islands\n" % island_count) + message("\tRun time %s\n" % (timeformat(time.time() - lap))) # Identify common intersection nodes between lines (e.g. streams) -def match_nodes(): - - global delete_count - - message ("Identify intersections...\n") - - lap = time.time() - node_count = len(nodes) - delete_count = 0 - - # Create set of common nodes for segment intersections - - for segment in segments: - if segment['used'] > 0: - nodes.add(segment['coordinates'][0]) - nodes.add(segment['coordinates'][-1]) - - for feature in features: - if feature['type'] == "LineString": - nodes.add(feature['coordinates'][0]) - nodes.add(feature['coordinates'][-1]) - - if not no_node: - - # Create bbox for segments and line features + create temporary list of LineString features - - for segment in segments: - [ segment['min_bbox'], segment['max_bbox'] ] = get_bbox(segment['coordinates'], 0) - - # Loop streams to identify intersections with segments - - for feature in features: - if feature['type'] == "LineString" and feature['object'] == "ElvBekk": - [ feature['min_bbox'], feature['max_bbox'] ] = get_bbox(feature['coordinates'], 0) - - for segment in segments: - if (segment['used'] > 0 or debug) and \ - (feature['min_bbox'][0] <= segment['max_bbox'][0] and feature['max_bbox'][0] >= segment['min_bbox'][0] and \ - feature['min_bbox'][1] <= segment['max_bbox'][1] and feature['max_bbox'][1] >= segment['min_bbox'][1]): - - intersections = set(feature['coordinates']).intersection(set(segment['coordinates'])) - if len(intersections) == 0: - continue - - for node in intersections: - index1 = feature['coordinates'].index(node) - index2 = segment['coordinates'].index(node) - - # First check if stream node may be removed og slightly relocated - - if index1 not in [0, len(feature['coordinates']) - 1] and node not in nodes: - if feature['coordinates'][ index1 - 1] not in intersections and \ - feature['coordinates'][ index1 + 1] not in intersections: - feature['coordinates'].pop(index1) - else: - lon, lat = node - offset = 10 ** (- coordinate_decimals + 1) # Last lat/lon decimal digit - feature['coordinates'][index1] = ( lon + 4 * offset, lat + 2 * offset ) - # Note: New node used in next test here - - # Then check if segment node may also be removed - - if index2 not in [0, len(segment['coordinates']) - 1]: - if segment['coordinates'][ index2 - 1] not in intersections and \ - segment['coordinates'][ index2 + 1] not in intersections: - segment['coordinates'].pop(index2) - - # Else create new common node with "water" segments (or reuse existing common node) - elif segment['object'] in ["Kystkontur", "Innsjøkant", "InnsjøkantRegulert", "ElvBekkKant"]: - nodes.add( node ) - - # Loop auxiliary lines and simplify geometry - - for segment in segments: - if segment['used'] > 0 and segment['object'] == "FiktivDelelinje": - delete_count += len(segment['coordinates']) - 2 - segment['coordinates'] = [ segment['coordinates'][0], segment['coordinates'][-1] ] +def match_nodes(): - message ("\t%i common nodes, %i nodes removed from streams and auxiliary lines\n" % (len(nodes) - node_count, delete_count)) - message ("\tRun time %s\n" % (timeformat(time.time() - lap))) + global delete_count + + message("Identify intersections...\n") + + lap = time.time() + node_count = len(nodes) + delete_count = 0 + + # Create set of common nodes for segment intersections + + for segment in segments: + if segment["used"] > 0: + nodes.add(segment["coordinates"][0]) + nodes.add(segment["coordinates"][-1]) + + for feature in features: + if feature["type"] == "LineString": + nodes.add(feature["coordinates"][0]) + nodes.add(feature["coordinates"][-1]) + + if not no_node: + + # Create bbox for segments and line features + create temporary list of LineString features + + for segment in segments: + [segment["min_bbox"], segment["max_bbox"]] = get_bbox( + segment["coordinates"], 0 + ) + + # Loop streams to identify intersections with segments + + for feature in features: + if feature["type"] == "LineString" and feature["object"] == "ElvBekk": + [feature["min_bbox"], feature["max_bbox"]] = get_bbox( + feature["coordinates"], 0 + ) + + for segment in segments: + if (segment["used"] > 0 or debug) and ( + feature["min_bbox"][0] <= segment["max_bbox"][0] + and feature["max_bbox"][0] >= segment["min_bbox"][0] + and feature["min_bbox"][1] <= segment["max_bbox"][1] + and feature["max_bbox"][1] >= segment["min_bbox"][1] + ): + + intersections = set(feature["coordinates"]).intersection( + set(segment["coordinates"]) + ) + if len(intersections) == 0: + continue + + for node in intersections: + index1 = feature["coordinates"].index(node) + index2 = segment["coordinates"].index(node) + + # First check if stream node may be removed og slightly relocated + + if ( + index1 not in [0, len(feature["coordinates"]) - 1] + and node not in nodes + ): + if ( + feature["coordinates"][index1 - 1] + not in intersections + and feature["coordinates"][index1 + 1] + not in intersections + ): + feature["coordinates"].pop(index1) + else: + lon, lat = node + offset = 10 ** ( + -coordinate_decimals + 1 + ) # Last lat/lon decimal digit + feature["coordinates"][index1] = ( + lon + 4 * offset, + lat + 2 * offset, + ) + # Note: New node used in next test here + + # Then check if segment node may also be removed + + if index2 not in [0, len(segment["coordinates"]) - 1]: + if ( + segment["coordinates"][index2 - 1] + not in intersections + and segment["coordinates"][index2 + 1] + not in intersections + ): + segment["coordinates"].pop(index2) + + # Else create new common node with "water" segments (or reuse existing common node) + + elif segment["object"] in [ + "Kystkontur", + "Innsjøkant", + "InnsjøkantRegulert", + "ElvBekkKant", + ]: + nodes.add(node) + + # Loop auxiliary lines and simplify geometry + + for segment in segments: + if segment["used"] > 0 and segment["object"] == "FiktivDelelinje": + delete_count += len(segment["coordinates"]) - 2 + segment["coordinates"] = [ + segment["coordinates"][0], + segment["coordinates"][-1], + ] + + message( + "\t%i common nodes, %i nodes removed from streams and auxiliary lines\n" + % (len(nodes) - node_count, delete_count) + ) + message("\tRun time %s\n" % (timeformat(time.time() - lap))) # OLD: @@ -1258,71 +1486,73 @@ def match_nodes(): # and: https://wms.geonorge.no/skwms1/wps.elevation2?request=DescribeProcess&service=WPS&version=1.0.0&identifier=elevation # Note: Slow api (approx 1.4 calls/second), to be replaced by Kartverket end of 2020 -def get_elevation_old(node): - global elevations, ele_count, retry_count +def get_elevation_old(node): - ns_ows="http://www.opengis.net/ows/1.1" - ns_wps="http://www.opengis.net/wps/1.0.0" + global elevations, ele_count, retry_count - ns = { - 'ns0': ns_wps, - 'ns2': ns_ows - } + ns_ows = "http://www.opengis.net/ows/1.1" + ns_wps = "http://www.opengis.net/wps/1.0.0" - url = "https://wms.geonorge.no/skwms1/wps.elevation2?request=Execute&service=WPS&version=1.0.0&identifier=elevation" + \ - "&datainputs=lat=%f;lon=%f;epsg=4326" % (node[1], node[0]) + ns = {"ns0": ns_wps, "ns2": ns_ows} - if node in elevations: - return elevations[ node ] + url = ( + "https://wms.geonorge.no/skwms1/wps.elevation2?request=Execute&service=WPS&version=1.0.0&identifier=elevation" + + "&datainputs=lat=%f;lon=%f;epsg=4326" % (node[1], node[0]) + ) - request = urllib.request.Request(url, headers=header) + if node in elevations: + return elevations[node] - max_retry = 5 - for retry in range(1, max_retry + 1): - try: - file = urllib.request.urlopen(request) - if retry > 1: - message ("\n") - break + request = urllib.request.Request(url, headers=header) - except urllib.error.HTTPError as e: - if retry < max_retry: - message ("\tRetry #%i " % retry) - time.sleep(2 ** (retry-1)) # First retry after 1 second - retry_count += 1 - else: - message("\t*** ELEVATION API UNAVAILABLE\n") - tree = ET.parse(e) - root = tree.getroot() - exception_text = root[0][0].text + max_retry = 5 + for retry in range(1, max_retry + 1): + try: + file = urllib.request.urlopen(request) + if retry > 1: + message("\n") + break - if "exceptionCode" in root[0].attrib: - exception_code = root[0].attrib['exceptionCode'] - else: - exception_code = "" + except urllib.error.HTTPError as e: + if retry < max_retry: + message("\tRetry #%i " % retry) + time.sleep(2 ** (retry - 1)) # First retry after 1 second + retry_count += 1 + else: + message("\t*** ELEVATION API UNAVAILABLE\n") + tree = ET.parse(e) + root = tree.getroot() + exception_text = root[0][0].text - message ("\tHTTP Error %s: %s " % (e.code, e.reason)) - message ("\tException: [%s] %s\n" % (exception_code, exception_text)) - return None + if "exceptionCode" in root[0].attrib: + exception_code = root[0].attrib["exceptionCode"] + else: + exception_code = "" - tree = ET.parse(file) - file.close() - root = tree.getroot() + message("\tHTTP Error %s: %s " % (e.code, e.reason)) + message("\tException: [%s] %s\n" % (exception_code, exception_text)) + return None - ele_tag = root.find('ns0:ProcessOutputs/ns0:Output[ns2:Identifier="elevation"]/ns0:Data/ns0:LiteralData', ns) - ele = float(ele_tag.text) + tree = ET.parse(file) + file.close() + root = tree.getroot() - if math.isnan(ele): - ele = None - if ele is None: - message (" *** NO ELEVATION: %s \n" % str(node)) + ele_tag = root.find( + 'ns0:ProcessOutputs/ns0:Output[ns2:Identifier="elevation"]/ns0:Data/ns0:LiteralData', + ns, + ) + ele = float(ele_tag.text) - elevations[ node ] = ele # Store for possible identical request later - ele_count += 1 + if math.isnan(ele): + ele = None + if ele is None: + message(" *** NO ELEVATION: %s \n" % str(node)) - return ele + elevations[node] = ele # Store for possible identical request later + ele_count += 1 + return ele # NEW: @@ -1330,431 +1560,518 @@ def get_elevation_old(node): # Note: Still testing this api # Note: Slow api, approx 4 calls / second + def get_elevation(node): - global elevations, ele_count, retry_count + global elevations, ele_count, retry_count - url = "https://ws.geonorge.no/hoydedata/v1/punkt?nord=%f&ost=%f&geojson=false&koordsys=4258" % (node[1], node[0]) -# url = "https://wstest.geonorge.no/hoydedata/v1/punkt?nord=%f&ost=%f&geojson=false&koordsys=4258" % (node[1], node[0]) -# url = "https://wstest.geonorge.no/hoydedata/v1/punkt?nord=60.1&koordsys=4258&geojson=false&ost=11.1" -# message ("%s\n" % url) + url = ( + "https://ws.geonorge.no/hoydedata/v1/punkt?nord=%f&ost=%f&geojson=false&koordsys=4258" + % (node[1], node[0]) + ) + # url = "https://wstest.geonorge.no/hoydedata/v1/punkt?nord=%f&ost=%f&geojson=false&koordsys=4258" % (node[1], node[0]) + # url = "https://wstest.geonorge.no/hoydedata/v1/punkt?nord=60.1&koordsys=4258&geojson=false&ost=11.1" + # message ("%s\n" % url) - if node in elevations: - return elevations[ node ] + if node in elevations: + return elevations[node] - request = urllib.request.Request(url, headers=header) - file = urllib.request.urlopen(request) + request = urllib.request.Request(url, headers=header) + file = urllib.request.urlopen(request) - result = json.load(file) - file.close() -# message ("Result: %s\n" % str(result)) -# ele = result['features'][0]['geometry']['coordinates'][2] - ele = result['punkter'][0]['z'] - elevations[ node ] = ele # Store for possible identical request later - ele_count += 1 + result = json.load(file) + file.close() + # message ("Result: %s\n" % str(result)) + # ele = result['features'][0]['geometry']['coordinates'][2] + ele = result["punkter"][0]["z"] + elevations[node] = ele # Store for possible identical request later + ele_count += 1 - if ele is None: - message (" *** NO ELEVATION: %s \n" % str(result)) -# ele = 0.0 + if ele is None: + message(" *** NO ELEVATION: %s \n" % str(result)) + # ele = 0.0 - return ele + return ele # Turn streams which have uphill direction -def fix_stream_direction(): - - global elevations, ele_count, retry_count - - elevations = {} # Already fetched elevations from api - ele_count = 0 # Numbr of api calls during program execution - retry_count = 0 # Number of retry to api - max_error = 1.0 # Meters of elevation difference - - lap = time.time() - - message ("Load elevation data from Kartverket and reverse streams...\n") - - stream_count = 0 - for feature in features: - if feature['object'] == "ElvBekk" and feature['type'] == "LineString": - stream_count += 1 - - message ("\t%i streams\n" % stream_count) - api_count = stream_count - reverse_count = 0 - - # Loop all streams and check elevation difference between first and last nodes - - for feature in features: - if feature['object'] == "ElvBekk" and feature['type'] == "LineString": - api_count -= 1 - - message ("\r\t%i " % api_count) - - ele_start = get_elevation(feature['coordinates'][0]) - if ele_start is None: - continue - - ele_end = get_elevation(feature['coordinates'][-1]) - if ele_end is None: - continue +def fix_stream_direction(): - # Reverse direction of stream if within error margin + global elevations, ele_count, retry_count - if ele_end - ele_start >= max_error: - feature['coordinates'].reverse() - reverse_count += 1 - feature['extras']['reversert'] = "%.2f" % (ele_end - ele_start) - else: - feature['extras']['bekk'] = "%.2f" % (ele_end - ele_start) + elevations = {} # Already fetched elevations from api + ele_count = 0 # Numbr of api calls during program execution + retry_count = 0 # Number of retry to api + max_error = 1.0 # Meters of elevation difference - if abs(ele_end - ele_start) < 2 * max_error: - feature['tags']['fixme'] = "Please check direction (%.1fm elevation)" % (ele_end - ele_start) + lap = time.time() - if retry_count > 0: - message ("\t%i retry to api\n" % retry_count) + message("Load elevation data from Kartverket and reverse streams...\n") - duration = time.time() - lap - message ("\r\t%i streams, %i reversed, %i api calls\n" \ - % (stream_count, reverse_count, ele_count)) - message ("\tRun time %s, %.2f streams per second\n" % (timeformat(duration), stream_count / duration)) + stream_count = 0 + for feature in features: + if feature["object"] == "ElvBekk" and feature["type"] == "LineString": + stream_count += 1 + message("\t%i streams\n" % stream_count) -# Load place names from SSR within given bbox + api_count = stream_count + reverse_count = 0 -def get_ssr_name (feature, name_categories): + # Loop all streams and check elevation difference between first and last nodes - global ssr_places, name_count + for feature in features: + if feature["object"] == "ElvBekk" and feature["type"] == "LineString": + api_count -= 1 - if feature['type'] == "Point": - bbox = get_bbox(feature['coordinates'], 500) # 500 meters perimeter to each side - else: - bbox = get_bbox(feature['coordinates'], 0) + message("\r\t%i " % api_count) - if type(feature['coordinates'][0]) is tuple: - polygon = feature['coordinates'] - else: - polygon = feature['coordinates'][0] + ele_start = get_elevation(feature["coordinates"][0]) + if ele_start is None: + continue - found_names = [] - names = [] + ele_end = get_elevation(feature["coordinates"][-1]) + if ele_end is None: + continue - # Find name in stored file + # Reverse direction of stream if within error margin - for place in ssr_places: - if bbox[0][0] <= place['coordinate'][0] <= bbox[1][0] and bbox[0][1] <= place['coordinate'][1] <= bbox[1][1] and \ - place['tags']['ssr:type'] in name_categories and place['tags']['name'] not in names and \ - (feature['type'] == "Point" or inside_polygon(place['coordinate'], polygon)): - found_names.append(place) - names.extend(place['tags']['name'].replace(" - ", ";").split(";")) # Split multiple languages + variants + if ele_end - ele_start >= max_error: + feature["coordinates"].reverse() + reverse_count += 1 + feature["extras"]["reversert"] = "%.2f" % (ele_end - ele_start) + else: + feature["extras"]["bekk"] = "%.2f" % (ele_end - ele_start) - # Sort and select name + if abs(ele_end - ele_start) < 2 * max_error: + feature["tags"][ + "fixme" + ] = "Please check direction (%.1fm elevation)" % (ele_end - ele_start) - if found_names: + if retry_count > 0: + message("\t%i retry to api\n" % retry_count) - found_names.sort(key=lambda name: name_categories.index(name['tags']['ssr:type'])) # place_name_rank(name, name_categories)) + duration = time.time() - lap + message( + "\r\t%i streams, %i reversed, %i api calls\n" + % (stream_count, reverse_count, ele_count) + ) + message( + "\tRun time %s, %.2f streams per second\n" + % (timeformat(duration), stream_count / duration) + ) - alt_names = [] - sea = ("Sjø" in found_names[0]['tags']['ssr:type']) - for place in found_names: - if not sea or "Sjø" in place['tags']['ssr:type']: - alt_names.append("%s [%s %s]" % (place['tags']['name'], place['tags']['ssr:type'], place['tags']['ssr:stedsnr'])) - # Name already suggested by NVE data, so get ssr:stedsnr and any alternative names - if "name" in feature['tags'] and feature['tags']['name'] in names: - name = feature['tags']['name'] - for place in found_names: - if name in place['tags']['name'].replace(" - ", ";").split(";"): - feature['tags'].update(place['tags']) - feature['tags']['name'] = name # Add back NVE name - feature['extras']['ssr:type'] = feature['tags'].pop("ssr:type", None) - if len(alt_names) > 1: - if name != place['tags']['name']: - alt_names.insert(0, "%s [NVE]" % name) - feature['tags']['fixme'] = "Alternative names: " + ", ".join(alt_names) - name_count += 1 - break +# Load place names from SSR within given bbox - # Only one name found, or only one name of preferred type - elif len(alt_names) == 1 or \ - "øyISjø" in found_names[0]['tags']['ssr:type'] and "øyISjø" not in found_names[1]['tags']['ssr:type'] or \ - "øy" in found_names[0]['tags']['ssr:type'] and "øy" not in found_names[1]['tags']['ssr:type'] or \ - "holme" in found_names[0]['tags']['ssr:type'] and "holme" not in found_names[1]['tags']['ssr:type']: - if "name" in feature['tags'] and (len(alt_names) > 1 or \ - feature['tags']['name'] not in found_names[0]['tags']['name'].replace(" - ", ";").split(";")): - alt_names.insert(0, "%s [NVE]" % feature['tags']['name']) +def get_ssr_name(feature, name_categories): - feature['tags'].update(found_names[0]['tags']) - feature['extras']['ssr:type'] = feature['tags'].pop("ssr:type", None) - if len(alt_names) > 1: - feature['tags']['fixme'] = "Alternative names: " + ", ".join(alt_names) - name_count += 1 + global ssr_places, name_count - else: - feature['tags']['fixme'] = "Consider names: " + ", ".join(alt_names) + if feature["type"] == "Point": + bbox = get_bbox( + feature["coordinates"], 500 + ) # 500 meters perimeter to each side + else: + bbox = get_bbox(feature["coordinates"], 0) - return found_names[0]['coordinate'] + if type(feature["coordinates"][0]) is tuple: + polygon = feature["coordinates"] + else: + polygon = feature["coordinates"][0] + + found_names = [] + names = [] + + # Find name in stored file + + for place in ssr_places: + if ( + bbox[0][0] <= place["coordinate"][0] <= bbox[1][0] + and bbox[0][1] <= place["coordinate"][1] <= bbox[1][1] + and place["tags"]["ssr:type"] in name_categories + and place["tags"]["name"] not in names + and ( + feature["type"] == "Point" + or inside_polygon(place["coordinate"], polygon) + ) + ): + found_names.append(place) + names.extend( + place["tags"]["name"].replace(" - ", ";").split(";") + ) # Split multiple languages + variants + + # Sort and select name + + if found_names: + + found_names.sort( + key=lambda name: name_categories.index(name["tags"]["ssr:type"]) + ) # place_name_rank(name, name_categories)) + + alt_names = [] + sea = "Sjø" in found_names[0]["tags"]["ssr:type"] + for place in found_names: + if not sea or "Sjø" in place["tags"]["ssr:type"]: + alt_names.append( + "%s [%s %s]" + % ( + place["tags"]["name"], + place["tags"]["ssr:type"], + place["tags"]["ssr:stedsnr"], + ) + ) + + # Name already suggested by NVE data, so get ssr:stedsnr and any alternative names + if "name" in feature["tags"] and feature["tags"]["name"] in names: + name = feature["tags"]["name"] + for place in found_names: + if name in place["tags"]["name"].replace(" - ", ";").split(";"): + feature["tags"].update(place["tags"]) + feature["tags"]["name"] = name # Add back NVE name + feature["extras"]["ssr:type"] = feature["tags"].pop( + "ssr:type", None + ) + if len(alt_names) > 1: + if name != place["tags"]["name"]: + alt_names.insert(0, "%s [NVE]" % name) + feature["tags"]["fixme"] = "Alternative names: " + ", ".join( + alt_names + ) + name_count += 1 + break + + # Only one name found, or only one name of preferred type + elif ( + len(alt_names) == 1 + or "øyISjø" in found_names[0]["tags"]["ssr:type"] + and "øyISjø" not in found_names[1]["tags"]["ssr:type"] + or "øy" in found_names[0]["tags"]["ssr:type"] + and "øy" not in found_names[1]["tags"]["ssr:type"] + or "holme" in found_names[0]["tags"]["ssr:type"] + and "holme" not in found_names[1]["tags"]["ssr:type"] + ): + + if "name" in feature["tags"] and ( + len(alt_names) > 1 + or feature["tags"]["name"] + not in found_names[0]["tags"]["name"].replace(" - ", ";").split(";") + ): + alt_names.insert(0, "%s [NVE]" % feature["tags"]["name"]) + + feature["tags"].update(found_names[0]["tags"]) + feature["extras"]["ssr:type"] = feature["tags"].pop("ssr:type", None) + if len(alt_names) > 1: + feature["tags"]["fixme"] = "Alternative names: " + ", ".join(alt_names) + name_count += 1 + + else: + feature["tags"]["fixme"] = "Consider names: " + ", ".join(alt_names) + + return found_names[0]["coordinate"] - else: - return None + else: + return None # Get place names for a category + def get_category_place_names(n50_categories, ssr_categories): - for feature in features: - if feature['object'] in n50_categories: - get_ssr_name(feature, ssr_categories) + for feature in features: + if feature["object"] in n50_categories: + get_ssr_name(feature, ssr_categories) # Get place names for islands, glaciers etc. # Place name categories: https://github.com/osmno/geocode2osm/blob/master/navnetyper.json -def get_place_names(): - - global ssr_places, name_count - global elevations, ele_count, retry_count - - message ("Load place names from SSR...\n") - - elevations = {} - ele_count = 0 - retry_count = 0 - - lap = time.time() - name_count = 0 - ssr_places = [] - - # Load all SSR place names in municipality - - url = "https://obtitus.github.io/ssr2_to_osm_data/data/%s/%s.osm" % (municipality_id, municipality_id) - request = urllib.request.Request(url, headers=header) - file = urllib.request.urlopen(request) - - tree = ET.parse(file) - file.close() - root = tree.getroot() - - for node in root.iter('node'): - entry = { - 'coordinate': (float(node.get('lon')), float(node.get('lat'))), - 'tags': {} - } - for tag in node.iter('tag'): - if "name" in tag.get('k') or tag.get('k') in ["ssr:stedsnr", "ssr:type"]: - if tag.get('k') == "fixme": - if "multiple name" in tag.get('v'): - entry['tags']['fixme'] = "Multiple name tags, please choose one and add the other to alt_name" - else: - entry['tags'][ tag.get('k') ] = tag.get('v') - - ssr_places.append(entry) - - message ("\t%s place names in SSR file\n" % len(ssr_places)) - - # Get island names - - name_category = ["øyISjø", "øygruppeISjø", "holmeISjø", "skjærISjø", "øy", "øygruppe", "holme", "skjær"] - for elements in [segments, features]: - for element in elements: - if "place" in element['tags'] and element['tags']['place'] in ["island", "islet"]: - get_ssr_name(element, name_category) - - # Get lake names + missing elevations for lakes - - if lake_ele and data_category == "Arealdekke": - message ("\tLoading lake elevations...\n") - - name_category = ["innsjø", "delAvInnsjø", "vann", "gruppeAvVann", "kanal", "gruppeAvTjern", "tjern", "lon", "pytt"] - lake_ele_count = 0 - - for feature in features: - if feature['object'] in ["Innsjø", "InnsjøRegulert"]: - - lake_node = get_ssr_name(feature, name_category) - area = abs(multipolygon_area(feature['coordinates'])) - feature['extras']['areal'] = str(int(area)) - - # Get lake's elevation if missing, and if lake is larger than threshold - - if lake_ele and "ele" not in feature['tags'] and (lake_node or area >= lake_ele_size): - - # Check that name coordinate is not on lake's island - if lake_node: - if not inside_multipolygon(lake_node, feature['coordinates']): - lake_node = None - else: - feature['extras']['elevation'] = "Based on lake name position" - - # If name coordinate cannot be used, try centroid - if lake_node is None: - lake_node = polygon_centroid(feature['coordinates'][0]) - feature['extras']['elevation'] = "Based on centroid" - if not inside_multipolygon(lake_node, feature['coordinates']): - lake_node = None - - # If all fail, just use coordinate of first node on lake perimeter - if lake_node is None: - lake_node = feature['coordinates'][0][0] - feature['extras']['elevation'] = "Based on first node" - - if lake_node: - ele = get_elevation(lake_node) - if ele: - feature['tags']['ele'] = str(int(round(ele))) - create_point(lake_node, "", "elevation %.1f" % ele) - lake_ele_count += 1 - message ("\r\t%i " % lake_ele_count) - - # Create lake centroid nodes for debugging - for feature in features: - if feature['object'] in ["Innsjø", "InnsjøRegulert"]: - centroid = polygon_centroid(feature['coordinates'][0]) - create_point(centroid, feature['gml_id'], "centroid") - - # Get glacier names - get_category_place_names(["SnøIsbre"], ["isbre", "fonn", "iskuppel"]) - - # Get wetland names - get_category_place_names(["Myr"], ["myr", "våtmarksområde"]) - - # Get cemetery names - get_category_place_names(["Gravplass"], ["gravplass"]) - - # Get piste names - get_category_place_names(["Alpinbakke"], ["alpinanlegg", "skiheis"]) - # Get dam names - get_category_place_names(["Steintipp"], ["dam"]) - - # Get waterfall (point) - get_category_place_names(["Foss"], ["foss", "stryk"]) +def get_place_names(): - message ("\r\t%i place names found\n" % name_count) - if lake_ele: - message ("\t%i lake elevations found\n" % lake_ele_count) - message ("\tRun time %s\n" % (timeformat(time.time() - lap))) + global ssr_places, name_count + global elevations, ele_count, retry_count + + message("Load place names from SSR...\n") + + elevations = {} + ele_count = 0 + retry_count = 0 + + lap = time.time() + name_count = 0 + ssr_places = [] + + # Load all SSR place names in municipality + + url = "https://obtitus.github.io/ssr2_to_osm_data/data/%s/%s.osm" % ( + municipality_id, + municipality_id, + ) + request = urllib.request.Request(url, headers=header) + file = urllib.request.urlopen(request) + + tree = ET.parse(file) + file.close() + root = tree.getroot() + + for node in root.iter("node"): + entry = { + "coordinate": (float(node.get("lon")), float(node.get("lat"))), + "tags": {}, + } + for tag in node.iter("tag"): + if "name" in tag.get("k") or tag.get("k") in ["ssr:stedsnr", "ssr:type"]: + if tag.get("k") == "fixme": + if "multiple name" in tag.get("v"): + entry["tags"][ + "fixme" + ] = "Multiple name tags, please choose one and add the other to alt_name" + else: + entry["tags"][tag.get("k")] = tag.get("v") + + ssr_places.append(entry) + + message("\t%s place names in SSR file\n" % len(ssr_places)) + + # Get island names + + name_category = [ + "øyISjø", + "øygruppeISjø", + "holmeISjø", + "skjærISjø", + "øy", + "øygruppe", + "holme", + "skjær", + ] + for elements in [segments, features]: + for element in elements: + if "place" in element["tags"] and element["tags"]["place"] in [ + "island", + "islet", + ]: + get_ssr_name(element, name_category) + + # Get lake names + missing elevations for lakes + + if lake_ele and data_category == "Arealdekke": + message("\tLoading lake elevations...\n") + + name_category = [ + "innsjø", + "delAvInnsjø", + "vann", + "gruppeAvVann", + "kanal", + "gruppeAvTjern", + "tjern", + "lon", + "pytt", + ] + lake_ele_count = 0 + + for feature in features: + if feature["object"] in ["Innsjø", "InnsjøRegulert"]: + + lake_node = get_ssr_name(feature, name_category) + area = abs(multipolygon_area(feature["coordinates"])) + feature["extras"]["areal"] = str(int(area)) + + # Get lake's elevation if missing, and if lake is larger than threshold + + if ( + lake_ele + and "ele" not in feature["tags"] + and (lake_node or area >= lake_ele_size) + ): + + # Check that name coordinate is not on lake's island + if lake_node: + if not inside_multipolygon(lake_node, feature["coordinates"]): + lake_node = None + else: + feature["extras"]["elevation"] = "Based on lake name position" + + # If name coordinate cannot be used, try centroid + if lake_node is None: + lake_node = polygon_centroid(feature["coordinates"][0]) + feature["extras"]["elevation"] = "Based on centroid" + if not inside_multipolygon(lake_node, feature["coordinates"]): + lake_node = None + + # If all fail, just use coordinate of first node on lake perimeter + if lake_node is None: + lake_node = feature["coordinates"][0][0] + feature["extras"]["elevation"] = "Based on first node" + + if lake_node: + ele = get_elevation(lake_node) + if ele: + feature["tags"]["ele"] = str(int(round(ele))) + create_point(lake_node, "", "elevation %.1f" % ele) + lake_ele_count += 1 + message("\r\t%i " % lake_ele_count) + + # Create lake centroid nodes for debugging + for feature in features: + if feature["object"] in ["Innsjø", "InnsjøRegulert"]: + centroid = polygon_centroid(feature["coordinates"][0]) + create_point(centroid, feature["gml_id"], "centroid") + + # Get glacier names + get_category_place_names(["SnøIsbre"], ["isbre", "fonn", "iskuppel"]) + + # Get wetland names + get_category_place_names(["Myr"], ["myr", "våtmarksområde"]) + + # Get cemetery names + get_category_place_names(["Gravplass"], ["gravplass"]) + + # Get piste names + get_category_place_names(["Alpinbakke"], ["alpinanlegg", "skiheis"]) + + # Get dam names + get_category_place_names(["Steintipp"], ["dam"]) + + # Get waterfall (point) + get_category_place_names(["Foss"], ["foss", "stryk"]) + + message("\r\t%i place names found\n" % name_count) + if lake_ele: + message("\t%i lake elevations found\n" % lake_ele_count) + message("\tRun time %s\n" % (timeformat(time.time() - lap))) # Get name from NVE Innsjødatabasen # API reference: https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer -def get_nve_lakes(): - - message ("Load lake data from NVE...\n") - - n50_lake_count = 0 - nve_lake_count = 0 - more_lakes = True - lakes = {} - - # Paging results (default 1000 lakes) - - while more_lakes: - -# Alternative url: -# url = "https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer/find?" + \ -# "searchText=%s&contains=true&searchFields=kommNr&layers=5&returnGeometry=false&returnUnformattedValues=true&f=pjson&resultOffset=%i&resultRecordCount=1000&orderByFields=areal_km2%20DESC" \ -# % (municipality_id, nve_lake_count) - -# url = "https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer/5/query?" + \ -# "where=kommNr%%3D%%27%s%%27&outFields=vatnLnr%%2Cnavn%%2Choyde%%2Careal_km2%%2CmagasinNr&returnGeometry=false&resultOffset=%i&resultRecordCount=1000&f=json" \ -# % (municipality_id, nve_lake_count) - - url = "https://nve.geodataonline.no/arcgis/rest/services/Innsjodatabase2/MapServer/5/query?" + \ - "where=kommNr%%3D%%27%s%%27&outFields=vatnLnr%%2Cnavn%%2Choyde%%2Careal_km2%%2CmagasinNr&returnGeometry=false&resultOffset=%i&resultRecordCount=1000&f=json" \ - % (municipality_id, nve_lake_count) - - - request = urllib.request.Request(url, headers=header) - file = urllib.request.urlopen(request) - lake_data = json.load(file) - file.close() - - for lake_result in lake_data['features']: - lake = lake_result['attributes'] - entry = { - 'name': lake['navn'], - 'ele': lake['hoyde'], - 'area': lake['areal_km2'], - 'mag_id': lake['magasinNr'] - } - lakes[ str(lake['vatnLnr']) ] = entry - nve_lake_count += len(lake_data['features']) - - if 'exceededTransferLimit' not in lake_data: - more_lakes = False - - # Update lake info - - for feature in features: - if "ref:nve:vann" in feature['tags']: - ref = feature['tags']['ref:nve:vann'] - if ref in lakes: - if lakes[ref]['name']: - feature['tags']['name'] = lakes[ref]['name'] - if lakes[ref]['ele'] and "ele" not in feature['tags']: - feature['tags']['ele'] = str(lakes[ref]['ele']) - if lakes[ref]['area'] > 1 and "water" not in feature['tags']: - feature['tags']['water'] = "lake" - if lakes[ref]['mag_id']: - feature['tags']['ref:nve:magasin'] = str(lakes[ref]['mag_id']) - feature['extras']['nve_areal'] = str(int(lakes[ref]['area'] * 1000000)) # Square meters - - if feature['object'] in ["Innsjø", "InnsjøRegulert"]: - n50_lake_count += 1 +def get_nve_lakes(): - message ("\t%i N50 lakes matched against %i NVE lakes\n" % (n50_lake_count, nve_lake_count)) + message("Load lake data from NVE...\n") + + n50_lake_count = 0 + nve_lake_count = 0 + more_lakes = True + lakes = {} + + # Paging results (default 1000 lakes) + + while more_lakes: + + # Alternative url: + # url = "https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer/find?" + \ + # "searchText=%s&contains=true&searchFields=kommNr&layers=5&returnGeometry=false&returnUnformattedValues=true&f=pjson&resultOffset=%i&resultRecordCount=1000&orderByFields=areal_km2%20DESC" \ + # % (municipality_id, nve_lake_count) + + # url = "https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer/5/query?" + \ + # "where=kommNr%%3D%%27%s%%27&outFields=vatnLnr%%2Cnavn%%2Choyde%%2Careal_km2%%2CmagasinNr&returnGeometry=false&resultOffset=%i&resultRecordCount=1000&f=json" \ + # % (municipality_id, nve_lake_count) + + url = ( + "https://nve.geodataonline.no/arcgis/rest/services/Innsjodatabase2/MapServer/5/query?" + + "where=kommNr%%3D%%27%s%%27&outFields=vatnLnr%%2Cnavn%%2Choyde%%2Careal_km2%%2CmagasinNr&returnGeometry=false&resultOffset=%i&resultRecordCount=1000&f=json" + % (municipality_id, nve_lake_count) + ) + + request = urllib.request.Request(url, headers=header) + file = urllib.request.urlopen(request) + lake_data = json.load(file) + file.close() + + for lake_result in lake_data["features"]: + lake = lake_result["attributes"] + entry = { + "name": lake["navn"], + "ele": lake["hoyde"], + "area": lake["areal_km2"], + "mag_id": lake["magasinNr"], + } + lakes[str(lake["vatnLnr"])] = entry + + nve_lake_count += len(lake_data["features"]) + + if "exceededTransferLimit" not in lake_data: + more_lakes = False + + # Update lake info + + for feature in features: + if "ref:nve:vann" in feature["tags"]: + ref = feature["tags"]["ref:nve:vann"] + if ref in lakes: + if lakes[ref]["name"]: + feature["tags"]["name"] = lakes[ref]["name"] + if lakes[ref]["ele"] and "ele" not in feature["tags"]: + feature["tags"]["ele"] = str(lakes[ref]["ele"]) + if lakes[ref]["area"] > 1 and "water" not in feature["tags"]: + feature["tags"]["water"] = "lake" + if lakes[ref]["mag_id"]: + feature["tags"]["ref:nve:magasin"] = str(lakes[ref]["mag_id"]) + feature["extras"]["nve_areal"] = str( + int(lakes[ref]["area"] * 1000000) + ) # Square meters + + if feature["object"] in ["Innsjø", "InnsjøRegulert"]: + n50_lake_count += 1 + + message( + "\t%i N50 lakes matched against %i NVE lakes\n" + % (n50_lake_count, nve_lake_count) + ) # Save geojson file for reviewing raw input data from GML file -def save_geojson(filename): - message ("Save to '%s' file...\n" % filename) +def save_geojson(filename): - json_features = { - 'type': 'FeatureCollection', - 'features': [] - } + message("Save to '%s' file...\n" % filename) - for feature_list in [features, segments]: - for feature in feature_list: - entry = { - 'type': 'Feature', - 'geometry': { - 'type': feature['type'], - 'coordinates': feature['coordinates'] - }, - 'properties': dict(feature['extras'].items() + feature['tags'].items() + { 'gml_id': gml_id }.items()) - } + json_features = {"type": "FeatureCollection", "features": []} - json_features['features'].append(entry) + for feature_list in [features, segments]: + for feature in feature_list: + entry = { + "type": "Feature", + "geometry": { + "type": feature["type"], + "coordinates": feature["coordinates"], + }, + "properties": dict( + feature["extras"].items() + + feature["tags"].items() + + {"gml_id": gml_id}.items() + ), + } - file = open(filename, "w") - json.dump(json_features, file, indent=2) - file.close() + json_features["features"].append(entry) - message ("\t%i features saved\n" % len(features)) + file = open(filename, "w") + json.dump(json_features, file, indent=2) + file.close() + message("\t%i features saved\n" % len(features)) # Indent XML output + def indent_tree(elem, level=0): - i = "\n" + level*" " + i = "\n" + level * " " if len(elem): if not elem.text or not elem.text.strip(): elem.text = i + " " if not elem.tail or not elem.tail.strip(): elem.tail = i for elem in elem: - indent_tree(elem, level+1) + indent_tree(elem, level + 1) if not elem.tail or not elem.tail.strip(): elem.tail = i else: @@ -1762,230 +2079,276 @@ def indent_tree(elem, level=0): elem.tail = i - # Save osm file + def save_osm(filename): - message ("Save to '%s' file...\n" % filename) - - osm_node_ids = {} # Will contain osm_id of each common node - relation_count = 0 - way_count = 0 - node_count = 0 - - osm_root = ET.Element("osm", version="0.6", generator="n50osm v"+version, upload="false") - osm_id = -1000 - - # Common nodes - - for node in nodes: - osm_id -= 1 - osm_node = ET.Element("node", id=str(osm_id), action="modify", lat=str(node[1]), lon=str(node[0])) - osm_root.append(osm_node) - osm_node_ids[ node ] = osm_id - node_count += 1 - - # Ways used by relations - - for segment in segments: - if segment['used'] > 0 or debug: # or segment['object'] == "Kystkontur": - osm_id -= 1 - osm_feature = ET.Element("way", id=str(osm_id), action="modify") - osm_root.append(osm_feature) - segment['osm_id'] = osm_id - segment['etree'] = osm_feature - way_count += 1 - - for node in segment['coordinates']: - if node in nodes: - osm_nd = ET.Element("nd", ref=str(osm_node_ids[ node ])) - else: - osm_id -= 1 - osm_node = ET.Element("node", id=str(osm_id), action="modify", lat=str(node[1]), lon=str(node[0])) - osm_root.append(osm_node) - osm_nd = ET.Element("nd", ref=str(osm_id)) - node_count += 1 - osm_feature.append(osm_nd) - - for key, value in iter(segment['tags'].items()): - osm_tag = ET.Element("tag", k=key, v=value) - osm_feature.append(osm_tag) - - if debug: - for key, value in iter(segment['extras'].items()): - osm_tag = ET.Element("tag", k=key.upper(), v=value) - osm_feature.append(osm_tag) - - # The main objects - - for feature in features: - - if feature['object'] == "Havflate": - continue - - if feature['type'] == "Point": - osm_id -= 1 - osm_feature = ET.Element("node", id=str(osm_id), action="modify", lat=str(feature['coordinates'][1]), lon=str(feature['coordinates'][0])) - osm_root.append(osm_feature) - node_count += 1 - - elif feature['type'] in "LineString": - osm_id -= 1 - osm_feature = ET.Element("way", id=str(osm_id), action="modify") - osm_root.append(osm_feature) - way_count += 1 - - for node in feature['coordinates']: - if node in nodes: - osm_nd = ET.Element("nd", ref=str(osm_node_ids [ node ])) - else: - osm_id -= 1 - osm_node = ET.Element("node", id=str(osm_id), action="modify", lat=str(node[1]), lon=str(node[0])) - osm_root.append(osm_node) - osm_nd = ET.Element("nd", ref=str(osm_id)) - node_count += 1 - osm_feature.append(osm_nd) - - elif feature['type'] == "Polygon": - - # Output way if possible to avoid relation - if len(feature['members']) == 1 and len(feature['members'][0]) == 1 and \ - not ("natural" in feature['tags'] and "natural" in segments[ feature['members'][0][0] ]['tags']): - osm_feature = segments[ feature['members'][0][0] ]['etree'] - - else: - osm_id -= 1 - osm_feature = ET.Element("relation", id=str(osm_id), action="modify") - osm_root.append(osm_feature) - relation_count += 1 - role = "outer" - - for patch in feature['members']: - for member in patch: - osm_member = ET.Element("member", type="way", ref=str(segments[member]['osm_id']), role=role) - osm_feature.append(osm_member) - role = "inner" - - osm_tag = ET.Element("tag", k="type", v="multipolygon") - osm_feature.append(osm_tag) - - else: - message ("\t*** UNKNOWN GEOMETRY: %s\n" % feature['type']) - - for key, value in iter(feature['tags'].items()): - osm_tag = ET.Element("tag", k=key, v=value) - osm_feature.append(osm_tag) - - if debug: - for key, value in iter(feature['extras'].items()): - osm_tag = ET.Element("tag", k=key.upper(), v=value) - osm_feature.append(osm_tag) - - - osm_root.set("upload", "false") - indent_tree(osm_root) - osm_tree = ET.ElementTree(osm_root) - osm_tree.write(filename, encoding='utf-8', method='xml', xml_declaration=True) - - message ("\t%i relations, %i ways, %i nodes saved\n" % (relation_count, way_count, node_count)) + message("Save to '%s' file...\n" % filename) + + osm_node_ids = {} # Will contain osm_id of each common node + relation_count = 0 + way_count = 0 + node_count = 0 + + osm_root = ET.Element( + "osm", version="0.6", generator="n50osm v" + version, upload="false" + ) + osm_id = -1000 + + # Common nodes + + for node in nodes: + osm_id -= 1 + osm_node = ET.Element( + "node", id=str(osm_id), action="modify", lat=str(node[1]), lon=str(node[0]) + ) + osm_root.append(osm_node) + osm_node_ids[node] = osm_id + node_count += 1 + + # Ways used by relations + + for segment in segments: + if segment["used"] > 0 or debug: # or segment['object'] == "Kystkontur": + osm_id -= 1 + osm_feature = ET.Element("way", id=str(osm_id), action="modify") + osm_root.append(osm_feature) + segment["osm_id"] = osm_id + segment["etree"] = osm_feature + way_count += 1 + + for node in segment["coordinates"]: + if node in nodes: + osm_nd = ET.Element("nd", ref=str(osm_node_ids[node])) + else: + osm_id -= 1 + osm_node = ET.Element( + "node", + id=str(osm_id), + action="modify", + lat=str(node[1]), + lon=str(node[0]), + ) + osm_root.append(osm_node) + osm_nd = ET.Element("nd", ref=str(osm_id)) + node_count += 1 + osm_feature.append(osm_nd) + + for key, value in iter(segment["tags"].items()): + osm_tag = ET.Element("tag", k=key, v=value) + osm_feature.append(osm_tag) + + if debug: + for key, value in iter(segment["extras"].items()): + osm_tag = ET.Element("tag", k=key.upper(), v=value) + osm_feature.append(osm_tag) + + # The main objects + + for feature in features: + + if feature["object"] == "Havflate": + continue + + if feature["type"] == "Point": + osm_id -= 1 + osm_feature = ET.Element( + "node", + id=str(osm_id), + action="modify", + lat=str(feature["coordinates"][1]), + lon=str(feature["coordinates"][0]), + ) + osm_root.append(osm_feature) + node_count += 1 + + elif feature["type"] in "LineString": + osm_id -= 1 + osm_feature = ET.Element("way", id=str(osm_id), action="modify") + osm_root.append(osm_feature) + way_count += 1 + + for node in feature["coordinates"]: + if node in nodes: + osm_nd = ET.Element("nd", ref=str(osm_node_ids[node])) + else: + osm_id -= 1 + osm_node = ET.Element( + "node", + id=str(osm_id), + action="modify", + lat=str(node[1]), + lon=str(node[0]), + ) + osm_root.append(osm_node) + osm_nd = ET.Element("nd", ref=str(osm_id)) + node_count += 1 + osm_feature.append(osm_nd) + + elif feature["type"] == "Polygon": + + # Output way if possible to avoid relation + if ( + len(feature["members"]) == 1 + and len(feature["members"][0]) == 1 + and not ( + "natural" in feature["tags"] + and "natural" in segments[feature["members"][0][0]]["tags"] + ) + ): + osm_feature = segments[feature["members"][0][0]]["etree"] + + else: + osm_id -= 1 + osm_feature = ET.Element("relation", id=str(osm_id), action="modify") + osm_root.append(osm_feature) + relation_count += 1 + role = "outer" + + for patch in feature["members"]: + for member in patch: + osm_member = ET.Element( + "member", + type="way", + ref=str(segments[member]["osm_id"]), + role=role, + ) + osm_feature.append(osm_member) + role = "inner" + + osm_tag = ET.Element("tag", k="type", v="multipolygon") + osm_feature.append(osm_tag) + + else: + message("\t*** UNKNOWN GEOMETRY: %s\n" % feature["type"]) + + for key, value in iter(feature["tags"].items()): + osm_tag = ET.Element("tag", k=key, v=value) + osm_feature.append(osm_tag) + + if debug: + for key, value in iter(feature["extras"].items()): + osm_tag = ET.Element("tag", k=key.upper(), v=value) + osm_feature.append(osm_tag) + + osm_root.set("upload", "false") + indent_tree(osm_root) + osm_tree = ET.ElementTree(osm_root) + osm_tree.write(filename, encoding="utf-8", method="xml", xml_declaration=True) + + message( + "\t%i relations, %i ways, %i nodes saved\n" + % (relation_count, way_count, node_count) + ) # Main program -if __name__ == '__main__': - - start_time = time.time() - message ("\n-- n50osm v%s --\n" % version) - - features = [] # All geometry and tags - segments = [] # Line segments which are shared by one or more polygons - nodes = set() # Common nodes at intersections, including start/end nodes of segments [lon,lat] - building_tags = {} # Conversion table from building type to osm tag - object_count = {} # Count loaded object types - - debug = False # Include debug tags and unused segments - n50_tags = False # Include property tags from N50 in output - json_output = False # Output complete and unprocessed geometry in geojson format - turn_stream = False # Load elevation data to check direction of streams - lake_ele = False # Load elevation for lakes - no_name = False # Do not load SSR place names - no_nve = False # Do not load NVE lake data - no_node = False # Do not merge common nodes at intersections - - # Parse parameters - - if len(sys.argv) < 3: - message ("Please provide 1) municipality, and 2) data category parameter.\n") - message ("Data categories: %s\n" % ", ".join(data_categories)) - message ("Options: -debug, -tag, -geojson, -stream, -ele, -noname, -nonve, -nonode\n\n") - sys.exit() - - # Get municipality - - municipality_query = sys.argv[1] - [municipality_id, municipality_name] = get_municipality_name(municipality_query) - if municipality_id is None: - sys.exit("Municipality '%s' not found\n" % municipality_query) - else: - message ("Municipality:\t%s %s\n" % (municipality_id, municipality_name)) - - # Get N50 data category - - data_category = None - for category in data_categories: - if sys.argv[2].lower() in category.lower(): - data_category = category - message ("N50 category:\t%s\n" % data_category) - break - if not data_category: - sys.exit("Please provide data category: %s\n" % ", ".join(data_categories)) - - # Get other options - - if "-debug" in sys.argv: - debug = True - if "-tag" in sys.argv or "-tags" in sys.argv: - n50_tags = True - if "-geojson" in sys.argv or "-json" in sys.argv: - json_output = True - if "-stream" in sys.argv or "-bekk" in sys.argv: - turn_stream = True - if "-ele" in sys.argv or "-høyde" in sys.argv: - lake_ele = True - if "-noname" in sys.argv: - no_name = True - if "-nonve" in sys.argv: - no_nve = True - if "-nonode" in sys.argv: - no_node = True - - if not turn_stream or not lake_ele: - message ("*** Remember -stream and -ele options before importing.\n") - - output_filename = "n50_%s_%s_%s" % (municipality_id, municipality_name.replace(" ", "_"), data_category) - - # Process data - - if data_category == "BygningerOgAnlegg": - load_building_types() - - load_n50_data(municipality_id, municipality_name, data_category) - - if json_output: - save_geojson(output_filename + ".geojson") - else: - split_polygons() - if data_category == "Arealdekke": - if turn_stream: - fix_stream_direction() # Note: Slow api - if not no_nve: - get_nve_lakes() - find_islands() # Note: "Havflate" is removed at the end of this process - if not no_name: - get_place_names() - match_nodes() - save_osm(output_filename + ".osm") - - duration = time.time() - start_time - message ("\tTotal run time %s (%i features per second)\n\n" % (timeformat(duration), int(len(features) / duration))) +if __name__ == "__main__": + + start_time = time.time() + message("\n-- n50osm v%s --\n" % version) + + features = [] # All geometry and tags + segments = [] # Line segments which are shared by one or more polygons + nodes = ( + set() + ) # Common nodes at intersections, including start/end nodes of segments [lon,lat] + building_tags = {} # Conversion table from building type to osm tag + object_count = {} # Count loaded object types + + debug = False # Include debug tags and unused segments + n50_tags = False # Include property tags from N50 in output + json_output = False # Output complete and unprocessed geometry in geojson format + turn_stream = False # Load elevation data to check direction of streams + lake_ele = False # Load elevation for lakes + no_name = False # Do not load SSR place names + no_nve = False # Do not load NVE lake data + no_node = False # Do not merge common nodes at intersections + + # Parse parameters + + if len(sys.argv) < 3: + message("Please provide 1) municipality, and 2) data category parameter.\n") + message("Data categories: %s\n" % ", ".join(data_categories)) + message( + "Options: -debug, -tag, -geojson, -stream, -ele, -noname, -nonve, -nonode\n\n" + ) + sys.exit() + + # Get municipality + + municipality_query = sys.argv[1] + [municipality_id, municipality_name] = get_municipality_name(municipality_query) + if municipality_id is None: + sys.exit("Municipality '%s' not found\n" % municipality_query) + else: + message("Municipality:\t%s %s\n" % (municipality_id, municipality_name)) + + # Get N50 data category + + data_category = None + for category in data_categories: + if sys.argv[2].lower() in category.lower(): + data_category = category + message("N50 category:\t%s\n" % data_category) + break + if not data_category: + sys.exit("Please provide data category: %s\n" % ", ".join(data_categories)) + + # Get other options + + if "-debug" in sys.argv: + debug = True + if "-tag" in sys.argv or "-tags" in sys.argv: + n50_tags = True + if "-geojson" in sys.argv or "-json" in sys.argv: + json_output = True + if "-stream" in sys.argv or "-bekk" in sys.argv: + turn_stream = True + if "-ele" in sys.argv or "-høyde" in sys.argv: + lake_ele = True + if "-noname" in sys.argv: + no_name = True + if "-nonve" in sys.argv: + no_nve = True + if "-nonode" in sys.argv: + no_node = True + + if not turn_stream or not lake_ele: + message("*** Remember -stream and -ele options before importing.\n") + + output_filename = "n50_%s_%s_%s" % ( + municipality_id, + municipality_name.replace(" ", "_"), + data_category, + ) + + # Process data + + if data_category == "BygningerOgAnlegg": + load_building_types() + + load_n50_data(municipality_id, municipality_name, data_category) + + if json_output: + save_geojson(output_filename + ".geojson") + else: + split_polygons() + if data_category == "Arealdekke": + if turn_stream: + fix_stream_direction() # Note: Slow api + if not no_nve: + get_nve_lakes() + find_islands() # Note: "Havflate" is removed at the end of this process + if not no_name: + get_place_names() + match_nodes() + save_osm(output_filename + ".osm") + + duration = time.time() - start_time + message( + "\tTotal run time %s (%i features per second)\n\n" + % (timeformat(duration), int(len(features) / duration)) + ) diff --git a/utm.py b/utm.py index a7e8202..d0744df 100644 --- a/utm.py +++ b/utm.py @@ -1,4 +1,4 @@ -''' +""" This module converts between lat long and UTM coordinates. Geographic coordinates are entered and displayed in degrees. @@ -17,8 +17,8 @@ Converted from javascript by Nenad Uzunovic Original source http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html -''' - +""" + import math @@ -26,455 +26,476 @@ sm_a = 6378137.0 sm_b = 6356752.314 sm_EccSquared = 6.69437999013e-03 - + UTMScaleFactor = 0.9996 def DegToFloat(degrees, minutes, seconds): - ''' + """ Converts angle in format deg,min,sec to a floating point number - ''' - if (degrees>=0): - return (degrees) + (minutes/60.0) + (seconds/3600.0) + """ + if degrees >= 0: + return (degrees) + (minutes / 60.0) + (seconds / 3600.0) else: - return (degrees) - (minutes/60.0) - (seconds/3600.0) + return (degrees) - (minutes / 60.0) - (seconds / 3600.0) def DegToRad(deg): - ''' + """ Converts degrees to radians. - ''' - return (deg / 180.0 * math.pi) + """ + return deg / 180.0 * math.pi def RadToDeg(rad): - ''' + """ Converts radians to degrees. - ''' - return (rad / math.pi * 180.0) + """ + return rad / math.pi * 180.0 def ArcLengthOfMeridian(phi): - ''' + """ Computes the ellipsoidal distance from the equator to a point at a given latitude. - + Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - + Inputs: phi - Latitude of the point, in radians. - + Globals: sm_a - Ellipsoid model major axis. sm_b - Ellipsoid model minor axis. - + Outputs: The ellipsoidal distance of the point from the equator, in meters. - ''' - + """ + # Precalculate n n = (sm_a - sm_b) / (sm_a + sm_b) - + # Precalculate alpha - alpha = ((sm_a + sm_b) / 2.0) \ - * (1.0 + (math.pow (n, 2.0) / 4.0) + (math.pow (n, 4.0) / 64.0)) - + alpha = ((sm_a + sm_b) / 2.0) * ( + 1.0 + (math.pow(n, 2.0) / 4.0) + (math.pow(n, 4.0) / 64.0) + ) + # Precalculate beta - beta = (-3.0 * n / 2.0) + (9.0 * math.pow (n, 3.0) / 16.0) \ - + (-3.0 * math.pow (n, 5.0) / 32.0) - + beta = ( + (-3.0 * n / 2.0) + + (9.0 * math.pow(n, 3.0) / 16.0) + + (-3.0 * math.pow(n, 5.0) / 32.0) + ) + # Precalculate gamma - gamma = (15.0 * math.pow (n, 2.0) / 16.0) \ - + (-15.0 * math.pow (n, 4.0) / 32.0) - + gamma = (15.0 * math.pow(n, 2.0) / 16.0) + (-15.0 * math.pow(n, 4.0) / 32.0) + # Precalculate delta - delta = (-35.0 * math.pow (n, 3.0) / 48.0) \ - + (105.0 * math.pow (n, 5.0) / 256.0) - + delta = (-35.0 * math.pow(n, 3.0) / 48.0) + (105.0 * math.pow(n, 5.0) / 256.0) + # Precalculate epsilon - epsilon = (315.0 * math.pow (n, 4.0) / 512.0) - + epsilon = 315.0 * math.pow(n, 4.0) / 512.0 + # Now calculate the sum of the series and return - result = alpha \ - * (phi + (beta * math.sin (2.0 * phi)) \ - + (gamma * math.sin (4.0 * phi)) \ - + (delta * math.sin (6.0 * phi)) \ - + (epsilon * math.sin (8.0 * phi))) - + result = alpha * ( + phi + + (beta * math.sin(2.0 * phi)) + + (gamma * math.sin(4.0 * phi)) + + (delta * math.sin(6.0 * phi)) + + (epsilon * math.sin(8.0 * phi)) + ) + return result def UTMCentralMeridian(zone): - ''' + """ Determines the central meridian for the given UTM zone. - + Inputs: zone - An integer value designating the UTM zone, range [1,60]. - + Outputs: The central meridian for the given UTM zone, in radians, or zero if the UTM zone parameter is outside the range [1,60]. Range of the central meridian is the radian equivalent of [-177,+177]. - ''' + """ return DegToRad(-183.0 + (zone * 6.0)) def FootpointLatitude(y): - ''' + """ Computes the footpoint latitude for use in converting transverse Mercator coordinates to ellipsoidal coordinates. - + Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - + Inputs: y - The UTM northing coordinate, in meters. - + Outputs: The footpoint latitude, in radians. - ''' - + """ + # Precalculate n (Eq. 10.18) n = (sm_a - sm_b) / (sm_a + sm_b) - + # Precalculate alpha_ (Eq. 10.22) # (Same as alpha in Eq. 10.17) - alpha_ = ((sm_a + sm_b) / 2.0) \ - * (1 + (math.pow (n, 2.0) / 4) + (math.pow (n, 4.0) / 64)) - + alpha_ = ((sm_a + sm_b) / 2.0) * ( + 1 + (math.pow(n, 2.0) / 4) + (math.pow(n, 4.0) / 64) + ) + # Precalculate y_ (Eq. 10.23) y_ = y / alpha_ - + # Precalculate beta_ (Eq. 10.22) - beta_ = (3.0 * n / 2.0) + (-27.0 * math.pow (n, 3.0) / 32.0) \ - + (269.0 * math.pow (n, 5.0) / 512.0) - + beta_ = ( + (3.0 * n / 2.0) + + (-27.0 * math.pow(n, 3.0) / 32.0) + + (269.0 * math.pow(n, 5.0) / 512.0) + ) + # Precalculate gamma_ (Eq. 10.22) - gamma_ = (21.0 * math.pow (n, 2.0) / 16.0) \ - + (-55.0 * math.pow (n, 4.0) / 32.0) - + gamma_ = (21.0 * math.pow(n, 2.0) / 16.0) + (-55.0 * math.pow(n, 4.0) / 32.0) + # Precalculate delta_ (Eq. 10.22) - delta_ = (151.0 * math.pow (n, 3.0) / 96.0) \ - + (-417.0 * math.pow (n, 5.0) / 128.0) - + delta_ = (151.0 * math.pow(n, 3.0) / 96.0) + (-417.0 * math.pow(n, 5.0) / 128.0) + # Precalculate epsilon_ (Eq. 10.22) - epsilon_ = (1097.0 * math.pow (n, 4.0) / 512.0) - + epsilon_ = 1097.0 * math.pow(n, 4.0) / 512.0 + # Now calculate the sum of the series (Eq. 10.21) - result = y_ + (beta_ * math.sin (2.0 * y_)) \ - + (gamma_ * math.sin (4.0 * y_)) \ - + (delta_ * math.sin (6.0 * y_)) \ - + (epsilon_ * math.sin (8.0 * y_)) - + result = ( + y_ + + (beta_ * math.sin(2.0 * y_)) + + (gamma_ * math.sin(4.0 * y_)) + + (delta_ * math.sin(6.0 * y_)) + + (epsilon_ * math.sin(8.0 * y_)) + ) + return result def MapLatLonToXY(phi, lambda_pt, lambda_ctr): - ''' + """ Converts a latitude/longitude pair to x and y coordinates in the Transverse Mercator projection. Note that Transverse Mercator is not the same as UTM; a scale factor is required to convert between them. - + Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - + Inputs: phi - Latitude of the point, in radians. lambda_pt - Longitude of the point, in radians. lambda_ctr - Longitude of the central meridian to be used, in radians. - + Outputs: xy - A 2-element array containing the x and y coordinates of the computed point. - ''' - + """ + # Precalculate ep2 - ep2 = (math.pow (sm_a, 2.0) - math.pow (sm_b, 2.0)) / math.pow (sm_b, 2.0) - + ep2 = (math.pow(sm_a, 2.0) - math.pow(sm_b, 2.0)) / math.pow(sm_b, 2.0) + # Precalculate nu2 - nu2 = ep2 * math.pow (math.cos (phi), 2.0) - + nu2 = ep2 * math.pow(math.cos(phi), 2.0) + # Precalculate N - N = math.pow (sm_a, 2.0) / (sm_b * math.sqrt (1 + nu2)) - + N = math.pow(sm_a, 2.0) / (sm_b * math.sqrt(1 + nu2)) + # Precalculate t - t = math.tan (phi) + t = math.tan(phi) t2 = t * t # tmp = (t2 * t2 * t2) - math.pow (t, 6.0) - + # Precalculate l l = lambda_pt - lambda_ctr - + # Precalculate coefficients for l**n in the equations below # so a normal human being can read the expressions for easting # and northing # -- l**1 and l**2 have coefficients of 1.0 l3coef = 1.0 - t2 + nu2 - + l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2) - - l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 \ - - 58.0 * t2 * nu2 - - l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 \ - - 330.0 * t2 * nu2 - + + l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2 + + l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2 + l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2) - + l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2) - + # Calculate easting (x) xy = [0.0, 0.0] - xy[0] = N * math.cos (phi) * l \ - + (N / 6.0 * math.pow (math.cos (phi), 3.0) * l3coef * math.pow (l, 3.0)) \ - + (N / 120.0 * math.pow (math.cos (phi), 5.0) * l5coef * math.pow (l, 5.0)) \ - + (N / 5040.0 * math.pow (math.cos (phi), 7.0) * l7coef * math.pow (l, 7.0)) - + xy[0] = ( + N * math.cos(phi) * l + + (N / 6.0 * math.pow(math.cos(phi), 3.0) * l3coef * math.pow(l, 3.0)) + + (N / 120.0 * math.pow(math.cos(phi), 5.0) * l5coef * math.pow(l, 5.0)) + + (N / 5040.0 * math.pow(math.cos(phi), 7.0) * l7coef * math.pow(l, 7.0)) + ) + # Calculate northing (y) - xy[1] = ArcLengthOfMeridian (phi) \ - + (t / 2.0 * N * math.pow (math.cos (phi), 2.0) * math.pow (l, 2.0)) \ - + (t / 24.0 * N * math.pow (math.cos (phi), 4.0) * l4coef * math.pow (l, 4.0)) \ - + (t / 720.0 * N * math.pow (math.cos (phi), 6.0) * l6coef * math.pow (l, 6.0)) \ - + (t / 40320.0 * N * math.pow (math.cos (phi), 8.0) * l8coef * math.pow (l, 8.0)) - + xy[1] = ( + ArcLengthOfMeridian(phi) + + (t / 2.0 * N * math.pow(math.cos(phi), 2.0) * math.pow(l, 2.0)) + + (t / 24.0 * N * math.pow(math.cos(phi), 4.0) * l4coef * math.pow(l, 4.0)) + + (t / 720.0 * N * math.pow(math.cos(phi), 6.0) * l6coef * math.pow(l, 6.0)) + + (t / 40320.0 * N * math.pow(math.cos(phi), 8.0) * l8coef * math.pow(l, 8.0)) + ) + return xy def MapXYToLatLon(x, y, lambda_ctr): - ''' + """ Converts x and y coordinates in the Transverse Mercator projection to a latitude/longitude pair. Note that Transverse Mercator is not the same as UTM; a scale factor is required to convert between them. - + Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. - + Inputs: x - The easting of the point, in meters. y - The northing of the point, in meters. lambda_ctr - Longitude of the central meridian to be used, in radians. - + Outputs: philambda - A 2-element containing the latitude and longitude in radians. - + Remarks: The local variables Nf, nuf2, tf, and tf2 serve the same purpose as N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect to the footpoint latitude phif. - + x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and to optimize computations. - ''' - + """ + # Get the value of phif, the footpoint latitude. - phif = FootpointLatitude (y) - + phif = FootpointLatitude(y) + # Precalculate ep2 - ep2 = (math.pow (sm_a, 2.0) - math.pow (sm_b, 2.0)) \ - / math.pow (sm_b, 2.0) - + ep2 = (math.pow(sm_a, 2.0) - math.pow(sm_b, 2.0)) / math.pow(sm_b, 2.0) + # Precalculate cos (phif) - cf = math.cos (phif) - + cf = math.cos(phif) + # Precalculate nuf2 - nuf2 = ep2 * math.pow (cf, 2.0) - + nuf2 = ep2 * math.pow(cf, 2.0) + # Precalculate Nf and initialize Nfpow - Nf = math.pow (sm_a, 2.0) / (sm_b * math.sqrt (1 + nuf2)) + Nf = math.pow(sm_a, 2.0) / (sm_b * math.sqrt(1 + nuf2)) Nfpow = Nf - + # Precalculate tf - tf = math.tan (phif) + tf = math.tan(phif) tf2 = tf * tf tf4 = tf2 * tf2 - + # Precalculate fractional coefficients for x**n in the equations # below to simplify the expressions for latitude and longitude. x1frac = 1.0 / (Nfpow * cf) - - Nfpow *= Nf # now equals Nf**2) + + Nfpow *= Nf # now equals Nf**2) x2frac = tf / (2.0 * Nfpow) - - Nfpow *= Nf # now equals Nf**3) + + Nfpow *= Nf # now equals Nf**3) x3frac = 1.0 / (6.0 * Nfpow * cf) - - Nfpow *= Nf # now equals Nf**4) + + Nfpow *= Nf # now equals Nf**4) x4frac = tf / (24.0 * Nfpow) - - Nfpow *= Nf # now equals Nf**5) + + Nfpow *= Nf # now equals Nf**5) x5frac = 1.0 / (120.0 * Nfpow * cf) - - Nfpow *= Nf # now equals Nf**6) + + Nfpow *= Nf # now equals Nf**6) x6frac = tf / (720.0 * Nfpow) - - Nfpow *= Nf # now equals Nf**7) + + Nfpow *= Nf # now equals Nf**7) x7frac = 1.0 / (5040.0 * Nfpow * cf) - - Nfpow *= Nf # now equals Nf**8) + + Nfpow *= Nf # now equals Nf**8) x8frac = tf / (40320.0 * Nfpow) - + # Precalculate polynomial coefficients for x**n. # -- x**1 does not have a polynomial coefficient. x2poly = -1.0 - nuf2 - + x3poly = -1.0 - 2 * tf2 - nuf2 - - x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 \ - - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2) - + + x4poly = ( + 5.0 + + 3.0 * tf2 + + 6.0 * nuf2 + - 6.0 * tf2 * nuf2 + - 3.0 * (nuf2 * nuf2) + - 9.0 * tf2 * (nuf2 * nuf2) + ) + x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2 - - x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 \ - + 162.0 * tf2 * nuf2 - + + x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2 + x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2) - + x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2) - + # Calculate latitude philambda = [0.0, 0.0] - philambda[0] = phif + x2frac * x2poly * (x * x) \ - + x4frac * x4poly * math.pow (x, 4.0) \ - + x6frac * x6poly * math.pow (x, 6.0) \ - + x8frac * x8poly * math.pow (x, 8.0) - + philambda[0] = ( + phif + + x2frac * x2poly * (x * x) + + x4frac * x4poly * math.pow(x, 4.0) + + x6frac * x6poly * math.pow(x, 6.0) + + x8frac * x8poly * math.pow(x, 8.0) + ) + # Calculate longitude - philambda[1] = lambda_ctr + x1frac * x \ - + x3frac * x3poly * math.pow (x, 3.0) \ - + x5frac * x5poly * math.pow (x, 5.0) \ - + x7frac * x7poly * math.pow (x, 7.0) - + philambda[1] = ( + lambda_ctr + + x1frac * x + + x3frac * x3poly * math.pow(x, 3.0) + + x5frac * x5poly * math.pow(x, 5.0) + + x7frac * x7poly * math.pow(x, 7.0) + ) + return philambda def LatLonToUTMXY(lat, lon, zone): - ''' + """ Converts a latitude/longitude pair to x and y coordinates in the Universal Transverse Mercator projection. - + Inputs: lat - Latitude of the point, in radians. lon - Longitude of the point, in radians. zone - UTM zone to be used for calculating values for x and y. If zone is less than 1 or greater than 60, the routine will determine the appropriate zone from the value of lon. - + Outputs: xy - A 2-element array where the UTM x and y values will be stored. - ''' - + """ + xy = MapLatLonToXY(lat, lon, UTMCentralMeridian(zone)) - + # Adjust easting and northing for UTM system. xy[0] = xy[0] * UTMScaleFactor + 500000.0 xy[1] = xy[1] * UTMScaleFactor - if (xy[1] < 0.0): + if xy[1] < 0.0: xy[1] = xy[1] + 10000000.0 - + return xy def UTMXYToLatLon(x, y, zone, southhemi): - ''' + """ Converts x and y coordinates in the Universal Transverse Mercator projection to a latitude/longitude pair. - + Inputs: x - The easting of the point, in meters. y - The northing of the point, in meters. zone - The UTM zone in which the point lies. southhemi - True if the point is in the southern hemisphere; false otherwise. - + Outputs: latlon - A 2-element array containing the latitude and longitude of the point, in radians. - ''' + """ x -= 500000.0 x /= UTMScaleFactor - + # If in southern hemisphere, adjust y accordingly. - if (southhemi): + if southhemi: y -= 10000000.0 - + y /= UTMScaleFactor - + cmeridian = UTMCentralMeridian(zone) latlon = MapXYToLatLon(x, y, cmeridian) - + return latlon def LatLonToUtm(lat, lon): - ''' + """ Converts lat lon to utm - + Inputs: lat - lattitude in degrees lon - longitude in degrees - + Outputs: xy - utm x(easting), y(northing) zone - utm zone hemi - 'N' or 'S' - ''' - - if ((lon < -180.0) or (180.0 <= lon)): - print ('The longitude you entered is out of range -', lon) - print ('Please enter a number in the range [-180, 180).') + """ + + if (lon < -180.0) or (180.0 <= lon): + print("The longitude you entered is out of range -", lon) + print("Please enter a number in the range [-180, 180).") return 0 - - if ((lat < -90.0) or (90.0 < lat)): - print ('The latitude you entered is out of range -', lat) - print ('Please enter a number in the range [-90, 90].') - + + if (lat < -90.0) or (90.0 < lat): + print("The latitude you entered is out of range -", lat) + print("Please enter a number in the range [-90, 90].") + # Compute the UTM zone. - zone = math.floor ((lon + 180.0) / 6) + 1 - + zone = math.floor((lon + 180.0) / 6) + 1 + # Convert - xy = LatLonToUTMXY (DegToRad(lat), DegToRad(lon), zone) - + xy = LatLonToUTMXY(DegToRad(lat), DegToRad(lon), zone) + # Determine hemisphere - hemi = 'N' - if (lat < 0): - hemi = 'S' - + hemi = "N" + if lat < 0: + hemi = "S" + return [xy, zone, hemi] def UtmToLatLon(x, y, zone, hemi): - ''' + """ Converts UTM coordinates to lat long - + Inputs: x - easting (in meters) y - northing (in meters) zone - UTM zone hemi - 'N' or 'S' - + Outputs: latlong - [lattitude, longitude] (in degrees) - ''' - if ((zone < 1) or (60 < zone)): - print ('The UTM zone you entered is out of range -', zone) - print ('Please enter a number in the range [1, 60].') + """ + if (zone < 1) or (60 < zone): + print("The UTM zone you entered is out of range -", zone) + print("Please enter a number in the range [1, 60].") return 0 - - if ((hemi != 'N') and (hemi != 'S')): - print ('The hemisphere you entered is wrong -', hemi) - print ('Please enter N or S') - + + if (hemi != "N") and (hemi != "S"): + print("The hemisphere you entered is wrong -", hemi) + print("Please enter N or S") + southhemi = False - if (hemi == 'S'): + if hemi == "S": southhemi = True - + # Convert latlon = UTMXYToLatLon(x, y, zone, southhemi) - + # Convert to degrees latlon[0] = RadToDeg(latlon[0]) latlon[1] = RadToDeg(latlon[1]) - + return latlon From 5d051d7ae0ff7154fa9d4e65804dcfca7536568e Mon Sep 17 00:00:00 2001 From: Harald Husum Date: Fri, 16 Sep 2022 20:12:05 +0200 Subject: [PATCH 2/2] Format with black --preview --- n50merge.py | 24 +++++--------------- n50osm.py | 64 ++++++----------------------------------------------- 2 files changed, 12 insertions(+), 76 deletions(-) diff --git a/n50merge.py b/n50merge.py index 2bfa3c3..95fffbc 100644 --- a/n50merge.py +++ b/n50merge.py @@ -17,7 +17,9 @@ overpass_api = "https://overpass-api.de/api/interpreter" # Overpass endpoint -import_folder = "~/Jottacloud/osm/n50/" # Folder containing import highway files (default folder tried first) +import_folder = ( # Folder containing import highway files (default folder tried first) + "~/Jottacloud/osm/n50/" +) n50_parts = ["coastline", "water", "wood", "landuse"] @@ -31,7 +33,6 @@ def message(text): - sys.stderr.write(text) sys.stderr.flush() @@ -40,7 +41,6 @@ def message(text): def timeformat(sec): - if sec > 3600: return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) elif sec > 60: @@ -53,7 +53,6 @@ def timeformat(sec): def get_municipality_name(query): - if query.isdigit(): url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query else: @@ -102,7 +101,6 @@ def get_municipality_name(query): def prepare_data(root, tree, nodes, ways, relations): - # Create dict of nodes # Each node is a tuple of (lon, lat), corresponding to GeoJSON format x,y @@ -179,7 +177,6 @@ def prepare_data(root, tree, nodes, ways, relations): def load_osm(): - global osm_root, osm_tree message("Load existing OSM elements from Overpass ...\n") @@ -187,7 +184,8 @@ def load_osm(): query = ( "[timeout:90];(area[ref=%s][admin_level=7][place=municipality];)->.a;" % municipality_id - + '( nwr["natural"](area.a); nwr["waterway"](area.a); nwr["landuse"](area.a); nwr["leisure"](area.a);' + + '( nwr["natural"](area.a); nwr["waterway"](area.a); nwr["landuse"](area.a);' + ' nwr["leisure"](area.a);' + 'nwr["aeroway"](area.a); nwr["seamark:type"="rock"](area.a); );' + "(._;>;<;);out meta;" ) @@ -230,7 +228,6 @@ def load_osm(): def load_n50(): - global n50_root, n50_tree message("Load N50 import file elements ...\n") @@ -260,9 +257,7 @@ def load_n50(): def filter_parts(part, n50_elements, found_elements): - for element_id, element in iter(n50_elements.items()): - all_tags = element["xml"].findall("tag") if all_tags != None: for tag in all_tags: @@ -295,7 +290,6 @@ def filter_parts(part, n50_elements, found_elements): def split_n50(): - message("Splitting N50 import file ...\n") for part in n50_parts: @@ -325,7 +319,6 @@ def split_n50(): if part in ["coastline", "water"]: for relation_id, relation in iter(n50_relations.items()): if relation["island"]: - for member in relation["members"]: if member in ways: relations.add(relation_id) @@ -372,7 +365,6 @@ def split_n50(): def merge_osm(): - message("Merge OSM ways ...\n") count = 0 @@ -406,7 +398,6 @@ def merge_osm(): def merge_tags(n50_xml, osm_xml): - n50_tags = {} for tag in n50_xml.findall("tag"): n50_tags[tag.attrib["k"]] = tag.attrib["v"] @@ -434,7 +425,6 @@ def merge_tags(n50_xml, osm_xml): def merge_n50(): - message("Merge N50 with OSM ...\n") lap_time = time.time() @@ -455,11 +445,9 @@ def merge_n50(): and not n50_way["incomplete"] and not osm_way["incomplete"] ): - # Check if conflicting tags, for example natural=water + natural=wood. Also update OSM tags. if merge_tags(n50_way["xml"], osm_way["xml"]): - # Swap ref if member in relations for parent_id in n50_way["parents"]: @@ -574,7 +562,6 @@ def indent_tree(elem, level=0): def save_osm(): - message("Saving file ...\n") # Merge remaining N50 tree into OSM tree @@ -596,7 +583,6 @@ def save_osm(): # Main program if __name__ == "__main__": - start_time = time.time() message("\nn50merge v%s\n\n" % version) diff --git a/n50osm.py b/n50osm.py index 6aa4969..8877d38 100644 --- a/n50osm.py +++ b/n50osm.py @@ -146,7 +146,6 @@ def tag_object(feature_type, geometry_type, properties, feature): - tags = {} missing_tags = set() @@ -276,7 +275,6 @@ def tag_object(feature_type, geometry_type, properties, feature): def message(output_text): - sys.stdout.write(output_text) sys.stdout.flush() @@ -285,7 +283,6 @@ def message(output_text): def timeformat(sec): - if sec > 3600: return "%i:%02i:%02i hours" % (sec / 3600, (sec % 3600) / 60, sec % 60) elif sec > 60: @@ -302,7 +299,6 @@ def timeformat(sec): def polygon_area(polygon): - if polygon[0] == polygon[-1]: lat_dist = math.pi * 6371009.0 / 180.0 @@ -327,14 +323,12 @@ def polygon_area(polygon): def multipolygon_area(multipolygon): - if ( type(multipolygon) is list and len(multipolygon) > 0 and type(multipolygon[0]) is list and multipolygon[0][0] == multipolygon[0][-1] ): - area = polygon_area(multipolygon[0]) for patch in multipolygon[1:]: inner_area = polygon_area(patch) @@ -353,7 +347,6 @@ def multipolygon_area(multipolygon): def polygon_centroid(polygon): - if polygon[0] == polygon[-1]: x = 0 y = 0 @@ -376,7 +369,6 @@ def polygon_centroid(polygon): def inside_polygon(point, polygon): - if polygon[0] == polygon[-1]: x, y = point n = len(polygon) @@ -404,14 +396,12 @@ def inside_polygon(point, polygon): def inside_multipolygon(point, multipolygon): - if ( type(multipolygon) is list and len(multipolygon) > 0 and type(multipolygon[0]) is list and multipolygon[0][0] == multipolygon[0][-1] ): - inside = inside_polygon(point, multipolygon[0]) if inside: for patch in multipolygon[1:]: @@ -430,7 +420,6 @@ def inside_multipolygon(point, multipolygon): def coordinate_offset(node, distance): - m = 1 / ((math.pi / 180.0) * 6378137.0) # Degrees per meter latitude = node[1] + (distance * m) @@ -444,7 +433,6 @@ def coordinate_offset(node, distance): def get_bbox(coordinates, perimeter): - if type(coordinates) is tuple: patch = [coordinates] elif type(coordinates[0]) is tuple: @@ -471,7 +459,6 @@ def get_bbox(coordinates, perimeter): def create_point(node, gml_id, note): - if debug: entry = { "object": "Debug", @@ -492,7 +479,6 @@ def create_point(node, gml_id, note): def parse_coordinates(coord_text): - global gml_id parse_count = 0 @@ -544,7 +530,6 @@ def parse_coordinates(coord_text): def get_municipality_name(query): - if query.isdigit(): url = "https://ws.geonorge.no/kommuneinfo/v1/kommuner/" + query else: @@ -593,7 +578,6 @@ def get_municipality_name(query): def load_building_types(): - # file = open("building_types.csv") url = "https://raw.githubusercontent.com/NKAmapper/building2osm/main/building_types.csv" request = urllib.request.Request(url, headers=header) @@ -626,7 +610,6 @@ def load_building_types(): def simple_length(coord): - length = 0 for i in range(len(coord) - 2): length += (coord[i + 1][0] - coord[i][0]) ** 2 + ( @@ -640,7 +623,6 @@ def simple_length(coord): def split_patch(coordinates): - for i in range(1, len(coordinates) - 1): first = coordinates.index(coordinates[i]) if first < i: @@ -662,7 +644,6 @@ def split_patch(coordinates): def get_property(top, ns_app): - properties = {} if ns_app in top.tag: tag = top.tag[len(ns_app) + 2 :] @@ -681,7 +662,6 @@ def get_property(top, ns_app): def load_n50_data(municipality_id, municipality_name, data_category): - global gml_id lap = time.time() @@ -765,7 +745,6 @@ def load_n50_data(municipality_id, municipality_name, data_category): properties = {} # Attributes provided from GML for app in feature[0]: - # Get app properties/attributes tag = app.tag[len(ns_app) + 2 :] @@ -913,7 +892,7 @@ def load_n50_data(municipality_id, municipality_name, data_category): message("\t\t%i\t%s\n" % (object_count[object_type], object_type)) if missing_tags: - message("\tNot tagged: %s\n" % (", ".join(missing_tags.difference("Havflate")))) + message("\tNot tagged: %s\n" % ", ".join(missing_tags.difference("Havflate"))) if stream_count > 0: message("\t%i streams\n" % stream_count) @@ -927,7 +906,6 @@ def load_n50_data(municipality_id, municipality_name, data_category): def create_border_segments(patch, members, gml_id, match): - # First create list of existing conncetions between coordinates i and i+1 of patch connection = [] @@ -956,7 +934,6 @@ def create_border_segments(patch, members, gml_id, match): start_index = 0 while start_index < len(patch) - 2: - while start_index < len(patch) - 1 and connection[start_index]: start_index += 1 @@ -985,7 +962,6 @@ def create_border_segments(patch, members, gml_id, match): def segment_position(segment_index, patch): - return patch.index(segments[segment_index]["coordinates"][1]) @@ -997,7 +973,6 @@ def segment_position(segment_index, patch): def split_polygons(): - message("Decompose polygons into segments...\n") # Create bbox for segments and line features @@ -1011,7 +986,6 @@ def split_polygons(): split_count = 0 for feature in features: - if feature["type"] == "Polygon": matching_polygon = [] @@ -1025,7 +999,6 @@ def split_polygons(): [patch_min_bbox, patch_max_bbox] = get_bbox(patch, 0) for i, segment in enumerate(segments): - if ( patch_min_bbox[0] <= segment["max_bbox"][0] and patch_max_bbox[0] >= segment["min_bbox"][0] @@ -1033,7 +1006,6 @@ def split_polygons(): and patch_max_bbox[1] >= segment["min_bbox"][1] and set(segment["coordinates"]) <= patch_set ): - # Note: If patch is a closed way, segment may wrap start/end of patch if len(segment["coordinates"]) == 2: @@ -1070,7 +1042,6 @@ def split_polygons(): "KantUtsnitt", ] ): - segment["used"] += 1 node1 = patch.index(segment["coordinates"][0]) node2 = patch.index(segment["coordinates"][1]) @@ -1120,7 +1091,6 @@ def split_polygons(): def find_islands(): - message("Find islands...\n") lap = time.time() @@ -1168,7 +1138,6 @@ def find_islands(): "FerskvannTørrfall", ]: for i in range(1, len(feature["members"])): - # Do not use patch with intermittent edge found = True @@ -1300,7 +1269,6 @@ def find_islands(): # Add island to features list if closed chain of ways if first_node == last_node: - members = [] coordinates = [first_node] for segment in island: @@ -1362,7 +1330,6 @@ def find_islands(): def match_nodes(): - global delete_count message("Identify intersections...\n") @@ -1384,7 +1351,6 @@ def match_nodes(): nodes.add(feature["coordinates"][-1]) if not no_node: - # Create bbox for segments and line features + create temporary list of LineString features for segment in segments: @@ -1407,7 +1373,6 @@ def match_nodes(): and feature["min_bbox"][1] <= segment["max_bbox"][1] and feature["max_bbox"][1] >= segment["min_bbox"][1] ): - intersections = set(feature["coordinates"]).intersection( set(segment["coordinates"]) ) @@ -1488,7 +1453,6 @@ def match_nodes(): def get_elevation_old(node): - global elevations, ele_count, retry_count ns_ows = "http://www.opengis.net/ows/1.1" @@ -1562,7 +1526,6 @@ def get_elevation_old(node): def get_elevation(node): - global elevations, ele_count, retry_count url = ( @@ -1598,7 +1561,6 @@ def get_elevation(node): def fix_stream_direction(): - global elevations, ele_count, retry_count elevations = {} # Already fetched elevations from api @@ -1668,7 +1630,6 @@ def fix_stream_direction(): def get_ssr_name(feature, name_categories): - global ssr_places, name_count if feature["type"] == "Point": @@ -1707,7 +1668,6 @@ def get_ssr_name(feature, name_categories): # Sort and select name if found_names: - found_names.sort( key=lambda name: name_categories.index(name["tags"]["ssr:type"]) ) # place_name_rank(name, name_categories)) @@ -1754,7 +1714,6 @@ def get_ssr_name(feature, name_categories): or "holme" in found_names[0]["tags"]["ssr:type"] and "holme" not in found_names[1]["tags"]["ssr:type"] ): - if "name" in feature["tags"] and ( len(alt_names) > 1 or feature["tags"]["name"] @@ -1781,7 +1740,6 @@ def get_ssr_name(feature, name_categories): def get_category_place_names(n50_categories, ssr_categories): - for feature in features: if feature["object"] in n50_categories: get_ssr_name(feature, ssr_categories) @@ -1792,7 +1750,6 @@ def get_category_place_names(n50_categories, ssr_categories): def get_place_names(): - global ssr_places, name_count global elevations, ele_count, retry_count @@ -1828,9 +1785,10 @@ def get_place_names(): if "name" in tag.get("k") or tag.get("k") in ["ssr:stedsnr", "ssr:type"]: if tag.get("k") == "fixme": if "multiple name" in tag.get("v"): - entry["tags"][ - "fixme" - ] = "Multiple name tags, please choose one and add the other to alt_name" + entry["tags"]["fixme"] = ( + "Multiple name tags, please choose one and add the other to" + " alt_name" + ) else: entry["tags"][tag.get("k")] = tag.get("v") @@ -1878,7 +1836,6 @@ def get_place_names(): for feature in features: if feature["object"] in ["Innsjø", "InnsjøRegulert"]: - lake_node = get_ssr_name(feature, name_category) area = abs(multipolygon_area(feature["coordinates"])) feature["extras"]["areal"] = str(int(area)) @@ -1890,7 +1847,6 @@ def get_place_names(): and "ele" not in feature["tags"] and (lake_node or area >= lake_ele_size) ): - # Check that name coordinate is not on lake's island if lake_node: if not inside_multipolygon(lake_node, feature["coordinates"]): @@ -1953,7 +1909,6 @@ def get_place_names(): def get_nve_lakes(): - message("Load lake data from NVE...\n") n50_lake_count = 0 @@ -1964,7 +1919,6 @@ def get_nve_lakes(): # Paging results (default 1000 lakes) while more_lakes: - # Alternative url: # url = "https://gis3.nve.no/map/rest/services/Innsjodatabase2/MapServer/find?" + \ # "searchText=%s&contains=true&searchFields=kommNr&layers=5&returnGeometry=false&returnUnformattedValues=true&f=pjson&resultOffset=%i&resultRecordCount=1000&orderByFields=areal_km2%20DESC" \ @@ -2031,7 +1985,6 @@ def get_nve_lakes(): def save_geojson(filename): - message("Save to '%s' file...\n" % filename) json_features = {"type": "FeatureCollection", "features": []} @@ -2083,7 +2036,6 @@ def indent_tree(elem, level=0): def save_osm(filename): - message("Save to '%s' file...\n" % filename) osm_node_ids = {} # Will contain osm_id of each common node @@ -2147,7 +2099,6 @@ def save_osm(filename): # The main objects for feature in features: - if feature["object"] == "Havflate": continue @@ -2187,7 +2138,6 @@ def save_osm(filename): osm_feature.append(osm_nd) elif feature["type"] == "Polygon": - # Output way if possible to avoid relation if ( len(feature["members"]) == 1 @@ -2246,7 +2196,6 @@ def save_osm(filename): # Main program if __name__ == "__main__": - start_time = time.time() message("\n-- n50osm v%s --\n" % version) @@ -2273,7 +2222,8 @@ def save_osm(filename): message("Please provide 1) municipality, and 2) data category parameter.\n") message("Data categories: %s\n" % ", ".join(data_categories)) message( - "Options: -debug, -tag, -geojson, -stream, -ele, -noname, -nonve, -nonode\n\n" + "Options: -debug, -tag, -geojson, -stream, -ele, -noname, -nonve," + " -nonode\n\n" ) sys.exit()