Capturing deprived areas using unsupervised machine learning and open data: a case study in São Paulo, Brazil

ABSTRACT

Managing the rapid growth of deprived areas (commonly known as slums, informal settlements, etc.) in cities of Low- to Middle-Income Countries (LMICs) demands detailed and consistent information that is often unavailable. Recent Earth Observation (EO) mapping approaches with supervised classification models overlook the diversity of deprived areas and require resource-intensive training sets. In this study, we analyse the potential of unsupervised machine learning (ML) models to capture intra-urban diversity of deprived areas in São Paulo, using solely open geodata. We provide a workflow of characterising deprivation at a city scale with a disaggregated approach, offering scalability and transferability potential. First, we extract a pool of spatial features from open geospatial datasets to characterise the morphological and environmental conditions of the study area. After input preparation, we train and optimise a k-means model, including a coupled feature importance tool. Four cluster types emerged with different deprivation aspects such as higher and lower accessibility to services and infrastructure, sparser and denser occupation; regular and complex morphology; flat and steep terrain. This alternative methodology to capture diversity of deprived areas with open EO-based features can inform locally targeted, thus more efficient, urban policies and interventions.

KEYWORDS:

• Deprivation
• slums
• remote sensing
• unsupervised machine learning
• Clustering, low-to-middle-income countries

Introduction

One out of eight people in the world lives in slum-like settlements, marked by poor housing and service provision (UN-Habitat, 2016). Most of them are located in cities of Low-to-Middle-Income Countries (LMIC), where high urban growth rate and weak planning capacities coincide (Bastos da Cunha et al., 2015). This leads to high housing demands that cannot be accommodated by the local governments, thus, to fast development of deprived areas (Kohli et al., 2012).

The traditional concept of slums has been used widely in literature (UN-Habitat, 2003), and despite its irrefutable significance, it is limited to identify the deprived at the household level and overlooks the areal characteristics of slums (Lucci et al., 2018; Olthuis et al., 2015). This definition leads to entire areas being classified as slums, disregarding the area-level nuances and multiple manifestations (Thomson D R et al., 2020; Patel et al., 2014). Therefore, this study adopts the term “deprived area” as neighbourhoods facing deprivation at social, environmental and ecological dimensions that differ across and within cities, and are highly dependent on their local contexts (Abascal et al., 2021).

With the current international interest in slum improvement, there are multiple governmental efforts to upgrade these areas, requiring consistent and up to date spatial information about their living conditions to support effective pro-poor plans and interventions (Kohli et al., 2013; Thomson et al., 2020). Yet, acquiring such information is a challenge. The most traditional method to acquire data on deprived areas is census surveys. The fieldwork produces information at the household level typically aggregated at administrative units, which are likely to ignore the morphological local heterogeneities and stigmatise neighbourhoods by assigning average values to the entire unit (Martínez et al., 2016; Taubenböck et al., 2018).

Besides, the large temporal gaps and the resource-consuming process of household data collection (e.g. Census) are a huge challenge to LMICs (Owen & Wong, 2013). To overcome such issues, the increasing availability of Earth Observation (EO) data allows disaggregated information with high temporal resolution. The recent remote sensing (RS) studies extract spatial, spectral and textural features from satellite imagery and translate it into information on the living conditions in deprived areas (Leonita et al., 2018). However, even with the rising availability of low-cost drone imagery, the acquisition of Very-High-Resolution (VHR) multispectral data is still cost prohibitive, restraining the EO-based approaches to generate image-based information to small urban districts (Gram-Hansen et al., 2019; Merodio et al., 2021; Small, 2014).

In this context, the potential of open geodata sources – such as the United States Geological Survey (USGS), the Copernicus Program and Open Street Maps (OSM) platforms – to provide a rich repository for characterizing different urban areas have been systematically explored in recent years (Bhuyan et al., 2022; Geiß et al., 2017). Within the deprivation studies, however, only a few studies have incorporated such datasets due to the harnessing and integration challenges (Ajami et al., 2019; Mahabir et al., 2020). Aditionally, the quality and completeness of opengeodata, especially derived from volunteered mapping initiatives, should be inspected for uncertainties. Cartographic and thematic innacuracies, non-systematic updates and level of completeness across different areas of the globe are some of the limitations to consider when building an approach on such datasets (Eckle et al., 2016; Yeboah et al., 2021).

In addition to data constraints to acquire information on deprived areas, literature shows methodological shortcomings. EO-based studies present different purposes and modelling approaches. Most research is focused on mapping the settlements, i.e. delineation purposes, classifying deprived and non-deprived areas based on their physical characteristics extracted from VHR imagery (Gevaert et al., 2019; Kohli et al., 2016). There are different methods to identify deprived areas – visual interpretation, pixel-based image classification and object-based image analysis (Dos Santos et al., 2022) – but recent Machine Learning (ML) and Deep Learning (DL) methods have demonstrated high performance (Ajami et al., 2019; Duque et al., 2017; Luo et al., 2022).

Supervised ML models are developed to classify deprived settlements in (semi-)automatic ways, combining spectral and textural information (Kuffer et al., 2016); or using ML-OBIA techniques (Duque et al., 2017); or pixel-based classifiers such as in Mahabir et al. (2020). Despite the undeniable importance of the detection models, they overlook the diverse nature of deprived settlements, obstructing the development of more contextualised, hence, effective policies (Kuffer et al., 2020). Moreover, the few studies that do work on capturing the different characteristics of deprived areas (Ajami et al., 2019; Luo et al., 2022) rely on the acquisition of high-quality training and reference data for the supervised models, making adoption in LMICs difficult (Jochem et al., 2020). This task is not only cost-prohibitive due to high VHR costs, but also resource-intensive, which hampers the scalability of these models (Gevaert et al., 2019). Especially in the context of COVID-19 pandemic, there is a pressing demand for data-driven methods that can provide an alternative to field work (Brito et al., 2020). Modelling with unsupervised ML techniques could be posed as an alternative, but till this day, no gridded unsupervised ML approach to capture local urban deprivation was documented in the literature. Even though non-gridded approaches are more common in research, they are much more computationally demanding and require the availability of homogeneous units (e.g. land parcels or streets blocks) that are often inaccessible in LMICs (Kuffer et al., 2020).

