Estimating the Multidimensional Poverty Index in Bogotá (Colombia) Using Satellite Imagery

Abstract

This paper proposes an alternative approach to estimating poverty levels using satellite imagery from Bogotá, Colombia, and surrounding municipalities. Two modeling strategies are compared: Convolutional Neural Networks (CNNs) and a newly proposed extension of the Generalized Tensor Regression model based on spatially-weighted aggregations, denoted SWGTR. The models are evaluated using multiple performance metrics, including Average Accuracy (AA), Overall Accuracy (OA), Area Between Curves (ABC), Kappa Index (KI), Mean Absolute Error (MAE), and Cumulative Sums (CS), to determine their effectiveness in capturing the spatial distribution of the Multidimensional Poverty Index (MPI) at the pixel level. The application leverages PlanetScope (PS) satellite images alongside block-level data from the 2018 National Population and Housing Census (CNPV2018). The SWGTR approach demonstrates competitive performance across several metrics, particularly those accounting for the ordinal structure of poverty levels. In addition, we illustrate how the probabilistic outputs of SWGTR can be used to construct pixel-level uncertainty measures, which may inform data collection strategies by identifying regions where predictions are less certain. Overall, the findings suggest that combining spatial modeling with remote sensing data provides a promising avenue to complement traditional survey-based methods, while highlighting important limitations and directions for future research.

1. Introduction

If we conduct a general analysis of all the countries in the world, it is not difficult to conclude that one of the problems that affects all of them in a transversal way is the high incidence of poverty, both in its multidimensional and monetary dimensions (Oxford Poverty and Human Development Initiative and United Nations Development Programme 2022). Out of the total of 111 countries considered as developing, with a total population of 6.1 billion people, approximately 1.2 billion of them (equivalent to 19.1% of the total) live in multidimensional poverty (Alkire et al. 2021).

An alarming 83% of the population in multidimensional poverty is concentrated in regions such as Sub-Saharan Africa, with 579 million people, and South Asia, with 385 million people. In Sub-Saharan Africa, nearly 3.8 million people are in total poverty, which means that they suffer a total deprivation of all components included in the definition of multidimensional poverty (Lederer 2024). All of this indicates that, on a global scale, the situation is quite critical. This is also reflected in Latin America, where most countries are developing, and the outlook is not encouraging. Among these, the situation in Haiti stands out, as it is, unfortunately, the country with the highest percentage of population living in poverty. This widespread poverty situation in the less developed regions of the world is derived from the fact that, on the one hand, the efforts focused on policies to measure, control, and reduce poverty are not precise and strong enough, creating a gap when compared to developed economies; and, on the other hand, the criteria for wealth distribution policies tend to favor a few and disadvantage the majority (Berry 2003).

In Colombia, poverty measurement is not a recent issue. Since the late 1980s, various methodologies have been adopted for its evaluation and use in public policy design. Currently, there are two official poverty measurements: monetary poverty and multidimensional poverty. These indicators are governed by the country’s official methodologies (DANE-DPS-DNP 2012) and are periodically processed by the National Statistical Office of Colombia, commonly known as DANE, the acronym of its name in Spanish.

For the purposes of this work, it is relevant to focus on the second measurement of poverty, multidimensional poverty, which has led to the construction of the Multidimensional Poverty Index (MPI). This is based on the methodology of Alkire and Foster (2011) from the Oxford Poverty and Human Development Initiative (OPHI), in which the index is constructed taking into account the five dimensions that encompass aspects related to quality of life: household educational conditions, childhood and youth conditions, health, work, access to public services, and housing conditions. Each dimension has the same weight (20%), and these are analyzed and quantified through fifteen indicators, which have equal weight in each dimension.

The search for a measure that meets the need for a comprehensive and updated multidimensional poverty index is a longstanding concern, leading to the MPI being calculated using different methodologies throughout the country’s history. Among these, three complementary efforts are worth mentioning: (i) the calculation of the MPI based on the National Quality of Life Survey (ENCV, for its acronym in Spanish; Departamento Administrativo Nacional de Estadística [DANE] 2024); (ii) the calculation of the MPI using information from the National Population and Housing Census (CNPV, for its acronym in Spanish; Departamento Administrativo Nacional de Estadística [DANE] 2018); and (iii) an approximation of the MPI using satellite images and light-night intensity.

The first effort has achieved a comprehensive national geographic representation in the most recent editions of the ENCV survey. The second effort, based on census information, has aimed to achieve a higher level of geographical disaggregation. Finally, the third poverty measurement effort has sought to generate an estimate close to the official MPI poverty indicator, using as input satellite images and nighttime light intensity information. This is intended to reduce the high cost associated with the two first efforts, primarily due to (a) the complexity of the sampling schemes and operational logistics required to provide measurements at different levels of geographic disaggregation, and (b) the need to collect data at a high frequency to ensure updated statistics on the country’s poverty situation.

For this third effort, images from 2018 of the Sentinel-2 satellite, which can be downloaded from Google’s application: Google Earth Engine (GEE), were considered. The complete methodology can be found in Head et al. (2017) and Jean et al. (2016), but its objective can be summarized as a way to estimate poverty levels in five African countries using Convolutional Neural Networks (CNN). This was achieved by using nighttime light intensity information grouped into three levels (Low, Medium, and High) as the response variable and satellite images as covariates in the model, thus classifying each image into one of the three levels. In the end, the goal was to compare the results of this model with the official estimates.

Additionally, works such as Steele et al. (2017) and Heitmann and Buri (2019) have reported the use of similar approaches for measuring poverty in Bangladesh, Ghana, and Uganda, where official statistics derived from surveys or measurements like the Poverty Probability Index, are used as approximations of the dependent variable. In turn, the inputs range from household surveys, Remote Sensing (RS), Call Detail Records (CDR) at spatial scales that include the political-administrative divisions of each country, to Voronoi polygons constructed with the location and coverage of cell phone towers.

At this point, it is evident that there are promising opportunities in the generation of poverty-related statistics. On one hand, there is the possibility of combining various data sources, both from traditional sources and administrative records, as well as information extracted from satellite images. On the other hand, given the high cost associated with conducting exhaustive surveys, increasing the sample size of a survey to cover all geographic areas and different population groups might not be a practical option.

Therefore, the relevance and urgency of having accurate and timely poverty statistics at more detailed disaggregation levels make it necessary for these surveys to be conducted at least every two or three years. In this sense, the incorporation of non-traditional data sources has the potential to overcome the inherent limitations of conventional data sources in terms of detail and timeliness in poverty estimation.

Recent work has demonstrated that deep learning applied to Earth observation and street-level imagery can recover fine-scale socioeconomic patterns, including poverty and vulnerability, intra-urban population distribution, and neighborhood-level income and education. Convolutional Neural Network (CNN)-based models trained on satellite or street-view imagery have been used to predict a variety of census-like indicators at neighborhood or census-sector scales, often achieving strong agreement with survey-based labels and showing that visual urban form encodes social and economic conditions (Ben Abbes et al. 2024; Machicao et al. 2022; Runfola et al. 2024). Parallel research on population mapping and spatial disaggregation has leveraged fully convolutional or tabular deep architectures and gradient-boosted trees to downscale or impute counts for multiple demographic or socioeconomic groups, but typically over regular grids or units of more uniform size and for relatively small sets of targets (Liu et al. 2023a; Runfola et al. 2024). Multimodal approaches that combine satellite with street-level imagery, points of interest, mobility, or text data further improve predictions of income and other socioeconomic indicators in cities, often via contrastive or attentional fusion mechanisms (Liu et al. 2023b; Machicao et al. 2022; Yong and Zhou 2024). However, much of this literature assumes fixed-size image patches and targets a narrow set of outcomes, limiting applicability when administrative units vary dramatically in spatial extent or when many census variables must be estimated simultaneously. Runfola et al. (2024) address this by introducing a multi-glimpse recurrent attention architecture designed for extreme scope variance, simultaneously predicting fifty-two census variables across Mexican municipalities whose areas span several orders of magnitude. Related work on post-hoc spatial-imprecision adjustment, deep-learning reviews in economics, and open-source frameworks for wealth estimation from satellite imagery further suggests a broader paradigm in which remotely sensed imagery and advanced deep learning can provide scalable, multi-indicator estimates of socioeconomic conditions across heterogeneous spatial units (Baier and Runfola 2025; Ben Abbes et al. 2024; Dell 2025; Zheng et al. 2023).

