<jupyter_start><jupyter_text>WITH OVERPASS<jupyter_code>import geojson
import requests
from shapely.geometry import box, Polygon, MultiPolygon, GeometryCollection, shape
from shapely import wkt
from shapely.ops import transform
import shapely
from functools import partial
import pyproj
import time
import json

n_processes = 10
def simplify(s, delta=0.05):
    
    while not check_length(s):
        s = s.simplify(delta, False)
        delta = delta + 0.05

    return s

def check_length(s, threshold=3000):
    
    return len(str(s)) < threshold

def get_area(s):
    s = shape(s)
    proj = partial(pyproj.transform, 
                   pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857'))

    return transform(proj, s).area / 1e6 # km

def threshold_func(g, value):
    
    return get_area(g) < value

def katana(geometry, threshold_func, threshold_value, count=0):
    """Split a Polygon into two parts across it's shortest dimension"""
    bounds = geometry.bounds
    width = bounds[2] - bounds[0]
    height = bounds[3] - bounds[1]
    if threshold_func(geometry, threshold_value) or count == 250:
        # either the polygon is smaller than the threshold, or the maximum
        # number of recursions has been reached
        return [geometry]
    if height >= width:
        # split left to right
        a = box(bounds[0], bounds[1], bounds[2], bounds[1]+height/2)
        b = box(bounds[0], bounds[1]+height/2, bounds[2], bounds[3])
    else:
        # split top to bottom
        a = box(bounds[0], bounds[1], bounds[0]+width/2, bounds[3])
        b = box(bounds[0]+width/2, bounds[1], bounds[2], bounds[3])
    result = []
    for d in (a, b,):
        c = geometry.intersection(d)
        if not isinstance(c, GeometryCollection):
            c = [c]
        for e in c:
            if isinstance(e, (Polygon, MultiPolygon)):
                result.extend(katana(e, threshold_func, threshold_value, count+1))
    if count > 0:
        return result
    # convert multipart into singlepart
    final_result = []
    for g in result:
        if isinstance(g, MultiPolygon):
            final_result.extend(g)
        else:
            final_result.append(g)
    return final_result




def to_geojson(x, full=False):
    g = geojson.Feature(geometry=x, properties={}).geometry
    if full:
        return g
    if g['type'] == 'MultiPolygon':
        return g['coordinates'][0][0]
    else:
        return g['coordinates'][0]

def swipe(x):
    return [[c[1], c[0]] for c in x]
    
def flatten(l):
    return [str(round(item, 4)) for sublist in l for item in sublist]


def to_overpass_coords(x):
    
    coords = to_geojson(x)
    coords = swipe(coords)
    coords = flatten(coords)
    coords = ' '.join(coords)
    return coords

def overpass_request(s):
    
    print(s['region_slug'])
    print(s['partition'])
    
    overpass_url = "http://overpass-api.de/api/interpreter"
    overpass_query = """
    [out:json];
    way["highway"]["highway"!~"^(construction|cycleway|footway|path|proposed|service|track|pedestrian|living_street|platform|steps)$"]["highway"](poly:"%s");
    for (t["highway"])
    {
       make stat highway=_.val,
           count=count(ways),length=sum(length());
       out;
    }""" % s['overpass_coords']
    
    for i in range(3):
        try:
            response = requests.get(overpass_url, 
                                    params={'data': overpass_query})

            json.dump(response.json(), open(RAW_PATH / 'osm_region_length' / (s['region_slug'] + '_' + str(s['partition'])), 'w'))
            break
        except Exception as e:
            if i == 2:
                print(e)
                print(response.text)
            print('retrying')
            continue
#     return response
url = 'https://docs.google.com/spreadsheets/d/197fccIfjwlsT-oeisHHLn2SuOvbtJjC0inoNuxfqrwI/gviz/tq?tqx=out:csv&sheet=metadata'

metadata = pd.read_csv(url)
metadata = metadata[[c for c in metadata.columns if 'Unnamed' not in c]]

metadata['region_shapefile_wkt'] = metadata['region_shapefile_wkt'].apply(lambda x: wkt.loads(x.replace('"', '')))

# From MultiPolygon to Polygon
metadata['region_shapefile_wkt'] = metadata['region_shapefile_wkt']\
    .apply(lambda x: max(x, key=lambda a: a.area) if isinstance(x, shapely.geometry.multipolygon.MultiPolygon) else x)

# metadata['region_shapefile_wkt'].apply(simplify)
metadata
threshold_value = 1000000
def partition_geometry(_df, threshold_value):
    
    geos = katana(_df['region_shapefile_wkt'], threshold_func, threshold_value)

    final = []
    for i, g in enumerate(geos):
        _df['partition'] = i
        _df['shapefile_partition'] = g
        final.append(deepcopy(_df.to_frame().T))

    return pd.concat(final)
metadata_part = pd.DataFrame()
for i, row in metadata.iterrows():
    print(row['region_slug'])
    metadata_part = pd.concat([metadata_part , partition_geometry(row, threshold_value)])
metadata_part = metadata_part[metadata_part['region_slug'].isin(missing)]
metadata_part['shapefile_partition'] = metadata_part['shapefile_partition'].apply(simplify)
metadata_part['overpass_coords'] = metadata_part['shapefile_partition'].apply(to_overpass_coords)
rows = [row for i, row in metadata_part.iterrows()]
with Pool(1) as p:
    p.map(overpass_request, rows)
regions = set(metadata['region_slug'])
def get_length(j):
    a = pd.DataFrame([r['tags'] for r in j['elements']])
    if len(a):
        return round(a['length'].astype(float).sum(), 2)

paths = (RAW_PATH / 'osm_region_length')
names = [p.name for p in paths.glob('*')]
length = [get_length(json.load(open(s, 'r'))) for s in paths.glob('*')]
length = pd.DataFrame(list(zip(names, length)), columns=['region_slug', 'total_length'])
length['region_slug'] = length['region_slug'].apply(lambda x: '_'.join(x.split('_')[:-1]))
length = length.groupby('region_slug').sum().sort_values(by='total_length')
missing = regions.difference(set(length.index))
missing
!rm -r /home/joaom/projects/waze_coronavirus/data/raw/osm_region_length/.ipynb_checkpoints
# No partition
length.sum()
# threshold_value = 1000000
length.sum()
# threshold_value = 100000
length.sum()
a = set(metadata['region_slug'].unique())
b = set(length.dropna()['region_slug'].unique())
length.dropna().sort_values(by='total_length')
geometry = metadata[metadata['region_slug'] == 'country_brazil']['region_shapefile_wkt'].iloc[0]
daily = pd.read_csv('https://bit.ly/idb-traffic-daily')
daily
import folium
m = folium.Map(zoom_start=1, tiles='cartodbpositron')
for i in a:
    m.add_children(folium.GeoJson(i))
m
paths = (RAW_PATH / 'osm_region_length')
names = [p.name for p in paths.glob('country_brazil_*')]
pd.DataFrame([], columns=['region_shapefile_wkt'])
import osm_road_length
a = metadata['region_shapefile_wkt'].iloc[0]
import pyarrow
import awswrangler as wr
 wr.s3.to_parquet(
    df=pd.DataFrame({
        'col': [1, 2, 3],
        'col2': ['A', 'A', 'B']
    }),
    path='s3://bucket/prefix',
    dataset=True,
    partition_cols=['col2']
)
import geobr
year = 2010
types = {
    'States': geobr.read_state(code_state='all', year=year),
    'Meso Regions': geobr.read_meso_region(code_meso="all", year=year),
    'Micro Regions': geobr.read_micro_region(code_micro='all', year=2017),
    'Municipalities': geobr.read_municipality(code_muni='all', year=year)
}
brasilio = pd.read_csv('https://data.brasil.io/dataset/covid19/caso.csv.gz')
last = brasilio[brasilio['is_last']]
last.head(1)
df.merge(last, left_on='code_muni', right_on='city_ibge_code')
fig, axis = plt.subplots(nrows=1, ncols=1,
                      figsize=(15, 15))

df.merge(last, left_on='code_muni', right_on='city_ibge_code')\
.plot(facecolor='#2D3E50', edgecolor='#FEBF57', ax=axis)
        
axis.set_title("Municípios com casos de Covid19", fontsize=20)
axis.axis('off')

fig.tight_layout()
nrows = 2
ncols = 2
fontsize = 20

fig, ax = plt.subplots(nrows=nrows, ncols=ncols,
                      figsize=(15, 15), dpi=300)

types_iter = iter(types.items())

for r in range(nrows):
    for c in range(ncols):
        
        axis = ax[r][c]
        name, df = next(types_iter)
        
        df.plot(facecolor='#2D3E50', edgecolor='#FEBF57', ax=axis)
        
        axis.set_title(name, fontsize=fontsize)
        axis.axis('off')

# General
fig.tight_layout()
fig.savefig(f'imgs/1.geobr_logo.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14748244256636114, 'fraction_numerical': 0.017021783121057017, 'mean_word_length': 3.639425731639978, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '255829', 'license': 'MIT', 'n_tokens_mistral': 3037, 'n_tokens_neox': 2905, 'n_words': 588, 'path': 'notebooks/Add length to metatada.ipynb', 'repo_name': 'EL-BID/Covid-19-Traffic-Impact-Dashboard', 'size': '8401'}