Motivated by these challenges, the present study proposes a workflow, exploring the potential unsupervised ML models to capture intra-urban diversity of deprived areas at city scale, using solely open data sources. The city of São Paulo, Brazil, serves as a case study.

The structure of the remainder of this paper is as follows. First part, we describe the research methods and materials, covering the case study area, the data collection, the data processing, the unsupervised machine learning model and the validation process. Then, we present the results, followed by the discussion. Finally, we summarise the results and provide recommendations for further research.

Materials and methods

The general workflow to capture deprived areas using a gridded unsupervised approach is presented in Figure 1, based on four main steps as follows: (i) we select a pool of indicators and collect open data sources, considering a specific set of requirements and the physical aspects of deprivation of the study area, (ii) define a sensible unit of analysis, extract spatial features and integrate them, (iii) employ unsupervised ML model and feature importance techniques and (iv) assess the model results, including multiple validation procedures. After the subsection about the case study, the following four subsections describe each one of the four steps presented in the methodology flowchart.

Figure 1. Overview of the proposed methodology.

Case study: São Paulo

São Paulo, capital of the State of São Paulo, is the most populated city in Brazil (12.3 million) and covers 1.521 km² according to the latest census (IBGE, 2010). We selected the city to develop the methodology, considering three main reasons: (1) The city’s complex layout and unprecedented growth rate of deprived settlements towards the peripheral areas. The population living in the informal areas grew four times the formal population, which poses an urgent and pressing challenge to policy makers (Pasternak & D’Ottaviano, 2016). (2) In 2019, the Agency of Urban Development and the Federal University of ABC Region in São Paulo developed a mapping and classification study for the Metropolitan Area of the State of São Paulo using Logistic Regression (da Fonseca Feitosa et al., 2021). The study was applied in different cities, but till the present date, no research has been released in the capital city. (3) The availability of a recent national data set on the spatial extent of deprived areas. The existence of an open-source GIS layer delineating the deprived areas is a prime requirement of the present workflow. The Brazilian Institute of Geography and Statistics (IBGE) released the new census layer in 2019, mapping deprived areas as Aglomerados Subnormais - abbreviated as AGSN and directly translated as Substandard Settlements (Figure 2). The concept, used since 1991, describes settlements with at least 51 houses in a substandard urban setting, deprived of legal ownership, regular arrangement and service provision (IBGE, 2010).

Figure 2. Study area Sao Paulo with substandard settlements in 2019. Source: (IBGE, 2019).

The previous AGSN layer, from the 2010 census, has omission and topological issues. Our visual inspection proved similar results to Ferreira and Feitosa (2020) that investigated such inconsistencies on the new layer and found better spatial coverage and fixed topological errors.

Step 1 - data collection

EO-based studies have shown that the physical characteristics of deprived areas can be derived from RS imagery that can be acquired worldwide (Ghaffarian et al., 2018; Leonita et al., 2018). However, in view of their underlying heterogeneity, some local adjustments are necessary in order to choose the spatial features that can translate most of the socioeconomic information of the study area (Olthuis et al., 2015; Wurm et al., 2017). We structure the data collection process by summarising the observable features that can indicate the spatial variations of deprived areas on the ground (Figure 3). Figure 3 shows the built-up morphology domain that encompasses the differences in building geometry (size, shape), density and layout (texture) characteristics, as well as accessibility to services and infrastructure; and the land morphology domain including topography, presence of green areas and hydrography.

Figure 3. Diagram of the domains of the diversity of deprived areas.

Besides the conceptual structure, the data collection process was guided by the following set of requirements, relying on the work of Mahabir et al. (2020): all features must be spatial, quantitative, available for the entire study area, recurrently updated, provided from open sources and manifest the local nuances of the study area.

Given the concepts and requirements outlined above, we compiled a list of 32 spatial features identified in literature to capture deprivation to build an open geodatabase for São Paulo (Table 1). In addition, we collected two auxiliary datasets to support the validation of the results: Census data (2010), acquired from the IBGE web portal and aggregated at the census tract level – the 10-year gap remains only due to the suspension of the 2020 census activities of 2020 due to pandemic – and the municipal land use layer Considering the replicability potential of the model and the unprecedented growth rate of deprived areas, there is a risk of erroneous assessment of the mapping results. For this reason, we focused this statistical validation on analysing the separability capacity of the clusters between socioeconomic aspects of deprivation.

Table 1. List of EO-based features found in literature.

Dimension	Feature	Abbreviation	Rationale	Data source	Specification
Texture	GLCM-texture measures	MEAN VARIANCE HOMOGENEITY CONTRAST DISSIMILARITY ENTROPY 2^ND_MOMENT	Asymmetry of size and shape of rooftops can distinguish deprived areas (Duque et al., 2017; Owen & Wong, 2013)	Landsat-8 imagery (BOA reflectance)	Band: Panchromatic/Resolution: 15 m/Acquired in 2020
Density	Population count	POP_COUNT	Building layout can vary depending on the stage of development of the deprived area and can be a proxy of adequate living space (Kohli et al., 2016; Sliuzas et al., 2010)	WorldPop Constrained Population Count layer	Resolution: 3 arc-second/UN adjusted/2020
	Built up density	PATCH_DENS		World Settlement Footprint (WSF) layer	Resolution: 10 m/2015