In addition to these approaches, our work focuses on two innovative methodologies that are comparatively evaluated to estimate a poverty indicator at the household level. This indicator is based on the use of the aforementioned non-traditional data sources, such as satellite images, in conjunction with traditional sources such as surveys, censuses, and geographic data. In this context, the challenge addressed involved the formulation of two models with a common goal: using unstructured information (images) to characterize the population at different poverty levels, considering that the underlying methodology in each of these models is different. In conclusion, a detailed comparison of both models was conducted to determine their advantages and disadvantages, in order to guide their future implementation effectively. The following is an overview of each model.

The first model uses Convolutional Neural Networks (CNN) for pixel classification or image segmentation and has as response variable a set of binary masks or images with dimensions of 320 × 320 × 5. In this scenario, 320 refers to the mask’s dimensions, that is, the height and width values, which were chosen to ensure the quality of the resolution of each chip. On the other hand, the value five represents the number of poverty levels considered in this study. In these masks, each pixel can only take values of zero or one, indicating whether a pixel belongs to a certain level of poverty.

As for the explanatory or independent variables, chips or crops from each of the satellite images of Planet Scope are used, where each crop has a dimension of 320 × 320 × 4. In this case, four channels are generated, materializing in the usual three RGB (R: red, G: green, B: blue) channels, plus a near-infrared (NIR) channel. The inclusion of the RGB channels contributes to capturing detailed information about the color spectrum in the satellite images, allowing for a more complete analysis of the visualized information. It is worth noting that each original Planet Scope image covers approximately 120 km².

The second model proposal is based on the methodology presented by Liu et al. (2019), which introduces a Generalized Tensor Regression (GTR) approach for the classification of hyperspectral images. This approach originates from a simple and efficient classifier known as the ridge regression (RR) model, whose aim is to classify each unit of analysis into one of the categories of interest. These categories can be organized in a label matrix where each row indicates the category to which each individual belongs. In this way, the GTR model incorporates CANDECOMP/PARAFAC (CP) decomposition to sparsely represent the information available via satellite images as tensors (these concepts are detailed in Section 2), which reduces the number of parameters to estimate (Kolda and Bader 2009). Furthermore, this model takes advantage of existing knowledge about tensor properties, particularly non-negativity, to enhance the model’s discriminative capacity.

Among the advantages that distinguish this proposed approach, the following stand out: firstly, prior knowledge of tensor characteristics is considered, and non-negativity is exploited to promote the parsimony of the GTR model. Second, GTR gathers all the categories and trains them together. This feature is different from other rank regression methods to which GTR belongs. For example, the approaches proposed by Guo et al. (2011) and Hou et al. (2012) use binary classifiers. In Hou et al. (2012), multiple classification is achieved by training each class separately.

Since the nature of the poverty variable is ordinal, this work focuses on contributing to the field of ordinal classification by adjusting algorithms and/or metrics present in the existing literature. In this line, the proposed contribution focuses on adjusting existing nominal classification algorithms in multiple categories to incorporate ordinal classification more effectively. By addressing this issue, we aim to improve the ability of the proposed models to handle classification problems where the labels have an inherent order, which can have relevant applications in various fields, such as poverty estimation, agriculture, land use, risk management in rural areas, or data science in general.

In summary, the main contributions of this work are as follows:

• This work enhances the performance of poverty estimation models by assessing the added value of methodologies incorporating non-traditional data sources. In this way, the intention is to provide a possibly more cost-effective perspective to complement traditional survey-based methodologies. This ultimately allows for greater detail in poverty statistics.

• It also expands the literature on Generalized Tensor Regression (GTR) by introducing a spatially weighted variant that accounts for the different models trained for each subregion (chip) within the segmentation of the total area under analysis. This variant allows GTR to be competitive in terms of classification power to complex CNN architectures. As GTR usually provides a predictive model for each chip it is trained for, the proposal in this work is to construct spatially-dependent weights to combine the information of nearby models into the predictions of the target chip. The details are presented in Subsubsection 5.2.1.

• It provides a comprehensive set of existing and newly proposed mechanisms to assess ordinal classification tasks. Previous works have introduced different tools to assess the performance of ordinal classification models. This work gathers them altogether to assess the quality of the proposed models, and introduces the ABC indicator, analogous to the AUC indicator for nominal classification, but adapted to account for ordinality.

The rest of this manuscript is structured as follows: Sections 2 and 3 present the conceptual framework, detailing the underlying Tensor Regression Models, including Ridge Regression for multi-label classification and Generalized Tensor Regression (GTR), as well as Convolutional Neural Networks (CNN) with a focus on architectural variants like U-Net, U-Net Attention, and DeepLabV3+, along with the loss functions employed. Section 4 outlines the data sources utilized in the study. Section 5 describes the methodology, including the configurations for CNN and Spatially-Weighted GTR (SWGTR) models, and the metrics used for evaluating multi-category classification performance. Section 6 provides an in-depth analysis and comparison of the results, contrasting different CNN architectures and comparing CNN with SWGTR strategies using both a single test chip and the entire testing dataset. Finally, Section 7 presents the conclusions and suggests avenues for future research.

2. Conceptual Framework: Tensor Regression Models

This section presents the tensor regression concepts used for this work. The notions and definitions presented below were extracted and adapted from Liu et al. (2019).

2.1. Ridge Regression for Multi-Label Classification (RR)

Ridge Regression (RR) is a promising solution for the multiclass classification problem with a vector of attributes in ${\mathbb{R}}^{I_1}$ . We consider a multichannel image with $C$ classes and a training sample of $I_2$ units, represented as $\mathbf{X}\in {\mathbb{R}}^{I_1\times I_2}$ , where $I_1$ is the number of channels and $I_2$ is the number of observations. Each training sample ${\mathbf{x}}_{:,i}\in {\mathbb{R}}^{I_1}$ is assigned to one of the $C$ classes, with its label denoted by the vector ${\mathbf{t}}_{:,i}=[0,\dots ,1,\dots ,0{]}^T\in {\mathbb{R}}^C$ , indicating with a “1” the position of the corresponding class. The goal of RR is to estimate the coefficients $\mathbf{B}\in {\mathbb{R}}^{C\times I_1}$ that model the linear relationship between ${\mathbf{x}}_{:,i}$ and ${\mathbf{t}}_{:,i}$ by solving the following optimization problem

$$\underset{\mathbf{B}}{\min }\parallel \mathbf{T}-\mathbf{BX}{\parallel }_F^2+\lambda \parallel \mathbf{B}{\parallel }_F^2,$$

(1)

where the first term represents the fidelity term and the second term is the regularization term. Here, $\mathbf{T}:=[{\mathbf{t}}_{:1},\dots ,{\mathbf{t}}_{:I_2}]\in {\mathbb{R}}^{C\times I_2}$ , $\lambda >0$ is a regularization parameter, and $\parallel ·{\parallel }_F$ denotes the Frobenius norm. The problem in Equation (1) is convex and smooth, and its solution can be expressed in closed form as $\hat{\mathbf{B}}={\mathbf{TX}}^T{({\mathbf{XX}}^T+\lambda \mathbf{I})}^{-1}$ , where $\mathbf{I}$ denotes the identity matrix. We can then determine the label of a test sample ${\mathbf{x}}^{*}\in {\mathbb{R}}^{I_1}$ as expressed in Equation (2), where ${\delta }_c(·)$ denotes the $c$ -th element of the vector,

$$\hat{t}=\arg \underset{c\in \{1,2\dots ,C\}}{\max }{\delta }_c\left(\hat{\mathbf{B}}{\mathbf{x}}^{*}\right).$$

(2)

2.2. Generalized Tensor Regression (GTR)

The model in Equation (1) can be extended from a vector space for the features into a tensor space by solving the following optimization problem for each class

$$\underset{\{{\mathcal{B}}_{c:}\}}{\min }\sum _{c\in \{1,2\dots ,C\}}\sum _{j=1}^N{(t_{\mathrm{cj}}-〈{\mathcal{X}}_{:j},{\mathcal{B}}_{c:}〉)}^2+\lambda \sum _{c\in \{1,2\dots ,C\}}\parallel {\mathcal{B}}_{c:}{\parallel }_F^2,$$

(3)

