<jupyter_start><jupyter_text>Create multi-label chips for multi-output regression sampled from the whole GBCreate chips and store proportions of signature types within as an input of multi-output regression problem<jupyter_code>import geopandas
import tobler
import pyogrio
import pygeos
import numpy
import pandas
import dask_geopandas
import rasterio
from scipy.sparse import coo_matrix

import dask
import dask.bag
from dask.distributed import Client, LocalCluster
%%time
df = geopandas.read_parquet("/home/jovyan/work/chips_gb/chip_bounds_32/").reset_index(drop=True)
df
signatures = pyogrio.read_dataframe(        
        '/home/jovyan/work/urbangrammar_samba/spatial_signatures/'
        'signatures/'
        'signatures_combined_levels_simplified.gpkg'
)
bds = df.total_bounds
signatures = signatures.cx[bds[0]:bds[2], bds[1]:bds[3]]
%%time

ids_src, ids_tgt = df.sindex.query_bulk(signatures.geometry, predicate="intersects")
od_matrix = pandas.DataFrame(dict(ids_src=ids_src, ids_tgt=ids_tgt))
%%time 
sjoined = df.set_geometry(df.centroid).sjoin(signatures[["signature_type", "geometry"]], how="left", predicate="within")
od_matrix["signature_type"] = sjoined.signature_type.iloc[ids_tgt].values
od_matrix.signature_type.value_counts().plot.bar()
od_matrix.signature_type.value_counts()
cap = 20000
keep = []
counts = {k:0 for k in od_matrix.signature_type.unique()}
for i, t in od_matrix.signature_type.sample(len(od_matrix)).iteritems():
    counts[t] += 1
    if counts[t] < cap:
        keep.append(i)
limited = od_matrix.iloc[keep]
limited.signature_type.value_counts().plot.bar()
limited.shape
%%time
chips = df.geometry.values[limited.ids_tgt.values].data
sig = signatures.geometry.values[limited.ids_src.values].data

# this could be parallelised
r = [
    pygeos.clip_by_rect(
        g, *chip
    ) for g, chip in zip(sig, pygeos.bounds(chips))
]

areas = pygeos.area(r)

table = coo_matrix(
    (
        areas,
        (limited.ids_src.values, limited.ids_tgt.values),
    ),
    shape=(signatures.shape[0], df.shape[0]),
    dtype=numpy.float32,
)

table = table.tocsr()
%%time
signatures["signature_type"] = pandas.Categorical(signatures["signature_type"])
%%time
categorical = {}
for value in unique:
    mask = signatures["signature_type"] == value
    categorical[value] = numpy.asarray(
        table[mask].sum(axis=0)
    )[0]
    print(value)
%%time
categorical = pandas.DataFrame(categorical)
%time areas = df.area.values
%time categorical = categorical.div(areas, axis="rows")
(categorical.sum(axis=1) == 0).sum()
geom = df[df.geometry.name].reset_index(drop=True)
ests = geopandas.GeoDataFrame(categorical, geometry=geom, crs=signatures.crs)
limited
ests = ests.iloc[limited.ids_tgt.values]
ests.drop(columns="geometry").sum().plot.bar()
ests.to_parquet("/home/jovyan/work/chips_gb/chip_proportions_32.pq")
ests.centroid.plot(figsize=(20, 20), markersize=.1)<jupyter_output><empty_output><jupyter_text>Create chips<jupyter_code>specs = {
    'chip_size': 32,
    'bands': [1, 2, 3], #RGB
    'mosaic_p': (
        '/home/jovyan/work/urbangrammar_samba/'
        'ghs_composite_s2/GHS-composite-S2.vrt'
    ),
}
bounds = geopandas.read_parquet("/home/jovyan/work/chips_gb/chip_proportions_32.pq")
centroid = bounds.centroid
bounds['X'] = centroid.x.astype(int)
bounds['Y'] = centroid.y.astype(int)
client = Client(
    LocalCluster(n_workers=16, threads_per_worker=1)
)
client
import numpy as np

