#!/bin/env python

# Extract Bhutan L3s (Chiwogs) from OSM & write as Sahana import file
# Designed to be run in Web2Py context with the BT_L1.csv & BT_L2.csv already imported

# L1 (Dzonkhags) are admin_level = '4' in OSM
# L2 (Gewogs) are admin_level = '6' in OSM
# L3 (Chiwogs) are admin_level = '8' in OSM

from json import loads as json_loads

from pydriosm.reader import GeofabrikReader

from shapely.geometry import point
from shapely.geometry import shape as geojson_loads
from shapely.wkt import dumps as wkt_dumps
from shapely.wkt import loads as wkt_loads

geofabrik_reader = GeofabrikReader(data_dir = None)
# NB This takes a while!
pbf_raw = geofabrik_reader.read_osm_pbf("Bhutan")
multipolygons = pbf_raw["multipolygons"]["multipolygons"]

db = current.db
table = current.s3db.gis_location
base_query = (table.level.belongs(("L1", "L2"))) & \
             (table.deleted == False)

output = "BT_L3.csv"
outputFile = open(output, "w", encoding="utf-8")
write = outputFile.write
write("Country,L1,L2,L3,WKT\n")
for feature in multipolygons:
    feature = json_loads(feature)
    p = feature["properties"]
    if p["admin_level"] != "8":
        continue
    #name = '"%s"' % p["name"].split(" Chiwogs")[0].split(" Chiwog")[0]
    name = p["name"].split(" Chiwogs")[0].split(" Chiwog")[0]
    geom = geojson_loads(feature["geometry"])
    centroid = geom.centroid
    lat = centroid.y
    lon = centroid.x
    test = point.Point(lon, lat)
    query = base_query & (table.lat_min < lat) & \
                         (table.lat_max > lat) & \
                         (table.lon_min < lon) & \
                         (table.lon_max > lon)
    rows = db(query).select(table.name,
                            table.level,
                            table.wkt,
                            )
    L1 = L2 = None
    for row in rows:
        shape = wkt_loads(row.wkt)
        ok = test.intersects(shape)
        if ok:
            if row.level == "L1":
                L1 = row.name
                if L2:
                    break
            elif row.level == "L2":
                L2 = row.name
                if L1:
                    break
    WKT = '"%s"' % wkt_dumps(geom)
    write("BT,%s,%s,%s,%s\n" % (L1, L2, name, WKT))

outputFile.close()