where ${\mathcal{X}}_{:j}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_M}$ contains the input variables of the $j$ -th training unit, and ${\left\{{\mathcal{B}}_{c:}\right\}}_c$ is the set of weight tensors of the same dimension as ${\mathcal{X}}_{:j}$ used to learn classification patterns for all the classes $c\in \{1,2,\dots ,C\}$ . In this case, the dimension of the parameter space is much larger than the size of the training set, that is, $\underset{m=1}{\overset{M}{\Pi }}{I}_m>>N$ . As a result, it is not possible to obtain estimates in the fidelity term of Equation (3) due to the over-dimensionality problem that arises when there are more parameters than sampling units. To simplify the model and leverage the structural information in ${\mathcal{X}}_{:j}$ , it is assumed that the weight tensor ${\mathcal{B}}_{c:}$ for each class is of low rank. Specifically, ${\mathcal{B}}_{c:}$ is assumed to be decomposable using a rank-one CANDECOMP/PARAFAC (CP) decomposition (Kolda and Bader 2009) as ${\mathcal{B}}_{c:}={\mathbf{b}}_{c:}^{\left(1\right)}°{\mathbf{b}}_{c:}^{\left(2\right)}°\cdots °{\mathbf{b}}_{c:}^{\left(M\right)}$ , with ${\mathbf{b}}_{c:}^{\left(m\right)}\in {\mathbb{R}}^{I_m}$ and ° denoting the exterior product. This reduces the number of parameters in problem (3) to $C\sum _{m=1}^M{I}_m$ , but each class has its own rank-one tensor of weights.

At this point, it is evident that both the $\mathrm{RR}$ and its tensor version (Equations (1) and (3), respectively) combine multiple binary classifiers using a one-versus-all strategy. In other words, their training process can be treated as training $C$ classifiers for $C$ classes separately. One way to jointly link all separate classifiers is to use the factorization ${\mathcal{B}}_{c:}={\mathbf{b}}_{c:}^{\left(1\right)}°{\mathbf{b}}^{\left(2\right)}°\cdots °{\mathbf{b}}^{\left(M\right)}$ , with only the first vector, ${\mathbf{b}}_{c:}^{\left(1\right)}$ , being class-dependent; thereby further reducing the number of parameters to ${\mathrm{CI}}_1+\sum _{m=2}^M{I}_m$ . For problem (3), we have the following formulation by substituting the CP decomposition into the fidelity and regularization terms as shown in Equations (4) and (5), respectively.

$$\begin{array}{c}\multicolumn{1}{c}{f\left({\mathcal{B}}_{c:}\right)=\sum _{c,j}{(t_{\mathrm{cj}}-〈{\mathcal{X}}_{:j},{\mathcal{B}}_{c:}〉)}^2=\sum _{c,j}{(t_{\mathrm{cj}}-〈{\mathcal{X}}_{:j},{\mathbf{b}}_{c:}^{\left(1\right)}°{\mathbf{b}}^{\left(2\right)}°\cdots °{\mathbf{b}}^{\left(M\right)}〉)}^2}\\ \multicolumn{1}{c}{=\sum _{c,j}{(t_{\mathrm{cj}}-{\mathcal{X}}_{:j}{\times }_1{\mathbf{b}}_{c:}^{\left(1\right)}\underset{m=2}{\overset{M}{\Pi }}{\times }_m{\mathbf{b}}^{\left(m\right)})}^2}\\ \multicolumn{1}{c}{=\parallel \mathcal{T}-\mathcal{X}\underset{m=1}{\overset{M}{\Pi }}{\times }_m{\mathbf{B}}^{\left(m\right)}{\parallel }_F^2.}\end{array}$$

(4)

$$\begin{array}{c}\multicolumn{1}{c}{g\left({\mathcal{B}}_{c:}\right)=\sum _c\parallel {\mathcal{B}}_{c:}{\parallel }_F^2=\sum _c\parallel {\mathbf{b}}_{c:}^{\left(1\right)}°{\mathbf{b}}^{\left(2\right)}°\cdots °{\mathbf{b}}^{\left(M\right)}{\parallel }_F^2}\\ \multicolumn{1}{c}{=\sum _c\sum _{i_1,i_2,\dots ,i_M}{(b_{{\mathrm{ci}}_1}^{\left(1\right)}{b}_{i_2}^{\left(2\right)}\dots b_{i_M}^{\left(M\right)})}^2}\\ \multicolumn{1}{c}{=\sum _c\parallel {\mathbf{b}}_{c:}^{\left(1\right)}{\parallel }_2^2·\parallel {\mathbf{b}}^{\left(2\right)}{\parallel }_2^2\dots \parallel {\mathbf{b}}^{\left(M\right)}{\parallel }_2^2}\\ \multicolumn{1}{c}{=\underset{m=1}{\overset{M}{\Pi }}\parallel {\mathbf{B}}^{\left(m\right)}{\parallel }_F^2.}\end{array}$$

(5)

In this formulation, $\mathcal{T}\in {\mathbb{R}}^{C\times 1\times 1\times N}$ is a tensor of order $(M+1)$ gathering all the $\left\{t_{\mathrm{cj}}\right\}$ elements, $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_M\times N}$ contains the input features of all the units in the training set. The $c$ -th row of ${\mathbf{B}}^{\left(1\right)}\in {\mathbb{R}}^{C\times I_1}$ is equal to ${\mathbf{b}}_{c:}^{\left(1\right)}$ , and ${\mathbf{B}}^{\left(m\right)}\in {\mathbb{R}}^{1\times \mathrm{Im}}$ is the transpose of ${\mathbf{b}}^{\left(m\right)}$ , $m\in 2,\dots ,M$ . The symbol ${\times }_m$ denotes the $m-$ mode product. Consequently, Equation (3) can be expressed as

$$\underset{{\mathbf{B}}^{\left(m\right)}}{\min }\parallel \mathcal{T}-\mathcal{X}\underset{m=1}{\overset{M}{\Pi }}{\times }_m{\mathbf{B}}^{\left(m\right)}{\parallel }_F^2+\lambda \underset{m=1}{\overset{M}{\Pi }}\parallel {\mathbf{B}}^{\left(m\right)}{\parallel }_F^2.$$

(6)

In the context of the model from Equation (6), it is necessary to impose a constraint on the possible values that ${\mathbf{B}}^{\left(m\right)}$ can take. In this case, ${\mathbf{B}}^{\left(m\right)}$ is restricted to only positive values. However, this constraint leads to the issue of “sparsity,” referring to a matrix that contains a significant number of zeros. This sparsity can be advantageous, as it reduces the model’s complexity and decreases computation times. In the case of tensor classification, the non-negativity assumption on the weight matrices ${\mathbf{B}}^{\left(m\right)}$ enables the model to capture generalized relationships in the input tensors, particularly those related to rank-one components obtained from the CP decomposition.

After running an optimization routine similar to the RR model (Liu et al. 2019), a test sample ${\mathcal{X}}^{*}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_M}$ , can be classified analogously by

$$\hat{t}=\arg \underset{c\in \{1,2\dots ,C\}}{\max }{\delta }_c\left({\mathcal{X}}^{*}\underset{m=1}{\overset{M}{\Pi }}{\times }_m{\hat{\mathbf{B}}}^{\left(m\right)}\right),$$

(7)

where ${\hat{\mathbf{B}}}^{\left(1\right)}={[{\hat{\mathbf{b}}}_{1:}^{\left(1\right)},\dots ,{\hat{\mathbf{b}}}_{C:}^{\left(1\right)}]}^T$ and ${\hat{\mathbf{B}}}^{\left(m\right)}={\left({\hat{\mathbf{b}}}^{\left(m\right)}\right)}^T$ , $m\in \{2,\dots ,M\}$ .

3. Conceptual Framework: Convolutional Neural Networks (CNN)

The second modeling approach considers the use of Convolutional Neural Networks (CNNs), which are a particular instance of artificial neural networks (ANNs) suitable for learning tasks involving high-dimension objects, such as images.

The article by O’Shea and Nash (2015) offers a succinct introduction to ANNs, which are computational systems inspired by the functioning of biological nervous systems like the human brain. ANNs are composed of numerous interconnected units known as artificial neurons. Unlike biological neurons that connect through dendrites and synapses in a complex, often arbitrary, manner; ANNs use weighted connections defined by parameters learned during training. Biological neurons process input signals collectively to produce an output, and ANNs mimic this mechanism through weighted links between artificial neurons. These adjustable weights allow the network to learn patterns and adapt to specific tasks. While structurally different from biology, the core principle remains: learning and task execution based on interconnected processing units.

The convolution operation is a core component of CNNs, as it enables the extraction of relevant features from images through the use of filters or kernels. Conceptually, convolution acts as a weighting function applied over an input signal to emphasize local patterns. When applied to multidimensional arrays, such as images, the kernel slides across the input and computes local linear combinations at each position. This operation can be defined in both continuous and discrete domains and naturally extends to 3D tensors, like RGB images. By applying multiple convolutional layers with various kernels, the network is capable of constructing increasingly abstract feature representations.