Geometry	Built up shape index	SHAPE_MN	Shape and layout of buildings can indicate morphological patterns of deprived areas (Leonita et al., 2018; Taubenböck et al., 2018)	World Settlement Footprint (WSF) layer	Resolution: 10 m/2015
	Built up mean area	AREA_MN
Spectral response	Raw spectral bands 2–7 and 11–12	BAND “number”	Raw image bands may be sensitive to morphological variations (Mahabir et al., 2020)	Sentinel-2A imagery	Bands: 2, 3,4,5,6,7,8, 11,12/Resolution: 10 m/Acquired in 2020
Infrastructure	Distance to major roads	DIST_ROADS	Deprived areas have lower accessibility (Lilford et al., 2019)	OpenStreetMap (OSM) layers	Shapefile/2020
	Unpaved road density	UNPAVED	Proximity to low-quality roads can indicate low economic status (Nobrega et al., 2008)
	Dead end nodes density	DEAD_END	Slums often have low internal road access (Taubenböck et al., 2018)
	VIIRS night-time light	VIIRS	Energy access can measure poverty rates (Van Dijk et al., 2019)	WorldPop resampled VIIRS night-time lights layer	Resampled from 15-arc second/2016
Services	Distance to schools	DIST_EDU	Dwellers of deprived areas often have lower access to service provision and job opportunities (Kohli et al., 2013; Kuffer et al., 2016)	Humanitarian OpenStreetMap Team (HOTOSM) layers	Shapefile/2020
	Distance to health facilities	DIST_HEALTH
	Distance to finance facilities	DIST_FINAN
Greenery	NDVI	NDVI	Vegetation can indicate different stage of development and location of slums (Kuffer et al., 2013)	Sentinel-2A imagery	Bands 4 and 8/Resolution: 10 m/Acquired in 2020
	Proportion of green areas	DIST_GREEN	Sparse deprived settlements are often close to green amenities.(Ebert et al., 2009)	World Database on Protected Areas (WDPA) layer	Shapefile/2020
Topography	Slope	SLOPE	High exposure to hazard is common to deprived areas; they are often located in low-lying or steep slope regions (Olthuis et al., 2015)	ASTER Global Digital Elevation Model (GDEM)	Resolution: 30 m/2019
	DEM	DEM
Water	Distance to nearest waterbody	DIST_WATER	Deprived areas are often concentrated near water banks and/or more exposed to flood risk (Ghaffarian et al., 2018; Kohli et al., 2016)	OpenStreetMap (OSM) layers	Shapefile/2020
	River network density	DENS_WATER

Step 2 - data processing

After data collection, we conduct a series of processing steps to extract the spatial features described in Table 1 and integrate them as the model input.

Spatial unit of analysis

Considering the benefits of providing gridded outputs with fuzzy continuous boundaries, we followed a pixel-based approach, and hence required the definition of a sensibly sized grid. This is a very important task, but rarely documented in literature. The few studies that do explain their decision process adopt a trial-and-error approach mostly considering: (1) the spatial resolution of the input datasets – here ranging from 10 to 100 m – and the morphology of the settlements (Taubenböck & Kraff, 2014; Wang et al., 2019); (2) the homogeneity of the sampled pixels (Owen & Wong, 2013). Thus, after inspecting the 2019 AGSN layer and correcting clear digitation errors, e.g. settlement polygons extrapolating into water bodies, we analysed the area of the 1,575 polygons through descriptive statistics, which indicate high variability in size (Table 2). Based on previous studies on mapping and characterisation of deprived areas, we created four regular grids of 10 × 10, 20 × 20, 50 × 50 and 100 × 100 m² cells and overlaid them with the AGSN layer to choose the one that best depicts the urban structure of the deprived areas in São Paulo (Duque et al., 2017; Kit et al., 2012; Mahabir et al., 2020).

Table 2. Area descriptive statistics of AGSN polygons.

Max	869,888 m^2
Min	745 m^2
Sum	61,724,714 m^2
Std. Deviation	70,890 m^2
Mean	39,190 m^2

Next, we excluded non-homogeneous pixels, i.e. where deprived and non-deprived areas exist, mainly found on the boundaries of the polygons. We performed a macro data-driven approach to establish the sampling threshold, indicating the minimal percentage of each pixel to be intersected by the AGSN layer, thus, the percentage of homogeneity of the deprived pixels. We assess three thresholds based on quartiles, comparing their total grid area with the total area of the AGSN polygons.

Feature extraction process

After the selection of the spatial unit of analysis, we derived the 32 GIS-and RS-based input features listed in Table 1. Figure 4 illustrates the workflow of seven main operations, some only encompassing standard processing steps – clip to extent, project, resample to chosen resolution – other followed by individualised extraction processes. All the operations shown in Figure 4 are performed in an open GIS platform with inbuilt functions, except from the fourth and fifth, that require specific processing steps to define the moving window size from which the features are calculated.

Figure 4. Processing steps of feature extraction structured in seven main operations.

Operation 1 extracted population count and night-time lights from WorldPop datasets. Operation 2 derived the Digital Elevation Model (DEM) and respective slope values from the Aster GDEM layer. Operation 3 included processing of spectral features derived from the mean values of Sentinel 2A imagery and the calculation of the Normalised Difference Vegetation Index (NDVI).