def bag_of_chips(chip_bbs, specs, npartitions):
    '''
    Load imagery for `chip_bbs` using a Dask bag
    ...
    
    Arguments
    ---------
    chip_bbs : GeoDataFrame
               Geo-table with bounding boxes of the chips to load
    specs : dict
            Metadata dict, including, at least:
            - `bands`: band index of each band of interest
            - `chip_size`: size of each chip size expressed in pixels
            - `mosaic_p`: path to the mosaic/file of imagery
    npartitions : int
                  No. of partitions to split `chip_bbs` before sending to
                  Dask for distributed computation
    Returns
    -------
    chips : ndarray
            Numpy tensor of (N, chip_size, chip_size, n_bands) dimension 
            with imagery data   
    '''
    # Split chip_bbs
    thr = np.linspace(0, chip_bbs.shape[0], npartitions+1, dtype=int)
    chunks = [
        (chip_bbs.iloc[thr[i]:thr[i+1], :], specs) for i in range(len(thr)-1)
    ]
    # Set up the bag
    bag = dask.bag.from_sequence(
        chunks, npartitions=npartitions
    ).map(chip_loader)
    # Compute
    chips = np.concatenate(bag.compute())
    return chips


def chip_loader(pars):
    '''
    Load imagery for `chip_bbs`
    ...
    
    Arguments (wrapped in `pars`)
    -----------------------------
    chip_bbs : GeoDataFrame
               Geo-table with bounding boxes of the chips to load
    specs : dict
            Metadata dict, including, at least:
            - `bands`: band index of each band of interest
            - `chip_size`: size of each chip size expressed in pixels
            - `mosaic_p`: path to the mosaic/file of imagery
    Returns
    -------
    chips : ndarray
            Numpy tensor of (N, chip_size, chip_size, n_bands) dimension 
            with imagery data
    '''
    chip_bbs, specs = pars
    b = len(specs['bands'])
    s = specs['chip_size']
    chips = np.zeros((chip_bbs.shape[0], b, s, s))
    with rasterio.open(specs['mosaic_p']) as src:
        for i, tup in enumerate(chip_bbs.itertuples()):
            img, transform = rasterio.mask.mask(
                src, [tup.geometry], crop=True, all_touched=True
            )
            img = img[:b, :s, :s]
            for ban, (l_min, l_max) in enumerate([(350, 1600), (500, 1600), (600, 1800)]):
                img[ban][img[ban] > l_max] = l_max
                img[ban][img[ban] < l_min] = l_min
                a_std = (img[ban] - l_min) / (l_max - l_min)
                img[ban] = a_std * 255
            chips[i, :, :, :] = img
    chips = np.moveaxis(chips, 1, -1)
    return chips.astype(rasterio.uint8)
chips = bag_of_chips(bounds, specs, 16)
numpy.save('../../chips_gb/multilabel_chip_32.npy', chips)<jupyter_output>distributed.worker - WARNING - Unmanaged memory use is high. This may indicate a memory leak or the memory may not be released to the OS; see https://distributed.dask.org/en/latest/worker.html#memtrim for more information. -- Unmanaged memory: 6.37 GiB -- Worker memory limit: 7.85 GiB
distributed.worker - WARNING - Unmanaged memory use is high. This may indicate a memory leak or the memory may not be released to the OS; see https://distributed.dask.org/en/latest/worker.html#memtrim for more information. -- Unmanaged memory: 5.57 GiB -- Worker memory limit: 7.85 GiB
distributed.worker - WARNING - Unmanaged memory use is high. This may indicate a memory leak or the memory may not be released to the OS; see https://distributed.dask.org/en/latest/worker.html#memtrim for more information. -- Unmanaged memory: 6.28 GiB -- Worker memory limit: 7.85 GiB
distributed.worker - WARNING - Unmanaged memory use is high. This may indicate a memory leak or the memory may not be released to the OS; see [...]
StarcoderdataJupyterStructuredCleanDedup