Two key concepts that refine the convolution process are padding and stride. Padding preserves the spatial dimensions of the input tensor by adding zeros around the edges, while stride controls the step size of the kernel, affecting the resolution of the output. In this work, we opt for a stride of one and no padding, which reduces the spatial dimensions progressively while ensuring full coverage of the input image to capture relevant patterns. The general form of convolution is expressed as a triple summation over spatial and depth indices, making the operation adaptable to complex multidimensional data.

To introduce non-linearity into the model, activation functions such as sigmoid, softmax, and ReLU are used. The sigmoid function is suitable for binary classification, while softmax is used in multiclass settings as it outputs normalized probability distributions. ReLU, the most widely adopted activation function, improves efficiency by activating only a subset of neurons—those with positive inputs—resulting in sparse representations. Additionally, pooling operations are applied to reduce the dimensionality of feature maps while retaining dominant characteristics. Max pooling and average pooling are the most common variants, contributing to invariance to small shifts and faster processing.

3.1. Architectures

An important step in an artificial neural network (ANN) is determining its architecture, which defines the number of hidden layers and their connections. The literature on ANN architecture design is quite extensive in the possible architectures that can be considered, suggesting that the definitive choice depends on the specific problem that we are facing. However, there are guidelines for creating optimal architectures. For example, Goodfellow et al. (2016) suggest defining the architecture based on the number of hidden layers (depth) and the width of each layer (width). Additionally, while high values for these two parameters can lead to good results, they also impose significant computational costs, making the architecture cumbersome to use. Therefore, selecting the structure should be done experimentally by monitoring the error on a test dataset to avoid overfitting.

This study employs three architectures, each designed for image segmentation. Their primary goal is to analyze images and detect objects of interest (e.g., bodies of water, forests, roof materials) by identifying their presence and precise location. A description of the architectures used is described below.

3.1.1. U-Net

Most convolutional networks assign a single label to an entire image, lacking pixel-level precision, which is an essential feature in areas like medical imaging. While early segmentation techniques existed, a major breakthrough came with U-Net (Ronneberger et al. 2015), which significantly improved medical image segmentation. U-Net is now widely applied in domains such as land use classification and autonomous driving.

U-Net has two main components. The contractive path (or “encoder”) mirrors a traditional CNN, extracting features through convolutions, ReLU activations, and max-pooling. The expansive path (or “decoder”) uses up-convolutions and concatenates features from the encoder to recover spatial resolution. It progressively refines the output before a final convolution produces a full-resolution segmented image. Its symmetric, “U”-shaped architecture connects both paths (see Figure 1a).

Figure 1.

Schematics of U-Net and U-Net attention architectures: (a) example of U-Net architecture and (b) example of attention architecture.

3.1.2. U-Net Attention

A key feature in image processing networks is the ability to focus on important objects while ignoring less relevant areas. Oktay et al. (2018) enhances the classic U-Net by introducing an Attention Gate (AG) module.

An AG module selects relevant features and discards irrelevant ones. In this architecture (see Figure 1b), each block in the expansive path includes an AG, filtering contractive path features before merging them with upsampled features. This process significantly improves segmentation performance without excessive computational cost, increasing the model’s sensitivity and precision.

Among attention mechanisms, additive attention is preferred for image processing due to its higher accuracy (Schlemper et al. 2019).

3.1.3. DeepLabV3+

Before introducing DeepLabV3+, it is essential to discuss atrous convolution or dilated convolution (Chen et al. 2017). In standard convolution, a fixed-size filter moves over an input image, multiplying filter values with corresponding image values to produce a single output.

In dilated convolution, the filter expands by inserting spaces between its values. The dilation rate $r$ determines the number of inserted pixels and is a tunable hyperparameter. When $r=1$ , dilated convolution is identical to standard convolution (see Figure 2a). A key advantage of dilated convolutions is that they retain more input feature details across layers without increasing the amount of trainable parameters. This expands the receptive field without adding complexity, improving model performance. For each location $i$ in the output feature map $y$ , a dilated convolution is applied to input $x$ using filter $w$ , where $r$ defines the inserted pixel gaps. This is equivalent to convolving $x$ with filters augmented by inserting $r-1$ zeros between consecutive filter values along each spatial dimension. Adjusting $r$ dynamically modifies the filter’s field of view.

Figure 2.

Schematics of atrous convolution and DeepLabV3+ architecture: (a) example of atrous convolution and (b) example of DeepLabV3+ architecture.

The DeepLabV3+ architecture consists of an encoder and a decoder, similar to U-Net. Its global structure is shown in Figure 2b. The encoder includes a backbone network, $\mathrm{ResNet}50$ in this case, for feature extraction (He et al. 2015), along with an Atrous Spatial Pyramid Pooling (ASPP) module. The decoder refines the output feature map by concatenating it with low-level features and applying bilinear interpolation to restore the input image size, producing the final segmentation.

The ASPP module employs parallel branches with atrous convolutions at different dilation rates to capture multi-scale information. As seen in DeepLabV3+, ASPP includes a $1\times 1$ convolution, three $3\times 3$ atrous convolutions with dilation rates $r=6,12,18$ , and a pooling branch. The outputs are concatenated and compressed via a final $1\times 1$ convolution for feature integration.

3.2. Loss Functions

This section defines the loss functions considered as objectives to be optimized during network training. First, the $\mathrm{DICE}$ loss function is introduced, which quantifies the overlap between predictions and ground truth values. Next, the $\mathrm{Tversky}$ loss and $\mathrm{FocalTversky}$ loss functions are presented; these incorporate parameters to penalize false positives (FP) and false negatives (FN), especially in highly imbalanced scenarios. Finally, the $\mathrm{Ordinal}$ Loss function is introduced, which penalizes the distance between classes, in this case, the five poverty levels; so that, for instance, if the true poverty level is “1” and the model predicts “5,” a higher penalty is applied in contrast to the case where the model predicts “2” instead.

Let us consider ${\mathbf{y}}_c,{\hat{\mathbf{y}}}_c$ , where each element of them takes values from the set $\{0,1\}$ . The first one denotes the ground truth of classification for all the $N$ sampling units regarding class $c$ , while the second one denotes the predictions for that same class. TP, FP, and FN represent the counts of true positives, false positives, and false negatives, respectively.

3.2.1. DICE Coefficient

A key feature of the DICE loss function is its ability to quantify the similarity or discrepancy between predictions and ground truth values, while maintaining a direct connection with the proportion of true positives, false positives, and false negatives. The DICE coefficient for class $c$ is defined as

$${\mathrm{DICE}}_c:=\mathrm{DICE}({\mathbf{y}}_c,{\hat{\mathbf{y}}}_c)=\frac{2·{\mathrm{TP}}_c}{2·{\mathrm{TP}}_c+{\mathrm{FP}}_c+{\mathrm{FN}}_c}.$$

(8)

A common approach to reduce the effects of class imbalance is to introduce a weight $w_c$ for each class $c$ , which is inversely proportional to the frequency of that level (Sudre et al. 2017). In that line of ideas, the overall DICE can be computed as a weighted average: $\mathrm{DICE}=\sum _c{w}_c·{\mathrm{DICE}}_c$ .

One limitation of the DICE loss function is that it weighs false positives (FP) and false negatives (FN) equally. In practice, FN counts should be weighted more heavily than FP in order to improve the recall rate (Abraham and Khan 2019).

3.2.2. Tversky Loss

Tversky loss adds two parameters, $\alpha >0$ and $\beta >0$ , with $\alpha +\beta =1$ ; to allow for different weighting of FP and FN. If $\alpha >\beta$ , more importance is put on the false negative rate. This is particularly useful in datasets with marked imbalance in terms of the class frequencies, allowing for a better control on how the model handles the two types of classification errors (Salehi et al. 2017). For class $c$ , it is defined as

$$T_{\alpha ,\beta ,c}:=T_{\alpha ,\beta }({\mathbf{y}}_c,{\hat{\mathbf{y}}}_c)=\frac{{\mathrm{TP}}_c}{{\mathrm{TP}}_c+\alpha {\mathrm{FP}}_c+\beta {\mathrm{FN}}_c}.$$

(9)

If $\alpha =\beta =0.5$ , ${\mathrm{DICE}}_c$ metric is retrieved. In an analogous way, the overall Tversky loss can be computed as a weighted average: $T_{\alpha ,\beta }=\sum _c{w}_c·T_{\alpha ,\beta ,c}$ .

3.2.3. Focal Tversky Loss