Operations 4 and 5 were based on the work of Hall-Beyer (2017) to compute a Principal Component Analysis (PCA) using Python packages¹ to support the window size decision, aka kernel size, by interpreting three main output metrics: correlation matrix, factor loadings and communality table. In operation 4, we extract the Gray-Level Co-occurrence Matrix (GLCM) from the Landsat 8 panchromatic band in R Studio to derive seven texture features – Mean, Variance, Homogeneity, Contrast, Dissimilarity, Entropy and Second Moment – from five different window sizes (3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11), which were further assessed by PCA metrics.

Operation 5 extracted the Landscape Metrics (LSM) from the World Settlement Footprint (WSF) layer. Relying on the study of Frazier & Kedron (2017), we calculated five metrics – Patch Area Mean, Shape Index, Fractal Dimension Index, Patch Density and Aggregation Index – using Fragstats with four window sizes (3 × 3, 5 × 5, 7 × 7 and 9 × 9), which are also evaluated with a PCA. The flowchart in Figure 5 illustrates the PCA experiments above in more detail. First, we assessed the metrics using all window sizes as input. If we could not derive a conclusion of which one maintains the highest PC loading scores, we ran the model with the features of one window size at a time.

Figure 5. Workflow of PCA iterations to choose the best moving window for hand-crafted features.

Operations 6 and 7 involve the extraction of the GIS-based features. After comparison, data from the Humanitarian OpenStreetMap (HOT) is used to calculate Euclidean distance to determine accessibility to services, while data from OSM is used to compute accessibility to infrastructure. The SQL queries expression to extract each attribute is described in Appendix 1. The WDPA layer is combined with OSM to measure the distance to green areas. Following the work of Mahabir et al. (2020), Kernel Density Estimation (KDE) is used to calculate the density-related features with a 1000 m search radius, which maintain the spatial structure of the features explicit even after the smoothing effect (Leonita et al., 2018).

Data Integration

After feature extraction, we integrated and prepared them as input for the unsupervised ML model. The chosen and sampled grid was converted to a point class layer, where the pixel centroids of each spatial feature was extracted and recorded. Given the size of the study area (149,947 cells, hence, points) and the processing capacity of GIS platforms, this process was performed in seven batches to encompass the entire extent.

Then, we merged the resulting attributes tables into a single georeferenced tabular file, export it to R Studio² and run Exploratory Data Analysis. Due to the differences in ranges among features values, we standardised the feature values (Z-Score) and analysed the descriptive statistics to detect outliers and avoid biases, check normal distribution and handle null values of population count as zero (Cai et al., 2010). Finally, we also compute the Pearson Correlation coefficient to investigate possible multicollinearity between features.

Step 3 - unsupervised machine learning model

After input preparation, we chose the best algorithm fit for the unsupervised ML task. Considering the instructions of Han et al. (2012), we selected the K-Means algorithm, because (1) it handles continuous numerical data such as the spatial features here used; (2) it performs well for high volume dataset with small processing time for regular machine capacities.³ Besides, when compared to other suitable algorithms, such as HDBSCAN, K-Means has simple parameters requirements, higher stability and compatibility with open programming languages (Campello et al., 2013; Madhulatha, 2012; McInnes et al., 2016). Considering the usability of the model in LMIC, simplicity and availability are fundamental technical criteria to the choice of the modelling approach in this study.

K-Means is one of the most popular unsupervised algorithms, falling into the category of portioning distance-based methods. The model initialises by portioning the data points based on “k” (number of) clusters, then it assigns the points to the closest cluster by calculating the mean-square distance (MSD) and comparing with the total within-cluster sum of square (WCSS), recomputing new cluster centroids accordingly and iteratively, until the clusters assignment no longer change. The WCSS, aka inertia or distortion score, sums the distance from each data point within a cluster to the cluster centroid (Chiang & Mirkin, 2010), hence, when the lowest WCSS score is constantly achieved, the highest internal coherence is achieved per cluster, so the model stops.

We ran the algorithm using Python⁴ and used the elbow method to choose the optimal k value. A plot of WCCS as a function of the number of clusters is used to determine the value of k when WCSS decreases exponentially and stabilises (Pedregosa et al., 2011). This breakpoint (“elbow”) suggests the optimal number of clusters. If more than one k value is indicated, it is suggested to train the models and assess the results (Han et al., 2012). We trained successive models, increasing the number of clusters from 3 to 10, computing and plotting all the total within-cluster sum of square (WCSS) results. The assessment considered the flattening of the elbow function with plots mentioned above and the computational costs (the more clusters, the longer the machine takes to derive them).

Optimisation techniques

In addition to EDA applied on the model input, literature highlights the usefulness of dimensionality reduction techniques for model optimisation (Martino et al., 2017). Clustering models like K-Means can be sensitive to the input, prone to overfitting and biassed results (Roy et al., 2020). There are several supervised feature selection algorithms, but little unsupervised algorithms are documented (Alelyani et al., n.d. Li et al., 2017). We included a feature importance tool in R, aka “FeatureImpCluster” (Pfaffel, 2021), that is documented, openly available in the Comprehensive R Archive Network (CRAN) and has not yet been applied to deprivation studies.

The tool calculates a permutation classification rate, i.e. the number of wrong cluster assignments divided by the number of iterations. When the algorithm changes a feature value and the model error elevates, the resulting permutation score is higher, indicating a higher relevance of the feature (Thrun et al., 2021). This tool works coupled with the k-means algorithm, computing the feature scores with simple statistical analysis, which is not computationally intensive like most of the feature importance techniques. The main by-products are an overall mean classification rate per feature and the rate aggregated by cluster type.

