import osmdigest.digest as digest

"""
Explanation: Buildings / Addresses
I want to understand how buildings and addresses are represented in OSM data.
Required reading: http://wiki.openstreetmap.org/wiki/Addresses
End of explanation
"""


import os
#filename = os.path.join("//media", "disk", "OSM_Data", "isle-of-wight-latest.osm.xz")
filename = os.path.join("..", "..", "..", "Data", "isle-of-wight-latest.osm.xz")

building_node_ids = []
addr_node_ids = []

for x in digest.parse(filename):
    if isinstance(x, digest.Node):
        if "building" in x.tags:
            building_node_ids.append(x)
        if any(key.startswith("addr:") for key in x.tags):
            addr_node_ids.append(x)

len(building_node_ids), building_node_ids[:5]

len(addr_node_ids), addr_node_ids[:5]

"""
Explanation: Nodes tagged as buildings / with addresses
For the Isle of Wight dataset, there are rather few nodes tagged as buildings.  There are quite a few more nodes which have address information, however.
End of explanation
"""


building_way_ids = []
addr_way_ids = []

for x in digest.parse(filename):
    if isinstance(x, digest.Way):
        if "building" in x.tags:
            building_way_ids.append(x)
        if any(key.startswith("addr:") for key in x.tags):
            addr_way_ids.append(x)

len(building_way_ids), building_way_ids[:5]

len(addr_way_ids), addr_way_ids[:5]

"""
Explanation: Same for ways
We expected buildings to mostly have a polygonal outline, and so be "ways".
The situation here is reversed: lots of buildings, and fewer addresses.  From eyeballing a few ways which have address information, but which are not buildings, we find that the "way" gives the total outline of, say, a school, which may contain a number of buildings, playing fields etc.
End of explanation
"""


building_rel_ids = []
addr_rel_ids = []

for x in digest.parse(filename):
    if isinstance(x, digest.Relation):
        if "building" in x.tags:
            building_rel_ids.append(x)
        if any(key.startswith("addr:") for key in x.tags):
            addr_rel_ids.append(x)

len(building_rel_ids), building_rel_ids[:5]

len(addr_rel_ids), addr_rel_ids[:5]

"""
Explanation: Finally for relations
Rather few.  These tend to be special cases: usually when the building is non-convex (say a stately home, with an inner courtyard) so a relation is required to specify the "inner" and "outer" ways.
End of explanation
"""


import numpy as np
import pandas as pd

gen = digest.parse(filename)
print(next(gen))
print(next(gen))

possible_address_tags = set()
for x in gen:
    for key in x.tags:
        if key.startswith("addr:"):
            possible_address_tags.add(key)
            
possible_address_tags

gen = digest.parse(filename)
osm = next(gen)
bounds = next(gen)

address_data = { key : [] for key in possible_address_tags }
address_data["osm_id"] = []

for x in gen:
    addr = {key : x.tags[key] for key in x.tags if key.startswith("addr:")}
    if len(addr) > 0:
        address_data["osm_id"].append(x.name+"/"+str(x.osm_id))
        for key in possible_address_tags:
            if key in addr:
                address_data[key].append(addr[key])
            else:
                address_data[key].append(np.nan)

data = pd.DataFrame(address_data)
data = data.set_index("osm_id")
data[:5]

"""
Explanation: Process to a pandas dataframe
In a separate notebook, we explore how to use geopandas to also store and manipulate geometry...
End of explanation
"""


import osmdigest.sqlite as sq
import os

filename = os.path.join("//tmp", "aaa", "illinois-latest.db")
#filename = os.path.join("..", "..", "..", "Data", "illinois-latest.db")

db = sq.OSM_SQLite(filename)

"""
Explanation: Do the same for USA / Chicago data
Let's look at a different country, namely the state of Illinois in the USA.

Download the Illinois data, in osm.bz2 format, from http://download.geofabrik.de/north-america/us/illinois.html
Run the script convert_to_db.py on this file, generating "illinois-latest.db"

We'll then load the data from the SQLite database, as this is far more memory efficient.  We'll find that we still end up with massive python data structures...
The overall pattern is the same, even though the dataset is much larger.
- Most buildings are "ways" and most addresses are "ways".
- There are a few point addresses stored as nodes, and complicated structures (like an airforce base) are stored as relations.
End of explanation
"""


def iterate_over_tags(iterator):
    buildings, addresses = [], []
    for element in iterator:
        if any(key.startswith("building") for key in element.tags):
            buildings.append(element)
        if any(key.startswith("addr") for key in element.tags):
            addresses.append(element)
    return buildings, addresses

building_nodes, address_nodes = iterate_over_tags(db.nodes())

len(building_nodes), building_nodes[:5]