The focal Tversky loss is motivated as an extension of the Tversky loss, improving its adaptability through an additional parameter $\gamma \in [1,3]$ . Salehi et al. (2017) state that its main advantage is the additional flexibility to adjust the importance of false positives and false negatives even more. By increasing $\gamma$ , this loss function focuses its attention on the most challenging classification cases. It is defined as

$${\mathrm{FT}}_{\alpha ,\beta ,\gamma ,c}:={\mathrm{FT}}_{\alpha ,\beta ,\gamma }({\mathbf{y}}_c,{\hat{\mathbf{y}}}_c)={\left(\frac{{\mathrm{TP}}_c}{{\mathrm{TP}}_c+\alpha {\mathrm{FP}}_c+\beta {\mathrm{FN}}_c}\right)}^{\gamma }.$$

(10)

In Abraham and Khan (2019), some characteristics and important considerations about this loss function are discussed. When $\gamma >1$ , the loss function focuses on penalizing more strictly those less precise predictions. In that same work, they suggest to use the optimized value of $\gamma =\frac{4}{3}$ , for that reason, that is the value that was used throughout this work. In an analogous way, the overall focal Tversky loss can be computed as a weighted average: ${\mathrm{FT}}_{\alpha ,\beta ,\gamma }=\sum _c{w}_c·T_{\alpha ,\beta ,\gamma ,c}$ .

3.2.4. Ordinal Loss

The ordinal loss function is utilized to quantify the difference between the predicted category and the real one, considering the distance between the two categories as a weighting factor within the overall loss function. The overall ordinal loss can be computed as

$$\mathrm{Ordinal}:=\frac{1}{N}\sum _{j=1}^N|t_j-{\hat{t}}_j|,$$

(11)

where $t_j$ is the observed label for the $j-$ th training unit, ${\hat{t}}_j$ is the predicted category, $|t_j-{\hat{t}}_j|$ corresponds to the amount of categories between the two in the natural order of the variable and $n$ is the size of the training set.

4. Material and Methods: Data Sources

This chapter describes the data sources. The main input for modeling consists of satellite images for Bogotá (The Capital District) and nearby municipalities, based on which a representative dataset is built. These images must meet two key criteria: (1) high resolution, since they must capture fine details such as rooftops and green areas, providing rich information; and (2) availability and recency, in the sense that they must be recent and cover the area of interest to accurately reflect the current situation.

The images used in this study were acquired from PlanetScope (PS) sensors, which are part of the largest satellite fleet in history with 130 satellites capturing daily images of Earth’s surface (Ghuffar 2018). The Geographic Institute for Colombia (Instituto Geográfico Agustín Codazzi, IGAC, for its acronym in Spanish) provided approximately 200 GB of PS data for academic and research purposes, conforming a large database covering the entire country from 2016 to 2019. Among PS satellites, some capture up to eight spectral bands, mainly for precision agriculture. However, this study focuses on four-band images, including red, green, blue (RGB), and near-infrared (NIR). PS images undergo pre-processing routines to ensure consistency, minimizing cloud interference and satellite orientation effects. Furthermore, each pixel has a resolution of 3 m (covering 9 m² in area), which makes them ideal for monitoring vegetation and urban structures in high detail.

On the other hand, data from the latest National Population and Housing Census (CNPV2018) was used to assign labels, based on the reported multidimensional poverty index, to the households covered by the satellite images (Departamento Administrativo Nacional de Estadística (DANE) 2018).

5. Material and Methods: Methodology

The selection process of the images available involved defining a rectangle encompassing Bogotá and surrounding municipalities (e.g., Mosquera, Soacha, La Calera), resulting in forty-seven images, each covering 120 to 210  $\mathrm{k}{\mathrm{m}}^2$ . Next, using the MPI data from CNPV2018 at the block level, blocks were classified into six levels (0–5). Each level has a specific meaning: Level 0 includes blocks lacking MPI data or where no residents were found during the census (not used in this study); Level 1 corresponds to an MPI between 1% and 10%; Level 2 ranges from 11% to 20%; Level 3 covers 21% to 33%; Level 4 spans 34% to 56%; and Level 5 includes blocks with an MPI above 57%. In conclusion, the higher the poverty level, the higher the MPI index, which means that the household is deprived of more basic needs.

Now, using as input the forty-seven images and the raster file containing the blocks classified into the five previously defined poverty levels (excluding level 0), the dataset required for the proposed methodology was created. This process involved creating “chips,” which are trimmed sections of the original images with a specific size. To achieve this, a function was applied to each image, taking the image, the raster file, and the parameters for size and step size as inputs. The choice of size and step size parameters determines the final dimensions of the chips and controls the level of overlap between them. For each of the forty-seven images, an “extent” was generated, representing the rectangle formed by the four corners of the image. In this work, the size parameter was set to 320 and the step-size to 320, meaning there is no overlap between adjacent chips. These values were carefully selected to ensure that each chip covers an appropriate area: neither too large to avoid including too many blocks, nor too small to avoid blurriness in the resulting chips. Consequently, each chip covers circa 0.92  ${\mathrm{m}}^2$ .

Out of the forty-seven initial images, a total of 18,222 chips were generated. However, 6,384 of these contained no useful information, as they corresponded to empty or blank areas in the images. These chips were excluded. Among the remaining 11,383 chips, a spatial intersection was performed to determine the amount of blocks contained in each chip. This step eliminated an additional 11,219 chips that did not contain relevant MPI-related information. An additional criterion was applied: ensuring that each chip contained blocks from all five poverty levels. This step aimed to prevent biases in the models due to missing categories. In conclusion, after removing the chips with blank spaces, those lacking MPI information and those that did not contain all five poverty levels, a total of 552 valid chips remained.

To conclude, Figure 3 presents the final results of the chip and mask generation process. The first row displays three examples of chips corresponding to Bogotá, while the second row shows the city blocks classified into the five respective poverty levels.

Figure 3.

Example of chip extraction (top row) and classification of blocks according to poverty levels (bottom row): (a) high poverty level, (b) low poverty level, (c) heterogeneous or no defined pattern of poverty level, (d) mostly high MPI blocks, (e) mostly low MPI blocks, and (f) heterogeneous MPI blocks.

The general strategy for fitting each of the above-mentioned models is shown in Figure 4. Steps 1 to 3 in that figure correspond to the pre-processing stages and were already described in this section. In step 4, the resulting 552 chips were split into training (402 chips, 72.8%), validation (93 chips, 16.8%), and test sets (57 chips, 10.4%) using a fixed seed for reproducibility. The two modeling approaches described in Sections 2 and 3 were then fitted to the data at hand. In the case of the CNN model, hyperparameters related to the problem under study, such as the corresponding loss function, the variants for the different weighting mechanisms for the categories, and the architectures (U-Net, U-Net Attention, DeepLabV3+) were exhaustively tried using both the training and validation sets. Other more common hyperparameters were taken from previous successful case studies, as described in the following subsection. The optimal configuration was then used to predict over the testing set. On the other hand, for the GTR model, each training-set chip required individual fitting of a model and predictions for the testing portion were issued by means of a weighted aggregation of the predictions issued by the neighboring chip models in the training set (around 8% of the chips in training). The weights were defined to be inversely proportional to the distance between the training-set chips and the target chip in the testing portion. Finally, in step 5, the two modeling approaches were compared in terms of their performance metrics when predicting the poverty level on the testing portion.

Figure 4.

Methodology summary.

5.1. CNN Model Configuration

This section outlines the CNN hyperparameters used for optimal model calibration, adopting values from Benedetti et al. (2022) to ensure reproducibility rather than conducting exhaustive hyperparameter optimization due to computational costs. Image augmentation techniques (cropping, zoom, x/y-axis rotation) were applied to address limited training data. The model used a learning rate of 0.001, trained for forty epochs with a batch size of sixteen. Each epoch processed thirty-two steps to ensure full data coverage, and early stopping was implemented with a patience of ten epochs to prevent overfitting. These parameters balanced computational efficiency with model performance while maintaining comparability with prior research.

5.1.1. Weighting for CNN Training

Different weighting mechanisms for the categories in each loss function were considered. These weighting strategies were essential due to the inherent class imbalance in the dataset. As described in this section, 552 chips were selected for model fitting, ensuring that each poverty level was represented. However, it was not possible to guarantee equal proportions of each class within every chip. In practice, poverty level 1 consistently appeared more frequently, followed by level 2, and so forth, with level 5 being the least frequent. This asymmetry reflects real-world class distributions and underscores the need to correct for potential learning biases during training. To address this, a weight vector $\mathbf{w}={(w_1,w_2,w_3,w_4,w_5)}^T\in {\mathbb{R}}^5$ was introduced to balance the contribution of each class to the loss function. This vector plays a critical role in ensuring that underrepresented categories, particularly those with fewer training units, are given higher importance proportionally. As a result, the model was encouraged to consider minority classes more attentively, leading to better generalization and performance on unseen data.