Following the instructions of Fränti and Sieranoja (2019), conducted seven experiments, by including and removing input features, to check the model robustness and how much information is changed or maintained accordingly. The findings of the PCA, EDA and feature importance were used to reason the selection of specific features.

Step 4 – model evaluation

To validate the refined model results, we used three qualitative assessments. First, we visually inspected the results using Google Earth and Street View images for internal validation, followed by a statistical assessment with auxiliary reference datasets (AGSN polygons, output of a model trained only with census-derived features and the municipal land use map). Next, we conducted an external validation with a local specialist through a semi-structured interview, where we had a guided discussion on the modelled spatial patterns, the major influencing factors and the applicability of the approach.

Results

In this section, we present the results of the proposed methodology of capturing intra-urban deprivation at city scale. We illustrate the selected spatial unit of analysis and moving window size, as well as the EDA findings. Then, we provide the outputs of the refined k-means model and the assessment of such findings.

Selection of spatial unit of analysis

The inspection of the AGSN polygons (Figure 6) indicates that 10 m and 100 m resolution only depict well the extreme settlement sizes, i.e. the minority of settlements. The medium settlement size is more appropriately depicted by 20 and 50 m grid sizes, however, 20 m performs better for elongated settlements. For the choice between 10 and 20 m grid size, we considered the reduction of computational time and the avoidance of unnecessary identifiable details of vulnerable communities. These considerations led to the conclusion that a small grid size (10 m) does not necessarily benefit the model, especially referring to the original resolution of the input datasets (Gram-Hansen et al., 2019). Thus, 20 × 20 m was chosen as the unit of analysis for São Paulo.

Figure 6. Exemplary AGSN polygons of deprived settlements overlayed with varying grid sizes: from left to right 100m, 50m, 20m, 10m.

Next, we compared the different sampling thresholds to maximise the pixels’ homogeneity. We chose to sample the cells that have at least 50% of their area intersected by the AGSN layer. Their total area is 59,978,800 m², which is the closest to the AGSN polygons area of 60,244,428 m². Figure 7 shows that the sampled grid of 20 × 20 m can depict even irregularly shaped deprived settlements, while avoiding most non-deprived information that could lead to misinterpretation of the spatial features.

Figure 7. Comparison between the grid base layer before and after the sampling process (threshold of 50%).

Selection of moving window sizes: PCA approach

As described in the Materials and Methods section, the feature extraction process created 35 GLCM texture features and 20 LSM features. Several PCAs were run to select moving window sizes; results are documented in detail at the open repository⁵ (Trento Oliveira, 2021). For the GLCM features, the PCA metrics indicated the 5 × 5 window size expressing less correlation among features and higher loading scores when compared to the other kernels. The metrics also show second moment and contrast features have the lowest communality values and account for a small proportion of the model variation.

For the LSM features the 5 × 5 moving window is chosen due to the fewer collinearity patterns when compared to the other window sizes. We also decide to remove the Aggregation Index (AI) and Fractal Index (FRAC) features of the model due to their very high correlation within each other and among the other features.

Exploratory data analysis

Exploratory data analysis consisted of generating boxplots and exclude outliers based on them, plotting histograms and correlation matrix.⁶ The results provide insights on input patterns such as that most data points are close to service and infrastructure facilities, but a few outliers are very far from the provision (Figure 8). It also shows the similar distribution and high collinearity of the spectral bands, which justify the model experiment without them. Considering the possible model improvements for dimensionality reduction, one experiment was also done excluding features with high collinearity (|r|>0.8).

Figure 8. EDA histograms created in R displaying input data patterns through descriptive statistics.

K-Means implementation

For model optimisation, we performed seven experiments, checking the model consistency and the refinement of the results. As shown in Figure 9, the tests include the following input: (1) only with GIS-based features; (2) only RS-based features; (3) both GIS- and RS-based features; (4) without Sentinel 2A spectral bands; (5) without Sentinel 2A but including Band 5; (6) only the selected features from Test 4; (7) including census-derived features.

Figure 9. Flowchart of model optimisation experiments.

The first two models tested the sensitivity of the model to the input, thus, they were initialised 20 times (Karmitsa et al., 2018). In each initialisation, a very similar result appeared – considering the number of clusters and spatial cluster patterns – which indicates good model robustness, especially when combined with feature importance scores. When GIS- and RS-based features were included in Model 3, the significance of the features remained similar to the results in Models 1 and 2. However, feature importance tools often indicate input features with high collinearity as highly significant to the model (Nugrahita & Surjandari, 2020), which can generate biassed results. For this reason, we trained Model 4 without all spectral bands and Model 5 with only band 5, as band 5 indicates an estimation of vegetation chlorophyll content and vegetation is present within deprived settlements in São Paulo. The illustration A in Figure 10 shows that both models have less pixelation effect, likely derived from the spectral information, but Model 4 also can detect non-residential land use types within the settlements which is very significant information.

Figure 10. Comparison of model outputs.

Model 6 tests the value of data redundancy on the K-Means algorithm, by removing the highly correlated features and the lowest misclassification rates. The result shows that the algorithm requires detailed information to detect finer differences between cluster types, hence highlighting morphological heterogeneities at the intra-urban level. Illustration B in Figure 10 shows how Model 4 captures the non-residential land uses that are not detected by Model 6, mainly due the removal of contrast and homogeneity texture features that identify fragmentation of the urban fabric.

Finally, we ran a model combining the input features of Model 4 with 12 census-derived features (see Appendix 2), to assess the impact of adding household-level information that are not depicted by the EO-derived features. Despite the temporal mismatch of the AGSN boundaries (2019) and the census data (2010), it adds important levels of deprivation, such as average income and homeownership. However, Model 7 cannot capture specific land uses and it detects the discrete census tract boundaries – creating spatial fallacies in the gridded approach (Illustration C - Figure 10).

