# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# + [markdown] toc="true"
# # Table of Contents
#  <p><div class="lev1 toc-item"><a href="#Overview" data-toc-modified-id="Overview-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Overview</a></div><div class="lev1 toc-item"><a href="#Load-LAU1-and-LAU2-codes" data-toc-modified-id="Load-LAU1-and-LAU2-codes-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Load LAU1 and LAU2 codes</a></div><div class="lev2 toc-item"><a href="#Create-the-hierarchy-of-the-LAU-codes" data-toc-modified-id="Create-the-hierarchy-of-the-LAU-codes-21"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>Create the hierarchy of the LAU codes</a></div><div class="lev2 toc-item"><a href="#Check-that-the-indexes/code-association-holds-at-each-level..." data-toc-modified-id="Check-that-the-indexes/code-association-holds-at-each-level...-22"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>Check that the indexes/code association holds at each level...</a></div><div class="lev2 toc-item"><a href="#Compute-the-population-for-each-shape-in-the-finest-level" data-toc-modified-id="Compute-the-population-for-each-shape-in-the-finest-level-23"><span class="toc-item-num">2.3&nbsp;&nbsp;</span>Compute the population for each shape in the finest level</a></div><div class="lev2 toc-item"><a href="#Compute-the-baricenter-of-each-LAU2" data-toc-modified-id="Compute-the-baricenter-of-each-LAU2-24"><span class="toc-item-num">2.4&nbsp;&nbsp;</span>Compute the baricenter of each LAU2</a></div><div class="lev2 toc-item"><a href="#Save-the-dataframes-and-mappings" data-toc-modified-id="Save-the-dataframes-and-mappings-25"><span class="toc-item-num">2.5&nbsp;&nbsp;</span>Save the dataframes and mappings</a></div>
# -

# %pylab inline
import pandas as pd
import os
import sys
import geopandas as pd
import overpy
import urllib2
import codecs

api = overpy.Overpass()

# +
searchByAreaName = False
countryOriginalName = "United-Kingdom"
searchByAreaID = True
areaID = "3600062149"

levelsToFetch = [6, 8, 10]
# -

for levelToFetch in levelsToFetch:
    if searchByAreaName:
        areaQuery = ('( area["{0}"="{1}"]; )->.a;'.format("name", countryOriginalName),
                     "type", "boundary", "admin_level", "%d" % levelToFetch)
    elif searchByAreaID:
        areaQuery = ("area(%s)->.a;" % areaID,
                     "type", "boundary", "admin_level", "%d" % levelToFetch)

    query = """
        [timeout:600];
        {0}
        (
              node["{1}"="{2}"]["{3}"="{4}"](area.a);
              way["{1}"="{2}"]["{3}"="{4}"](area.a);
              relation["{1}"="{2}"]["{3}"="{4}"](area.a);
        );
        out body;
        >;
        out skel qt;
        """.format(*areaQuery)
    
    query = query.replace("\n", "")
    query = query.replace("   ", "")
    query = query.replace("  ", "")
    
    print "Doing query: ", query
    
    result = urllib2.urlopen("https://overpass-api.de/api/interpreter", query)
    res = result.read()
    result.close()
    
    print "Query done with code: ", result.code, "\n\n"
    
    ofolder = os.path.join("resources/geoJsons", countryOriginalName)
    if not os.path.exists(ofolder): os.makedirs(ofolder)
    
    with codecs.open(os.path.join(ofolder, "dec_lvl%02d.osm" % (levelToFetch)), "w", "utf-8-sig") as f:
        f.write(res.decode("utf-8"))
    
    with open(os.path.join(ofolder, "enc_lvl%02d.osm" % (levelToFetch)), "w") as f:
        f.write(res)

target



# +
# Convert osm data to geojson and shapefile...
import glob
import geopandas as gpd
from shapely.geometry import MultiPolygon
from subprocess import Popen, PIPE

for osmdatafile in glob.glob(os.path.join(ofolder, "dec_*.osm")):
    target = osmdatafile.replace(".osm", ".json")
    with open(target, "w") as fout:
        args = "node --max_old_space_size=16384 /usr/local/bin/osmtogeojson -m %s" % osmdatafile
        args = args.split(" ")
        p = Popen(args, stdout=fout, stderr=PIPE, shell=False)
        out, err = p.communicate()
        print osmdatafile, err
    
    target_shp = target.replace(".json", ".shp")
    gdf = gpd.read_file(target)
    gdf = gdf[gdf.geometry.type.isin(["Polygon", "MultiPolygon"])]
    gdf.geometry = gdf.geometry.apply(lambda p: p if p.type == "MultiPolygon" else MultiPolygon([p]))
    gdf.to_file(target_shp, driver="ESRI Shapefile")
# -







StarcoderdataJupyterScriptsDedupFiltered