Weighting alternatives can be divided into two classes: (1) those in which the weights do not sum to 1, and (2) those in which the weights do sum to 1. Additionally, three weighting strategies are introduced: proportional, equal, and unweighted. These strategies were defined to assess the impact of different weighting schemes on poverty level analysis.

The proportional approach assigns weights according to the relative frequency of each poverty level, giving greater importance to less frequent levels. This allocation is based on the pixel-level distribution corresponding to each poverty level across the 552 chips analyzed. Second, the equal approach assigns the same weight to each poverty level, with the exact values depending on the specific method applied. Finally, to evaluate the impact of using weights, a third category is included in which no weighting is applied. In this case, the assigned weight is always equal to 1, serving as a baseline for comparing results across the other approaches.

5.2. Spatially-Weighted GTR (SWGTR) Model Configuration

This section describes the methodology for obtaining GTR model predictions on test-set chips using a spatially-weighted adaptation that incorporates the GTR model described in Section 2.

5.2.1. Distance-Weighted Predictions

Predictions were weighted using inverse distance weighting (IDW) following Tobler’s Law, which assigns higher weights to geographically proximate observations (Wikle et al. 2019, 79). Considering $u$ trained models in the vicinity of the target, the IDW predictor $t\left(s_0\right)$ for each $320\times 320$ pixel in the testing portion is defined as a weighted majority vote mechanism, that is, its prediction corresponds to

$$\hat{t}\left(s_0\right)=\arg \underset{c\in \{1,2\dots ,C\}}{\max }{\delta }_c\left(\hat{\mathit{\nu }}\left(s_0\right)\right),$$

(12)

where

$$\begin{array}{c}\multicolumn{1}{c}{\hat{\mathit{\nu }}\left(s_0\right)={\left({\hat{\mathit{\nu }}}_i\left(s_0\right)\right)}_{i=1,2,\dots ,5},}\\ \multicolumn{1}{c}{{\hat{\mathit{\nu }}}_i\left(s_0\right)=\sum _{j=1}^u{p}_j\left(s_0\right)·{\mathbb{I}}_{({\hat{t}}_j=i)}.}\end{array}$$

The weight terms, $p_j\left(s_0\right)$ , define the relative importance of the prediction issued by the model at location $s_j$ based on its distance to $s_0$ and can be computed as $p_j\left(s_0\right)=\frac{\widehat{p_j\left(s_0\right)}}{\sum _{l=1}^u\widehat{p_l\left(s_0\right)}}$ , where $\widehat{p_j\left(s_0\right)}=\frac{1}{d{(\left(s_j\right),\left(s_0\right))}^{\alpha }},$ and $d(\left(s_j\right),\left(s_0\right))$ is the “distance” between the training chip located at $s_j$ and the test chip centroid at $s_0$ , for all $j$ in the vicinity of it. The exponent $\alpha$ was employed as a smoothing parameter. While this parameter is typically determined through cross-validation techniques, in this study it was fixed at a value of 2; $\mathbb{I}$ denotes the indicator function that takes the value of $1$ when the prediction issued by model $j$ , ${\hat{t}}_j$ , matches category $i$ for MPI. In summary, the vector $\hat{\nu }\left(s_0\right)$ gathers the weighted votes in support of each of the MPI categories, $i=1,2,\dots ,5$ ; and $\hat{t}\left(s_0\right)$ corresponds to the category with the highest weighted voting score. For illustrative purposes, Figure 5 exemplifies the voting scheme. To prevent spatial oversmoothing (where all pixels might converge to a single poverty level), we implemented an 8% proximity threshold, using only the $u=32$ nearest training chips for each prediction. The thirty-two nearest chips were chosen based on the planar Euclidean distance between chip centroids. However, only the 402 chips belonging to the training portion were used as potential neighbors, avoiding the risk of spillover.

Figure 5.

Weighted-distance voting mechanism for prediction and cumulative sum curve for model evaluation: (a) schematics of weighted inverse-distance-based voting mechanism and (b) cumulative sums curve (CS) and the ABC indicator.

5.3. Classification Performance Metrics for Multi-Category Models

We evaluated model performance using established metrics: Overall Accuracy, Kappa Index, and Average Accuracy from Liu et al. (2019); Cumulative Sum (CS) curves adapted from Niu et al. (2016); and the Area between Curves (ABC) metric, proposed by the authors, as an adaptation of the multi-class ABC measure in Hand and Till (2001) computed for the CS curve, which accounts for the ordinal nature of the variable of interest. Although the ordinal variable does not have a quantitative nature, the literature also considers Mean Absolute Error (MAE) as a measure to assess the performance of ordinal classification (Niu et al. 2016). Its value can be interpreted as the average number of classes between the predicted class and the actual one.

Assuming that the confusion matrix of the classification problem can be denoted as $\mathbf{A}:={\left(a_{\mathrm{kl}}\right)}_{k,l}$ , where $a_{\mathrm{kl}}$ corresponds to the total of pixels belonging in reality to category $k$ and classified in category $l$ ; and $N_T$ denotes the total of pixels in the testing portion, each of the metrics is discussed in detail.

• Overall Accuracy: it corresponds to the proportion of correctly classified pixels. It measures the classifier’s general prediction accuracy and can be computed as

$$\mathrm{OA}=\frac{{\sum }_{k=1}^5{a}_{\mathrm{kk}}}{N_T}.$$

• Kappa Index: the Kappa index is a metric that assesses the agreement between a classifier’s predictions and the actual labels, accounting for the likelihood of random classifications. It measures the proportion of agreement between the classifier and the true labels while adjusting for expected chance agreement. The Kappa index is calculated using the following formula

$$\mathrm{KI}=\frac{\mathrm{OA}-\mathrm{Probability}\mathrm{of}\mathrm{Random}\mathrm{Agreement}}{1-\mathrm{Probability}\mathrm{of}\mathrm{Random}\mathrm{Agreement}}.$$

• Average Accuracy: the average accuracy is a metric that provides a more detailed view of a classifier’s performance by considering the accuracy for each individual class. Instead of evaluating only the overall accuracy, this metric calculates the accuracy for each class separately and then averages these values. This is especially useful when classes are not equally represented in the dataset. The formula for calculating the average accuracy is expressed as

$$\mathrm{AA}=\frac{1}{5}\sum _{k=1}^5\left(\frac{a_{\mathrm{kk}}}{\sum _{l=1}^5{a}_{\mathrm{kl}}}\right).$$

• Cumulative Sums Curve (CS): this curve maps the distance between observed and predicted labels and the cumulative frequency for that distance value. For a null distance, it adds the frequencies of exact matches, which is equivalent to OA. For a one-unit distance, it adds OA with the relative frequency of classification cases that are at most one category apart from the true label, that is, it considers the relative contribution of all $a_{\mathrm{kl}}$ such that $|k-l|\le 1$ . This process is carried out until the maximum distance between categories is reached. In general, for a given value of distance $\delta$ , $\mathrm{CS}\left(\delta \right)=\sum _{k,l:|k-l|\le \delta }\left(\frac{a_{\mathrm{kl}}}{N_T}\right)$ . The results are plotted with the x-axis showing the distance between prediction and ground truth $\delta$ , and the y-axis showing the cumulative percentage up to that distance, $\mathrm{CS}\left(\delta \right)$ .

• Figure 5b shows an example where the cumulative sum reaches 100% with at most a two-level difference. Therefore, model selection using this metric favors higher CS values at smaller prediction-error distances. This metric is particularly useful because it accounts for the ordinal nature of the variable of interest instead of penalizing all mismatches (the elements outside of the diagonal of matrix $\mathbf{A}$ ).

• Area Between Curves (ABC): the ABC metric evaluates classification performance by measuring the area between the model’s CS curve and the horizontal axis. A higher ABC value indicates better performance, reflecting the model’s ability to distinguish poverty levels. It is calculated via numerical integration of the CS curve. Due to its design, ABC offers a comprehensive assessment of model performance. Figure 5 illustrates the shaded area of interest, with a value of 95%. This score is scaled between 0 and 1, since the unscaled ABC can range from 0 to 4.

• Mean Absolute Error (MAE): although ordinal labels do not possess all mathematical properties of quantitative variables, MAE is commonly reported in the literature as a measure of prediction closeness in ordinal response models (Niu et al. 2016). It is computed as