It is important to state that in all experiments, the elbow method was calculated using two Python packages for consistency. The resulting plots for Model 4, chosen as the optimal one for further assessment, indicates three eligible k values (Appendix 3). We mapped and inspected each of them visually for a qualitative interpretation as well.

K-Means results

São Paulo is differentiated into four clusters; and the model provides two main outputs: (1) a map with the generated cluster types and violin plots with feature values per cluster type (Appendix 4), that we condense in radar graphs using the mean feature values for easier visual comparison (Figure 11); (2) the misclassification scores, providing information on the relevance of each feature per cluster type (Figure 11).

Figure 11. Radar graphs and clustering map of emerged cluster types.

Cluster 0 is allocated to fewer cells, predominantly located in the outskirts of the city, with low accessibility to financial, educational and health services. Cluster 0 occupies elevated areas close to dense and healthy vegetation, indicated by the higher DEM and NDVI values. Lowest night-time light (VIIRS) values indicate both less population density and lower economic status. We visually inspected the cells to check whether the lowest values of mean and variance (texture features) reveal the presence of bare soil or slum-like areas and found several construction sites, hence, bare soil. Due to their peripheral unpopulated location, the cells of Cluster 0 present very low geometry and density values, which indicate the absence of buildings captured by the WSF layer. These features could be depicting small, scattered settlements in the periphery. Thus, we label this cluster as “Infant Settlements in Open Spaces”.

Cluster 1, located near non-built-up and rural areas, reflects higher deprivation levels than Cluster 0, because of the higher population counts, the even further distance from main roadways, highest density of poor-quality roads and highest slope values. High levels of homogeneity and entropy are controversial, indicating both spatial uniformity and disorder, which has to do with the organic landscape layout. The geometry features present fragmentation of the urban fabric and irregular and sprawled structures, such as settlements in an early development stage. For this reason, we named the cluster “Poorly Consolidated Settlements”.

Cluster 2 can be found in various locations around the city, often near roadways and riverways, clearly closer to services and infrastructure, assuring basic provision. The proximity to waterways is mostly referring to channelized rivers, since we only used major water veins in the study, even so, it might indicate higher flood risk. The lower vegetation values and high night-time light values can indicate industrial areas included in the sampling. Texture-wise, the high mean and variance features, combined with the lowest homogeneity values and skewed distribution of patch mean area, suggest more pixel edges and sharpening, indicating a formal layout such as non-residential land uses. Therefore, we nominated Cluster 2 as “Less deprived settlements connected to non-residential areas”.

Finally, Cluster 3 corresponds to almost half of the total cells and are mostly centrally located. It presents the highest population and dead-end street density values. The highest homogeneity and lowest entropy, contrast and dissimilarity values indicate pixel uniformity. Combined with the high mean texture values, the cells depict denser and compacted urban areas. The skewed distribution of the shape index and large urban patches indicate more complex settlements. Based on this, we entitled the cluster “Densely urbanised and mature settlements”.

The output of the feature importance tool is also very relevant for the analysis of the cluster types (Appendix 5). Several RS-based features have prevailing importance, but it makes sense as they are 2/3 of the features included. In light of the above interpretation, the highest importance of the geometry and vegetation features to Cluster 0 are quite meaningful. For Cluster 1, service provision scores highest, together with high DEM and nighttime light suggesting peripheral areas with lower income in elevated terrain. The inverse influence of texture features to Cluster 2 and 3 state the built-up structure of each cluster type: to 2, it refers to the presence of non-residential areas and to 3, evinces a slum-like morphology.

K-Means validation

Visual assessment

The visual assessment of the clusters shows clear morphological differences. Figure 12 shows two deprived areas, each containing two of the cluster types. The first area is located in the far north, catching much poorer accessibility and service provision. View A presents a low-density settlement within a highly vegetated area (Cluster 0) and View B shows a hilly terrain with precarious building structures (Cluster 1). The second area shows the Paraisópolis district, the largest and most populated deprived settlement in the city, close to the centre, thus provided with basic infrastructure. View C exhibits larger non-residential building sizes and even proper sidewalks for Cluster 2, while View D illustrates Cluster 3 with dense residential fabric, with multi-storey buildings.

Figure 12. Visual assessment of clusters using satellite and street-view images. Source: (Google Maps, retrieval date: 24^th May 2021).

Statistical assessment

We compared the distribution of cluster polygon sizes with the areal statistics of the AGSN polygons (Table 3). We defined the assignment of a cluster type to a settlement considering that at least 50% of the settlement area belongs to a certain cluster. Table 3 shows that the two peripheral clusters (0 and 1) have a significant difference. Cluster 0 is mostly allocated in large settlements – encompassing commission errors – while Cluster 1 is often assigned to small ones. This suggests another level of area deprivation, as the settlements are not only located in periurban areas, but also take place in a scattered way.

Table 3. Number of settlements per cluster type according to the polygons areal statistics – Quartile 1 (Q1): 0.7 ha. Median (M): 1.9 ha. Quartile 3 (Q3): 4.2 ha.

	AGSN settlement sizes
	Cluster	Tiny (area<Q1)	Small (Q1<area<M)	Average (M<area<Q3)	Large (Q3<area<Max)	Huge (area>Max)
Settlements per cluster type	0	2	0	3	7	7
	1	26	47	46	34	34
	2	177	126	100	49	24
	3	186	201	200	135	43
	Sum	391	374	349	225	108

