<jupyter_start><jupyter_text>Urban blocks  Imports & Parameters<jupyter_code>import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import pandas as pd
import rasterio
from tqdm import tqdm_notebook as tqdm
from shapely.geometry import shape
# Add local module to the path
src = os.path.abspath('../src')
if src not in sys.path:
    sys.path.append(src)
from metadata import City, CITIES
from landsat import Scene
from raster import euclidean_distance
%matplotlib inline
%load_ext autoreload
%autoreload 2
def pretty(city_id):
    """formatmat city label."""
    name = city_id.title()
    name = name.replace('_', '-')
    return name
# Plotting style
plt.style.use(['seaborn-paper', 'seaborn-whitegrid'])
plt.rc('font', family='serif', serif='cmr10', size=12)
plt.rc('axes', titlesize=12)
plt.rc('axes', labelsize=12)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('legend', fontsize=12)
plt.rc('figure', titlesize=12)
plt.rc('text', usetex=True)<jupyter_output><empty_output><jupyter_text>Data availability  Road network<jupyter_code>road_network = pd.Series(index=CITIES)

for city_name in CITIES:
    
    city = City(city_name)
    roads = city.roads
    area_of_interest = shape(city.aoi)
    road_network.loc[city_name] = round(np.sum(roads.length) * 1e-3)
road_network.sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Only three case studies have a road network that does not reach 1000km (Gao, Saint-Louis and Chimoio).  Urban blocks<jupyter_code>availability = pd.Series(index=CITIES)

for city_name in CITIES:
    
    city = City(city_name)
    availability.loc[city_name] = len(city.blocks)
availability.sort_values(ascending=False)
availability.index = availability.index.map(pretty)
availability.plot(kind='barh', color='C0')
plt.xlabel('Number of blocks (polygons)');<jupyter_output><empty_output><jupyter_text>Block size<jupyter_code>surfaces = np.arange(1, 50, 1)  # In hectares
spectral_distance = pd.DataFrame(index=CITIES, columns=surfaces)
count = pd.DataFrame(index=CITIES, columns=surfaces)

#progress = tqdm(total=len(CITIES) * len(surfaces))

for city_name in CITIES:
    
    city = City(city_name)
    scene = Scene(city.landsat_dir)
    blocks_sizes = city.blocks_raster
    
    bu = city.reference == 1
    sign_ref = scene.signature(bu)
    
    for max_ha in surfaces:
        
        blocks = (blocks_sizes > 0) & (blocks_sizes <= max_ha)
        count.at[(city_name, max_ha)] = np.count_nonzero(blocks)
        
        if np.count_nonzero(blocks) == 0:
            #progress.update(1)
            continue
        
        sign_blocks = scene.signature(blocks)
        spectral_distance.at[(city_name, max_ha)] = euclidean_distance(sign_blocks, sign_ref)
        #progress.update(1)
    
#progress.close()
spectral_distance
count
f, axes = plt.subplots(1, 2, figsize=(8, 4))

axes[0].plot(spectral_distance.mean(axis=0))
axes[0].set_xlim(0, 50)
axes[0].set_xlabel('Maximum block surface (hectares)')
axes[0].set_ylabel('Spectral distance to built-up')
axes[0].axvline(x=3, color='C3', linestyle=':', label='Optimal\nthreshold')
axes[0].legend(frameon=True)

for city in ['chimoio', 'gao', 'katsina', 'saint_louis', 'windhoek']:
    axes[1].plot(count.loc[city], label=pretty(city))

axes[1].set_xlim(0, 50)
axes[1].set_xlabel('Maximum block surface (hectares)')
axes[1].set_ylabel('Number of samples (pixels)')
axes[1].axvline(x=3, color='C3', linestyle=':', label='Optimal\nthreshold')
axes[1].legend(frameon=True, ncol=2, loc='upper left')

f.tight_layout()
f.savefig('../paper/figures/urban_blocks.pdf')
f.savefig('../paper/figures/urban_blocks.png');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14039806347498654, 'fraction_numerical': 0.016675632060247445, 'mean_word_length': 5.0967213114754095, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '532456', 'license': 'MIT', 'n_tokens_mistral': 1426, 'n_tokens_neox': 1335, 'n_words': 265, 'path': 'notebooks/05-Urban_Blocks.ipynb', 'repo_name': 'yannforget/builtup-classification-osm', 'size': '3718'}