$$\mathrm{MAE}=\frac{1}{N_T}\sum _{l=1}^{N_T}|y_l-{\hat{y}}_l|,$$

•  where $N_T$ is the total number of observations in the testing portion, $y_l$ is the true label for pixel $l$ and ${\hat{y}}_l$ is the predicted poverty level for pixel $l$ .

6. Analysis, Comparison, and Discussion of Results

6.1. Comparison Between CNN Architectures

To begin this section, a summary is provided detailing the results of the CNN models applied to the chips of Bogotá, considering the different combinations of loss function (DICE, Tversky, Focal Tversky, and Ordinal), network architecture (UNET, Attention, and DeepLabV3), and weights (scaled and unscaled).

It is important to note that the CNN configurations explored in this study were intentionally constrained in terms of hyperparameter tuning, architectural complexity, and input design. This choice was made to maintain computational feasibility and reproducibility, and to provide an initial, controlled comparison with the alternative modeling strategy. As such, the results presented here should be interpreted as exploratory rather than as a fully optimized benchmark of CNN performance.

The criterion for selecting the models with the best metrics involves identifying the combination of: (a) loss function, (b) architecture, and (c) weights that display the highest values in the DICE indicator. This selection is intended to provide a set of competitive baseline models for comparison, rather than to establish an exhaustive or globally optimal CNN specification.

This approach allows for the identification of three specific models, which will be used for comparison with the results of the GTR model. The selected models are as follows: Model 1—CNN, which trains under the Focal Tversky Loss function and uses the UNET architecture and does not use weights. Model 2—CNN, which trains the DICE loss function, the Attention architecture and does not use weights; and Model 3—CNN, which applies the Focal Tversky Loss function, the DeepLabV3+ architecture, and uses unscaled proportional weights, specifically $(0.63,0.72,0.82,0.88,0.92)$ , for each category.

Additionally, the analysis showed that including scaled weights (i.e., those summing to one) did not provide significant benefits compared to unscaled weights. This was evident as all tested combinations of architecture, loss function, and weights resulted in DICE values under 0.5.

6.2. Comparison of CNN and SWGTR Strategies

This subsection presents the resulting predictions from both approaches, the CNN and the SWGTR model, applied to the test set chips. This involved computing confusion matrices for each selected model to obtain a comparison table summarizing their performance using the metrics defined in Subsection 5.3: Overall Accuracy (OA), Average Accuracy (AA), Kappa Index (KI), Area Between Curves (ABC), and Mean Absolute Error (MAE). Results are first shown for a single chip to visualize individual model predictions, followed by aggregated results across all test chips.

It is important to highlight that the training, validation, and testing splits were not constructed using spatial partitioning. As a result, neighboring chips may be present across these subsets. Given that the SWGTR methodology explicitly incorporates spatial proximity through weighted aggregation, this setup may favor its predictive performance relative to aspatial approaches such as CNNs.

6.2.1. Comparison Using a Single Test Chip

Figure 6 shows the results for a test chip located at coordinates (4.580609, −74.193558), in the San Mateo neighborhood in the municipality of Soacha. It includes the original chip, displaying the spatial layout of buildings, roads, and green areas, the ground truth mask showing the poverty level per pixel, and the predictions from the three selected CNN models and the SWGTR model. In the ground truth mask (see Figure 6b), the chip is predominantly classified as poverty level 1. However, the right side, particularly in mountainous regions, shows higher poverty levels. This contrast is attributed to the challenging terrain, which hinders housing development and limits access to basic services, contributing to higher poverty in those areas. For illustrative purposes, an inverse-distance weighting benchmark was also considered to isolate the ability of the spatial neighbors to predict poverty levels without the mediation of a model, in contrast to SWGTR, which combines information from models trained on neighboring chips. The neighboring chips were chosen as the top 8% of closest chips within the training set, mimicking the spatial voting mechanism incorporated within the SWGTR methodology.

Figure 6.

Model results compared to the original chip and ground truth: (a) original chip, (b) ground truth, (c) Model 1—CNN, (d) Model 2—CNN, (e) Model 3—CNN, and (f) SWGTR model.

A visual inspection suggests that none of the models fully replicates the ground truth mask. Among the three CNNs, Models 1 and 2 show similar predictions for poverty level 1 but struggle with other levels. For instance, Model 3 predicts mostly levels 1 and 2, with only a few pixels labeled as level 3 or higher, failing to capture any of the true level 5 areas. In contrast, the SWGTR model produces a distribution of level 1 pixels (in green) that closely resembles the ground truth. However, it tends to underestimate higher poverty levels. Specifically, pixels labeled as level 5 in the reference mask are mostly classified as level 4. While the SWGTR does not perfectly match each level, it provides a closer approximation of the actual poverty distribution in this area. At the same time, it is worth noting that the CNN-based models exhibit stronger performance in specific metrics and localized structures, suggesting that their limitations may be partly associated with the constrained training setup rather than an inherent inability to capture the underlying patterns.

The IDW model contributes more heterogeneity into the predictions, but it does not appear to learn the systematic patterns required to assign accurate labels to each pixel.

Further assessing model performance, Table 1 summarizes evaluation metrics for the previously analyzed chip. It includes results for all four models, Model 1, 2, and 3 (CNNs), the SWGTR model, and the IDW benchmark; across the following metrics: OA, AA, KI, ABC, and MAE.

• Overall Accuracy (OA): Model 2 achieved the highest OA (0.723), followed by Model 1 (0.683), SWGTR (0.556), and Model 3 (0.166). Thus, Model 2 shows the best overall classification performance.

• Average Accuracy (AA): Model 3 had the highest per-class accuracy (0.226), followed by Model 2 (0.200), Model 1 (0.193), and SWGTR (0.176).

• Kappa Index (KI): SWGTR achieved the highest Kappa (0.341), followed by Model 1 (0.131), Model 2 (0.000), and Model 3 (−0.076), making SWGTR the strongest performer by this measure.

• Cumulative Sums Curve (CS): as shown in Figure 7, Model 2 performs best at zero distance (exact matches). At distances of 1 to 2 levels, all models behave similarly. However, at a distance of 3, SWGTR outperforms the rest.

• Area Between Curves (ABC): This metric reflects the area under each model’s curve, with higher values indicating better performance. Model 1 achieved the highest ABC (0.825), followed by SWGTR (0.820), Model 2 (0.802), and Model 3 (0.766). Thus, Model 1 leads in this metric, closely trailed by SWGTR.

• Mean Absolute Error (MAE): Model 1 also had the lowest MAE (0.860), indicating the smallest average error. It was followed by Model 2 (0.922), SWGTR (0.945), and Model 3 (1.362). This highlights Model 1 as the top performer in prediction accuracy.

Table 1.

Performance Metrics for the Four Selected Models and the Benchmark for the Chip Shown in Figure 6.

Metric	Model 1—CNN	Model 2—CNN	Model 3—CNN	SWGTR model	IDW model
OA	0.683	0.723	0.166	0.556	0.252
AA	0.193	0.200	0.226	0.176	0.185
KI	0.131	0.000	−0.076	0.341	0.12
ABC	0.825	0.802	0.766	0.820	0.691
MAE	0.860	0.922	1.362	0.945	1.614

Figure 7.

Cumulative sums for all the four models and the benchmark selected for the chip shown in Figure 6.

Taken together, these results indicate that different models excel under different evaluation criteria. For instance, Model 2 achieves the highest overall accuracy, while Model 1 performs best in terms of ABC and MAE, and SWGTR attains the highest Kappa index and competitive ordinal performance. Rather than establishing a definitive ranking, these findings highlight trade-offs across models depending on the metric of interest. In particular, metrics such as ABC and MAE, which account for the ordinal nature of the response variable, suggest that SWGTR provides competitive performance in capturing the relative ordering of poverty levels, as confirmed by the analysis of the entire testing set of chips in the next subsection. However, CNN-based models remain competitive under nominal accuracy measures.

While CNN models face challenges, particularly in capturing higher poverty levels under the current configuration, the SWGTR model, although not flawless, tends to provide a closer approximation to the spatial distribution of poverty in this specific setting. This observation should be interpreted in light of the methodological choices described above.

The inferior performance of the IDW model in both Table 1 and Figure 7 showcases that, although considering spatial information in predicting poverty is important, there is a need for a model that extracts the important features from the images, rather than merely weighting in the information of the spatial neighboring chips to create a prediction. However, as part of the future work, a more elaborate experimental setting that incorporates spatially structured data splitting is important to better disentangle the contributions of spatial weighting and the GTR model.