Then, we trained a model with 12 census-derived variables (Appendix 2) to be used as reference data estimating the socioeconomic conditions of the settlements. Two clusters were derived from the elbow method, labelled “less deprived” and “more deprived” (workflow details in Appendix 6). Figure 13 illustrates the differences between the two clusters, where higher built-up density, average income and provision of sewage infrastructure and lower average household population refers to less deprived areas. Clusters 2 and 3 mainly overlay with less deprived cells, while half of Cluster 1 cells overlay the more deprived cluster. This indicates another level of precariousness to Cluster 1, already depicted as the most morphologically deprived among them.

Figure 13. Radar graph with mean features values of census-derived model and bar graph showing the distribution of each cluster type per census-based category.

Lastly, we overlaid the results with the land use layer. Initially with 30 land use classes, we aggregated to six: mix use, conservation areas, high-income residential, low-income residential, commerce/industry and transport corridors. Figure 14 shows the distribution of each cluster type per aggregated land use class, and it validates the results of our model. Cluster 0 mostly occurs in conservation areas but also in industrial areas (with bare soil). In relative values per cluster, Cluster 1 mostly present in low-income residential and conservation areas. Considering the hilly location, it might indicate more social vulnerability. Cluster 2 has a diverse distribution across land uses, which makes sense with its various locations, but the graph also indicates its significant predominance in non-residential classes.

Figure 14. Cells count per cluster type for each aggregated land-use class.

Expert evaluation

The results of the model were presented to a local expert on urban deprivation who confirmed the diversity of the deprived settlements in São Paulo and how well the resulting cluster types are qualifying these areas. She evinced the heterogeneity between central and peripheral areas, the existence of early-development stage settlements that face severe deprivation, depicted by Cluster 1. The fact that the model can capture these recently occupied regions with lack of service and infrastructure provision was found to be very important for local policies.

Another useful aspect of the clustering results is the differentiation of built-up and non-built-up areas, especially considering the errors of commission in the AGSN layer. She highlighted that the proximity of waterways and highways of Cluster 2 also reflect a contamination factor, e.g. open sewer and industrial pollution. The expert also stressed the existence of non-residential land uses within the AGSN boundaries and praised the ability of the model to identify these areas.

Discussion

Without discrediting the necessity of traditional binary approaches of deprived vs non-deprived areas as a starting point, we are proposing a step further by exploring the potential of machine learning models to capture their intra-urban diversity. We designed a methodology that uses a gridded unsupervised approach and solely open GIS- and RS-based data and software packages, dealing with scalability and transferability issues of state-of-the-art models.

By employing an unsupervised ML method, one does not require training and testing datasets which are computationally demanding and often unaffordable for LMICs, ensuring scalability potential. Using features derived from freely available datasets guarantees transferability potential, fighting current data poverty issues. This combination allows applicability to other contexts and provides a step forward in understanding deprivation at different scales. Yet, the model still requires tailoring from two main fronts: the selection and extraction of features and the spatial unit of analysis. Particularly to the use of open geodata, the cleaning and preparation process is quite demanding, once we inspected each dataset for quality issues in order to reduce modelling uncertainties.

This research throws light on the local manifestation of deprivation by developing a list of potential spatial features, using open data sources in an effort to bridge the gaps of census data and VHR imagery, considering spatial coverage, resolution and completeness. The eight dimensions conceptualised here to characterise the living conditions of neighbourhoods of São Paulo can be generalised to other study areas but require tailoring, not only regarding the inclusion or removal of a certain feature, but also their combined interpretability for planning decisions. For example, Cluster 3 shows a high density of dead-end streets and high built-up density which can increase costs for infrastructure implementation, while for Cluster 1 the intervention costs might be more affected by the local topography and sprawled distribution. However, it is important to state that the features extraction processes proposed here are time-consuming and not fully automated, which can hinder the practicality of the workflow.

The sensible selection of the unit of analysis is a major step in the methodology. In addition to avoiding ecological fallacies, it increases comparability with other reference datasets for decision-making, such as global gridded outputs or hazard maps for disaster risk plans. Besides the increasing number of recent slum studies employing gridded approaches, they rarely document how the unit of analysis was chosen. Moreover, most of the work uses VHR imagery and we present a promising result with 20 m resolution, sufficiently detailed for a city scale analysis. Nevertheless, the increasing availability of high-quality sources would be only beneficial for deprivation research.

It is important to stress that features interpretations can change depending on the adopted unit of analysis. Especially texture- and geometry-related features can provide very different depictions in a different spatial resolution. For example, with 1 m resolution, high entropy values indicate complex and dense morphology, while with a 20 m – one pixel encompassing 2 housing units on average – high entropy suggests higher variability between the pixel and its surrounding, which is associated with more formal areas.

Additionally, with respect to reproducibility, this research employs an efficient and simple algorithm, which results highly depend on the number of clusters and the input dimensionality. The outcomes here presented indicate that removing some features due to high collinearity or model significance can hamper the cluster’s separability. The usefulness of the feature importance tool is undeniable, especially for interpretability of model robustness and performance, but for dimensionality reduction the effects could be more considerable for larger scales, hence, larger datasets.

Conclusion

We have demonstrated that an unsupervised gridded ML approach, not yet exploited in literature for slum characterization studies, can capture the morphological dimensions of deprived areas and identify the intra-city differences. The multiple experiments and validation procedure show the good model performance and the significant contribution of input features for each resulting cluster type. The main advantage of the methodology is its transferability and scalability capacity, pushing the model application towards larger scales, thus, critical to global development goals.