len(address_nodes), address_nodes[:5]

"""
Explanation: Nodes
End of explanation
"""


building_ways, address_ways = iterate_over_tags(db.ways())

len(building_ways), building_ways[:5]

len(address_ways), address_ways[:5]

"""
Explanation: Ways
End of explanation
"""


building_rels, address_rels = iterate_over_tags(db.relations())

len(building_rels), building_rels[:5]

len(address_rels), address_rels[:5]

"""
Explanation: Relations
End of explanation
"""


features = []

def make_feature(el, centroid):
    return { "properties": {
                    "street": el.tags["addr:street"],
                    "housenumber": el.tags["addr:housenumber"],
                    "osm_id": "{}/{}".format(el.name, el.osm_id)
                },
                "geometry": { "type": "Point",
                    "coordinates": centroid } }

for el in db.search_node_tag_keys({"addr:street", "addr:housenumber"}):
    features.append(make_feature(el, [el.longitude, el.latitude]))
    
for el in db.search_way_tag_keys({"addr:street", "addr:housenumber"}):
    way = db.complete_way(el)
    features.append(make_feature(el, way.centroid()))

for el in db.search_relation_tag_keys({"addr:street", "addr:housenumber"}):
    rel = db.complete_relation(el)
    features.append(make_feature(el, rel.centroid()))

import geopandas as gpd

frame = gpd.GeoDataFrame.from_features(features)
#frame = frame.set_geometry("centroid")
frame[:5]

"""
Explanation: Process into a useful form
Ultimately, I want:

Elements which have a street name and a housenumber.  These are addr:street and addr:housenumber.
A single coordinate to use.  For a way we can take the centroid.  For a relation, we have to perform a bit of exploring: it is tempting to hope there is always one way with role "outer".  This doesn't work, so we settle on looking at each member and taking its centroid, working recursively if necessary.  This will not be particularly meaningful in some situations, but should suffice for now.

We construct a GeoJSON-like data structure, and then import it into geoPandas.
End of explanation
"""


unexpected_addresses = frame[~ frame.housenumber.map(lambda x : all(y>='0' and y<='9' for y in x))]
unexpected_addresses.head()

"""
Explanation: Explore inconsistencies
Many of the "housenumber" values are not integers.

From the 5 examples below, "2211 W" should be "2211"; "145 N." should be "145" but the street should be "North Curtis Avenue".  For the remaining 3, I do not know enough about North American addresses to be certain.
All these examples come from areas of data where we don't have house-level building data.  So in many ways, the data is not useful for my application anyway.
End of explanation
"""


unexpected_addresses[5:10]

"""
Explanation: The addresses with "#" in appear to be neighbours in a "trailer park", judging from satelite images.  This addressing scheme follows how apartments/flats are addressed.
The address "14543-59" is a church, i.e. a large building which probably occupied more than one notional address.
End of explanation
"""


import re

one = re.compile("^\\d+\\s+#\\d+$")
assert one.match("4262 #12") is not None
assert one.match("4 62 #12") is None
two = re.compile("^\\d+[NSEW]\\d+$")
assert two.match("19N479") is not None
assert two.match("19NS479") is None
assert two.match("19 479") is None
three = re.compile("^\\d+\\s*[a-zA-Z]\\.*$")
assert three.match("152A") is not None
assert three.match("152 A") is not None
assert three.match("152 A.") is not None
assert three.match("152Ac") is None
four = re.compile("^\\d+\\s*1/2")

matches = {one, two, three, four}

left = unexpected_addresses[~ unexpected_addresses.housenumber.map(lambda x : any(m.match(x) is not None for m in matches))]
left.head()

"""
Explanation: So let's filter out things which match.  (Stand-back, I know regular expressions).
- digits, space, #, digits  "\d+\s+#\d+"
- digits, one of NSEW, digits "\d+[NSEW]\d+"
- digits, maybe a space, a single letter. maybe a full stop "\d+\s[a-zA-Z]\."
- digits, maybe a space, "East"/"West" etc.
- things like "123 1/2"
End of explanation
"""


for way in db.search_way_tag_keys({"addr:interpolation"}):
    print(way)

"""
Explanation: This leaves a lot left over.
What interests me is whether there are any "way"s which "interpolate" addresses, see http://wiki.openstreetmap.org/wiki/Addresses
The answer is: "no".  Eyeballing these few in OSM doesn't show anything interesting.
End of explanation
"""


import fiona
fiona.supported_drivers

frame.to_file("illinois_building.json", driver="GeoJSON")

"""
Explanation: Finally, save the data out.  We can use any supported fiona driver.  Here I use GeoJSON, as it's human readable, and no less space efficient than a Shapefile.  A Shapefile can be imported into QGis etc., of course.
End of explanation
"""