Consequently, this single-chip analysis should be interpreted as illustrative of model behavior rather than as a conclusive model selection exercise.

6.2.2. Comparison Using the Entire Testing Set

This section presents the findings based on the complete set of test chips. Accordingly, to facilitate interpretation, Figure 8 displays boxplots summarizing the performance of all models across the evaluation metrics considered in this study. The visualizations reveal the following:

• Overall Accuracy (OA): the SWGTR model exhibits the highest mean and median overall accuracy, followed by Model 1, Model 2, and Model 3. Thus, SWGTR exhibits higher overall accuracy under the current experimental setup.

• Average Accuracy (AA): all models show comparable performance; however, SWGTR displays a higher presence of outliers, indicating greater variability in class-wise precision.

• Kappa Index (KI): the models exhibit similar and generally low Kappa values, rendering this metric inconclusive in the current context.

• Area Between Curves (ABC): while all models present similar trends, SWGTR attains the highest median ABC value, followed closely by Model 1, indicating stronger ordinal consistency.

• Mean Absolute Error (MAE): SWGTR shows the lowest median MAE, suggesting that it provides a closer match between predicted and actual poverty levels under this metric.

Figure 8.

Boxplots with the results of all fifty-seven chips in the test set: (a) boxplots for metrics AA, ABC, KI, OA and (b) boxplots for metric MAE.

On the other hand, Table 2 presents the average results across all fifty-seven test chips for each metric and model. This summary indicates that, on average, the SWGTR model attains higher values in OA and ABC metrics and lower MAE compared to the selected CNN configurations. However, these differences should be interpreted in light of the experimental design, including the spatial structure of the data split and the constrained CNN setup.

Table 2.

Average Performance Metrics for the Four Selected Models and All the Chips.

Metric	Model 1—CNN	Model 2—CNN	Model 3—CNN	SWGTR model
OA	0.392	0.364	0.349	0.419
AA	0.204	0.200	0.216	0.203
KI	0.066	0.000	0.010	0.012
ABC	0.821	0.792	0.816	0.828
MAE	1.064	1.249	1.109	0.996

The findings above are broadly consistent with the single-chip analysis, where the CNN models show difficulties in capturing the full range of poverty levels under the current configuration, while the SWGTR model tends to provide a closer approximation to the observed distribution. It is important to highlight that the first three metrics (OA, AA, and KI) focus on nominal classification performance and do not account for the “closeness” between predicted categories. For this reason, greater importance should be placed on the ABC and MAE metrics.

It should be noted that the SWGTR model performs more accurately in central areas of Bogotá, where an abundant number of neighboring chips is available. In contrast, as the location moves toward the outskirts, such as in the municipalities of Bojacá and Madrid, the model tends to misclassify high poverty levels as lower ones. This decline in accuracy is largely attributed to the limited availability of neighboring contextual information, which restricts the model’s ability to interpolate effectively in more remote regions. These observations motivate further investigation into more balanced experimental designs, including spatial cross-validation schemes and more extensive CNN optimization, in order to more clearly disentangle the effects of spatial information and model capacity on predictive performance.

6.3. Potential Use for Survey Targeting and Uncertainty-Guided Sampling

One of the motivations for this study is the potential to complement traditional survey-based approaches to poverty measurement by providing spatially explicit information derived from satellite imagery. While the primary focus of this work is on predictive performance, the SWGTR framework also enables the construction of uncertainty measures that can be useful for guiding data collection efforts.

In particular, for each pixel and predicted poverty class $c$ , an uncertainty measure can be derived from the voting mechanism underlying SWGTR. Let ${\hat{\mathit{\nu }}}_c\left(s_0\right)$ denote the estimated probability that a pixel at location $s_0$ belongs to the chosen class $c$ . Each pixel is assigned to the category that maximizes this probability. This quantity can be interpreted as the estimated mean of a Bernoulli random variable, allowing us to define its corresponding standard deviation as

$$s_c\left(s_0\right)=\sqrt{{\hat{\mathit{\nu }}}_c\left(s_0\right)(1-{\hat{\mathit{\nu }}}_c\left(s_0\right))}.$$

This measure captures the level of classification uncertainty, with higher values indicating lower confidence in the assigned class.

Figure 9 illustrates an example of such an uncertainty map for a selected chip. Areas with higher uncertainty tend to correspond to regions where class boundaries are less well defined or where the model exhibits greater ambiguity in assigning poverty levels.

Figure 9.

Example of pixel-level uncertainty derived from the SWGTR voting mechanism.

From a practical standpoint, these uncertainty estimates could be integrated into survey design. For instance, a probability-proportional-to-size (PPS) sampling scheme could be adapted such that the inclusion probability of a spatial unit is proportional to an aggregate measure of prediction uncertainty. This would prioritize areas where the model is less certain, allowing survey-based data collection to focus on regions where additional information is most valuable.

We emphasize that this application is illustrative and not fully developed in the current study. A more rigorous integration of model-based uncertainty into survey design remains an important direction for future research.

7. Conclusions, Limitations, and Future Work

This study investigated the use of satellite imagery to estimate the Multidimensional Poverty Index (MPI) in Bogotá, Colombia; with the goal of complementing traditional survey-based methods such as the ENCV and CNPV. By generating poverty maps from remote sensing data, this approach aims to provide timely and spatially detailed insights that are otherwise limited by the infrequency of national surveys.

Two modeling strategies were explored: Convolutional Neural Networks (CNNs) and a modified regression approach called Spatially-Weighted Generalized Tensor Regression (SWGTR). Model performance was evaluated using both numerical metrics derived from confusion matrices, and visual tools, including CS curves and a hybrid metric (ABC) that quantifies the area under these curves. These two approaches were compared under a common experimental framework designed to provide an initial assessment of their relative strengths and limitations when applied to high-resolution satellite imagery.

The empirical results indicate that both modeling strategies present strengths and limitations depending on the evaluation criterion considered. Under the configurations explored in this study, the CNN-based approaches show difficulties in capturing the full range of poverty levels, particularly at higher categories, although they remain competitive under certain nominal metrics such as overall accuracy.

The SWGTR approach, in turn, demonstrates competitive performance across several metrics, particularly those that account for the ordinal nature of the response variable, such as ABC and MAE. In this sense, SWGTR appears to provide a useful framework for capturing relative differences in poverty levels across space. At the same time, these findings should be interpreted with caution. The CNN models were implemented under constrained configurations in terms of hyperparameter tuning and computational resources, and the evaluation design did not enforce spatial separation between training, validation, and testing samples. As a result, part of the observed performance differences may reflect these design choices rather than intrinsic differences in model capability.

In general, the results highlight the importance of incorporating evaluation metrics that explicitly account for the ordinal structure of poverty levels, as they provide a more nuanced perspective on predictive performance beyond standard nominal accuracy measures.

From a data construction perspective, an important limitation of this study arises from the selection of chips that contain all poverty levels. While this criterion ensures that each training instance provides information across the full range of categories, it might also induce selection bias toward highly heterogeneous areas, which are more likely to correspond to dense urban environments. Consequently, the results presented here may not generalize to more homogeneous regions, such as rural or peri-urban areas, where poverty patterns exhibit less spatial variability.

In addition, the absence of a spatially explicit validation strategy implies that the SWGTR methodology explicitly might have a particular advantage as it leverages spatial proximity.

Finally, the CNN models were evaluated under simplified configurations, without extensive hyperparameter optimization, pretraining on large-scale remote sensing datasets, or advanced data augmentation strategies. While this choice was made to maintain computational feasibility and comparability, it opens a road to fully exploit CNNs’ representational capacity in this context.

Future work should address the limitations identified above through complementary directions. First, expanding the data construction process to include a broader range of spatial contexts, including more homogeneous chips, would allow for a more comprehensive assessment of model performance across different types of landscapes.

Second, implementing spatial cross-validation or geographically structured train-test splits would provide a more thorough evaluation of the contribution of spatial information, particularly for models such as SWGTR that explicitly rely on proximity-based mechanisms.

Third, further research should explore more advanced CNN configurations, including transfer learning from large-scale remote sensing datasets, more extensive hyperparameter tuning, and class-imbalanced training strategies, in order to more fully assess their potential in this application.

Finally, an important applied line of work would be to explore how model outputs could explicitly and effectively support the sampling design of poverty surveys conducted by National Statistical Offices (DANE in the Colombian context). In particular, areas where the models exhibit higher predictive uncertainty or systematic errors could be prioritized for targeted data collection, enabling a more efficient allocation of survey resources.