We have shown that physical characteristics of such settlements can be derived from open GIS- and RS-based data and how well they can generate relevant information on the nuances of deprivation within the city. Despite our results, future work should not drop spectral features automatically, because the model is context-based and a good cluster separability can come from land cover details, e.g. in areas where detection of bare soil is important. Local specialists are fundamental for the decision on the inclusion or exclusion of certain features. They can provide insights during both the development of indicators and interpretation of clustering results.

In this sense, recently available land cover products could also be beneficial to the model results. Considering the environment of certain deprived areas, high-resolution land cover maps can improve the model’s ability to diffentiate between bare soil and other land cover classes (Dos Santos et al., 2022), between bare soil and roofing materials (Kuffer et al., 2016; Owen & Wong, 2013) and the quality of road infrastructure (Nobrega et al., 2008).We acknowledge the importance of a basic understanding of the study area, while the proposed feature extraction workflow remains the same.

In terms of applicability, this study can inform decision-makers for locally tailored policy-making such as: the location of pocket slums located in conservation areas (Cluster 0); the urgency of service provision for settlements in steep locations (Cluster 1); the presence of non-residential land uses in bordering zones of the AGSN layer (Cluster 2); and the predominance of complex and dense urban fabric (Cluster 3).

Nevertheless, we see room for further investigation and improvement that should primarily focus on upscaling the model to regional level and the full automation of the workflow. Future study should automate the model in an API environment that can process all datasets, without relying on the computational capacity of the user. This would ensure the full applicability of this pilot research to the municipalities. To ensure the full applicability of this pilot research to the municipalities, GIS experts with knowledge on the local urban systems are required for tailoring the model features and unit of analysis.

Acknowledgments

The authors acknowledge the support of Alexandra Pedro on the validation procedure, Dr. Raian Maretto and Dr. Flávia Feitosa for the technical consultancy on machine learning models, and Dr. Caroline Gevaert for special assistance on uncertainties sources.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The entire workflow, datasets, model input, code scripts and model output are archived and available at GitHub platform (https://github.com/ltrentooliveira/MSc_Archive) to ensure maximal replicability.

Additional information

Funding

This research received no external funding.

Notes

¹ https://github.com/ltrentooliveira/MSc_Archive/blob/main/PCA_processing.ipynb.

² https://github.com/ltrentooliveira/MSc_Archive/blob/main/EDA.R.

³ Machine used for the present study: LENOVO _MT_20QU_BU_Think_FM_ThinkPad P1 Gen2 Intel Core i7-9750 H CPU, 16GB RAM, Windows 10.

⁴ https://github.com/ltrentooliveira/MSc_Archive/blob/main/kmeans_model.ipynb.

⁵ https://github.com/ltrentooliveira/MSc_Archive/blob/main/PCA_processing.ipynb.

⁶ https://github.com/ltrentooliveira/MSc_Archive/blob/main/EDA.R.

Appendix 1.

SQL expressions to GIS attributes selection

Table

Feature name	Data sources	SQL expression
Distance from major roads	OSM	“type” = “primary” OR “type” = “primary_link” OR “type” = “secondary” OR “type” = “secondary_link” OR “type” = “service” OR “type” = “services” OR “type” = “tertiary” OR “type” = “tertiary_link” OR “type” = “trunk” OR “type” = “trunk_link”
Density of poor-quality roads	OSM	“surface” = “paving_stones” OR “surface” = “unpaved” OR “surface” = “dirt” OR “surface” = “gravel” OR “surface” = “mud” OR “surface” = “pebblestone” OR “surface” = “pedras_para_pavimentação” OR “surface” = “sand” OR “surface” = “ground” OR “surface” = “earth” OR “surface” = “compacted” OR “surface” = “fine_gravel”
Distance to education facilities	HOTOSM	“amenity” = “school” OR “amenity” = “university” OR “type” = “amenity” OR “amenity” = “kindergarten”
Distance to health facilities	HOTOSM	“amenity” = “clinic” OR “amenity” = “clinica” OR “type” = “hospital”

Appendix 2.

Acquired by IBGE, we derive features 01 to 04 from table ‘Básico_SP’; feature 05 from table ‘DomicílioRenda_SP’ and feature 06 to 12 from table ‘Domicilio01_UF’

Table

Initial	Feature name	Description	Processing
F01	POP_DENS	Population density	V002/total tract area
F02	AV_HHPOP	Average number of residents per households	V003 (no processing)
F03	BU_DENS	Dwelling/building density	V001/total tract area
F04:	AV_INC	Average monthly income of head of household per dwelling unit (V005)	V005 (no processing)
F05	INC_MWAGE	Proportion of households with monthly income per capita of up to 1 minimum wage	(V005 + V006 + V007 + V008)/V001
F06	HH_APT	Proportion of apartment-type dwellings	V005/V002
F07	HH_PRIV	Proportion of private owned households	(V006 + V007)/V002
F08	HH_GARB	Proportion of dwelling units with adequate waste disposal	V035/V002
F09	HH_WATER	Proportion of dwelling units with adequate water provision	V012/V002
F10	HH_SEWAGE:	Proportion of dwelling units with adequate sewage network	V017/V002
F11	HH_PITLATRINE	Proportion of dwelling units with sewage discharged in pit latrine	V019/V002
F12	HH_NOSEWAGE:	Proportion of dwelling units with sewage discharged in septic in river or lake	V021/V002

Appendix 3.

Elbow method calculation from different implementation packages. Left: from sklearn package; Right: from yellowbrick.cluster package.

Appendix

Appendix 5.

Visualisation of features using FeatureImpCluster tool.

Appendix 6.

Visualisation of features using FeatureImpCluster tool.