Prediction of Socioeconomic Indicators in Vale do Ribeira using Deep Learning and Satellite Imagery

ISABELLA DE MELO SOUSA PREDICTION OF SOCIOECONOMIC INDICATORS IN VALE DO RIBEIRA USING DEEP LEARNING AND SATELLITE IMAGERY S˜ao Paulo 2022 ISABELLA DE MELO SOUSA PREDICTION OF SOCIOECONOMIC INDICATORS IN VALE DO RIBEIRA USING DEEP LEARNING AND SATELLITE IMAGERY Work presented to the Escola Polit´ecnica da Universidade de S˜ao Paulo in order to obtain the Bachelor’s Degree in Computer Engineering. S˜ao Paulo 2022 ISABELLA DE MELO SOUSA PREDICTION OF SOCIOECONOMIC INDICATORS IN VALE DO RIBEIRA USING DEEP LEARNING AND SATELLITE IMAGERY Work presented to the Escola Polit´ecnica da Universidade de S˜ao Paulo in order to obtain the Bachelor’s Degree in Computer Engineering. Advisor: Prof. Dr. Pedro Luiz Pizzigatti Corrˆea Co-Advisor: PhD Marina Jeaneth Machicao Justo S˜ao Paulo 2022 ACKNOWLEDGMENTS I would like to thank my family for their unconditional support throughout my graduation. As well as my dear friends for their encouragement, in special Gabriel Barbutti de Lima Baker, Megumi Tsuru and Vinicius Akira Imaizumi, for their support and insightful discussions about the project. I also would like to express my sincere gratitude to my supervisors, Prof. Dr. Pedro Luiz Pizzigatti Corrˆea and PhD Marina Jeaneth Machicao Justo, for guiding me in this difficult and rewarding journey that was developing this project. I thank them for all the help, attention and above all, their patience. Finally, I am also grateful to the PARSEC group and all its researchers for their support in conducting this project. RESUMO Medidas-chave dos indicadores socioeconˆomicos s˜ao essenciais para a tomada de decis˜oes pol´ıticas informadas, mas devido aos altos custos e dificuldades operacionais dos esfor¸cos tradicionais de coleta de dados, a obten¸c˜ao de dados socioeconˆomicos confi´aveis continua sendo um desafio, principalmente nos pa´ıses em desenvolvimento. Este trabalho apresenta uma metodologia de aprendizado profundo para estimar indicadores socioeconˆomicos utilizando imagens de sat´elite. O modelo de rede neural desenvolvido foi treinado na regi˜ao brasileira de S˜ao Paulo e Paran´a com o objetivo de analisar o indicador socioeconˆomico de renda na regi˜ao do Vale do Ribeira. O modelo resultou em um R-quadrado de 0,4016 e teve um desempenho significativamente melhor do que o modelo treinado apenas com as bandas RGB. Palavras-Chave – aprendizagem profunda, imagens de sat´elite, indicadores socioeconˆomicos. ABSTRACT Key measures of socioeconomic indicators are essential for making informed policy decisions, but due to the high costs and operational difficulties of traditional data collection efforts, obtaining reliable socioeconomic data remains a challenge, particularly in developing countries. This work presents a deep learning methodology to estimate socioeconomic indicators using satellite imagery. The neural network model developed was trained at the Brazilian region of S˜ao Paulo and Paran´a with the goal of analyzing the socioeconomic indicator of income in the Vale do Ribeira region. The model yielded a R-squared of 0.4016 and performed significantly better than the model trained only on RGB bands. Keywords – deep learning, satellite imagery, socioeconomic indicators LIST OF FIGURES 1 Multilayer neural network topology. . . . . . . . . . . . . . . . . . . . . . . 17 2 High-level general CNN architecture. . . . . . . . . . . . . . . . . . . . . . 18 3 A schematic display of 5-fold CV. A set of n observations is randomly split into five non-overlapping groups. Each of these fifths acts as a validation set (shown in beige), and the remainder as a training set (shown in blue). The test error is estimated by averaging the five resulting MSE estimates. . 22 4 Separation of the sample into training, validation and test subsets. . . . . . 23 5 Vale do Ribeira and its municipalities. . . . . . . . . . . . . . . . . . . . . 30 6 Workflowdiagram. ............................... 34 7 Boxplot of census areas by urban or rural types. . . . . . . . . . . . . . . . 35 8 Cluster boundaries of the municipality of Apiai, SP. . . . . . . . . . . . . . 36 9 Grid zoomed to the Vale do Ribeira region. . . . . . . . . . . . . . . . . . . 38 10 Bands plotted for block number 2916, centered at the coordinates (26.314129◦ S, 51.276913◦W). The color map refers to the normalized pixel values. . . 40 11 ResNet18 with preactivation adapted to accept as input multi-band satellite imagery with C channels and to do a regression instead of a classification. . 42 12 Training and validation curves for MS models. The red lines represent the checkpoints where the model obtained the lowest MSE. . . . . . . . . . . . 44 13 Training and validation curves for NL models. The red lines represent the checkpoints where the model obtained the lowest MSE. . . . . . . . . . . . 45 14 Example of leave-one-group-out cross-validation for the concatenated model MS+NL. In this case, lowest MSE = 0.001 and best alpha = 16 . . . . . . 46 15 Regression plot for the combined MS+NL model. The dotted line corresponds to the line of best fit. . . . . . . . . . . . . . . . . . . . . . . . . . . 48 16 Finalmetricsresults............................... 48 17 Correlation results for the income predictions of all models a) comparison of r²metric and b) comparison of rank metric . . . . . . . . . . . . . . . . 49 18 Regression plot after calculating the income for each municipality. The dotted line corresponds to the line of best fit. . . . . . . . . . . . . . . . . 50 19 Heatmaps for a) Real HDI Income indicator and b) Predicted HDI Income indicator. .................................... 51 20 Differences between the real and predicted labels for the municipalities of SPandPR. ................................... 51 21 Regression plot for the municipalities of Vale do Ribeira. The dotted line corresponds to the line of best fit. . . . . . . . . . . . . . . . . . . . . . . . 52 22 Heatmaps for a) Real HDI Income indicator for VR and b) Predicted HDI IncomeindicatorforVR............................. 53 23 Differences between the real and predicted labels for the municipalities from VR. ....................................... 53 14 Furthermore, this research is being developed as part of the PARSEC project [2] and it aims to share the results with the group, in which one of its goals is to promote adequate monitoring of the environmental situation in the various biomes in Brazil and to analyze the socioeconomic influence of protected areas on the cities that surround them. PART II THEORETICAL FOUNDATION AND STATE-OF-THE-ART 16 2 THEORETICAL FOUNDATION AND STATE-OF-THE-ART 2.1 Deep Learning Artificial Neural Networks are a computational model that shares some properties with the animal brain in which many simple units are working in parallel with no centralized control unit. The weights between the units are the primary means of long-term information storage in neural networks. Updating the weights is the primary way the neural network learns new information. (PATTERSON; GIBSON, 2017) [3]. The behavior of neural networks is shaped by its network architecture, which can be essentially defined by the following: number of neurons, number of layers and types of connections between layers. [3] The most well-known and simplest-to-understand neural network is the feedforward multilayer neural network. It has an input layer, one or many hidden layers, and a single output layer. Each layer can have a different number of neurons and each layer is fully connected to the adjacent layer. The connections between the neurons in the layers form an acyclic graph, as illustrated in Figure 1. [3] 17 Figure 1: Multilayer neural network topology. Source: Patterson and Gibson, (2017) [3]. Deep Learning refers to the use of Artificial Neural Networks with three or more layers [4]. Note that additional hidden layers can help optimize and refine the accuracy of the model, but on the other hand they increase its complexity and computational cost. 2.1.1 Convolutional Neural Network (CNN) A Convolutional Neural Network (CNN) is an a Deep Learning technique in which the goal is to learn higher-order features in the data via convolutions. According to Patterson and Gibson (2017) [3], although there are several variations in the architecture of a convolutional neural network, in general it is based on the pattern of layers shown in Figure 2, consisting essentially of three major groups: 1. Input layer: accepts three-dimensional input generally in the form spatially of the size (width ×height) of the image and has a depth representing the color channels (generally three for RGB color channels); 2. Feature-extraction (learning) layers: find a number of features in the images and progressively construct higher-order features; 3. Classification layers: consists in one or more fully connected layers to take the higher-order features and produce class probabilities or scores. 18 Figure 2: High-level general CNN architecture. Source: Patterson and Gibson, (2017) [3]. 2.2 Regression 2.2.1 Linear Regression Linear regression is a linear approach for modelling the relationship between a variable of interest and one or more explanatory variables (also known as dependent and independent variables). If there is just one explanatory variable, the method is known as simple linear regression; if there are more, it is known as multiple linear regression. The multiple linear regression model has the form: yi= ˆyi+ϵi(2.1) And ˆyi=β0+ p X j=1 β0xij, i = 1,2, ..., n (2.2) Where: •yirepresents the ith observed response value; •ˆyirepresents the ith response value that is predicted; •β0, β1, ..., βprepresents the coefficients to be estimated (also known as weights); 19 •xij represents the independent variable; •n represents the size of the sample; •p represents the number of independent variables. To estimate the regression coefficients, a common approach is the linear least squares function, that consists in minimizing the residual sum of squares (RSS), that can be defined as: RSS = n X i=1 ϵi.(2.3) Where ϵi=yi−ˆyirepresents the ith residual, which is the difference between ith observed response value and the ith response value that is predicted by the linear model. In this methodology, the model is not penalized for its choice of weights. This means that for some features the model may place weights that are too large, leading to overfitting. 2.2.2 Ridge Regression The Ridge Regression solves a regression model where the loss function is the linear least squares function and regularization is given by the l2-norm. It avoids overfitting because it shrinks the regression coefficients by imposing a penalty on their size and it is also very useful to deal with multicolinearity between the independent variables. The ridge coefficients minimize a penalized residual sum of squares, given by: RSS +λ p X j=1 β2 j(2.4) Here λ≥0 is a complexity parameter that controls the amount of shrinkage: the larger the value of λ, the greater the amount of shrinkage. [5]. 2.3 Evaluation Metrics In order to evaluate the performance of a statistical learning method on a given data set, it is necessary to measure how well its predictions actually match the observed data. This section presents some metrics used in this project. 20 2.3.1 Mean Squared Error (MSE) The Mean Squared Error (MSE) its one of the most common measures when fitting a regression model. It assesses the average squared difference between the observed and predicted values. When a model has no error, the MSE equals zero. As model error increases, its value increases. The MSE value is calculated by the expression: MSE(y, ˆy) = 1 n n X i=1 (yi−ˆyi)2(2.5) Where: •yirepresents the ith observed response value; •ˆyirepresents the ith response value that is predicted; •n represents the size of the sample. 2.3.2 Coefficient of Determination (R2) The coefficient of determination is used to identify the strength of the model, it captures how well the predictions match the observations, or how much of the variation in the observed data is explained by the predictions. Usually the R²varies between 0 and 1, being expressed as a percentage (the closer to 1, the better). According to James et al. (2013) [6], the R2can be calculated by the following equation: R2(y, ˆy)=1−RSS TSS (2.6) With TSS = n X i=1 (yi−¯yi)2(2.7) And ¯yi=1 n n X i=1 yi(2.8) Where: 21 •yirepresents the ith observed response value; •ˆyirepresents the ith response value that is predicted; •RSS represents residual sum of squares defined by Equation 2.3; •TSS represents the total sum of squares; •¯yirepresents the arithmetic mean of all values of the response variable in the sample; •n represents the size of the sample. 2.3.3 Person Correlation Coefficient (r) The Pearson correlation coefficient is the most common way of measuring a linear correlation. It is a number between –1 and 1 that measures the strength and direction of the relationship between two variables. r=Pn i=1(xi−¯x)(yi−¯y) pPn i=1(xi−¯x)2pPn i=1(yi−¯y)2(2.9) Where: •xi, yiare the individual sample points indexed with i; •¯x=1 nPn i=1 xi(sample mean) and analogously for ¯y; •n represents the size of the sample. 2.3.4 Spearman’s Rank Correlation (rank) The Spearman rank correlation coefficient is a nonparametric measure of the monotonicity of the relationship between two datasets. [7] rank =Pn i=1(R(xi)−R(¯x))(R(yi)−R(¯y)) pPn i=1(R(xi)−R(¯x))2pPn i=1(R(yi)−R(¯y))2(2.10) Where: •R(x) and R(y) are the ranks of the x and y variables; •R(x) and R(y) are the mean ranks; •n represents the size of the sample. 22 2.4 Cross-Validation Cross-validation is a technique used to evaluate the capacity of generalization of a prediction model. In this project, two Cross-Validation techniques are adopted: k-Fold Cross-Validation and Leave-one-group-out Cross-validation. 2.4.1 k-fold Cross-Validation This approach involves dividing the set of observations into k groups, or folds, of approximately equal size. The first fold is treated as a validation set, and the method is fit on the remaining k - 1 folds. The mean squared error, MSE, is then computed on the observations in the held-out fold. This procedure is repeated k times; each time, a different group of observations is treated as a validation set. At the end, you can estimate the model error by averaging the errors of all the validations performed. (JAMES et al., 2013 [6]). Figure 3: A schematic display of 5-fold CV. A set of n observations is randomly split into five non-overlapping groups. Each of these fifths acts as a validation set (shown in beige), and the remainder as a training set (shown in blue). The test error is estimated by averaging the five resulting MSE estimates. Source: JAMES et al., 2013 [6]. Note that in some approaches, besides the segregation of the data set into training and validation subsets, it is proposed the separation of a holdout subset (test). In such method, the test subset provides a final estimate of the machine learning model’s performance after it has been trained and validated. 23 Figure 4: Separation of the sample into training, validation and test subsets. Source: HASTIE, TIBSHIRANI e FRIEDMAN (2009) [5]. 2.4.2 Leave-one-group-out Cross-Validation This approach provides train/test indices to split data such that each training set is comprised of all samples except ones belonging to one specific group. Compared to the 5-Fold Cross Validation, it validates the machine learning model more times resulting in more precise metrics. However, it is more computationally expensive and time consuming. 2.5 Literature Review 2.5.1 Yeh, C. et al. Using publicly available satellite imagery and deep learning to understand economic well-being in Africa (2020) In this article, public satellite imagery was used to train a combined deep learning model to estimate socioeconomic conditions over space and time in several regions of sub-Saharan Africa. For this purpose, two ResNet-18 [8] models were trained separately on multispectral daytime imagery from 30m/pixel Landsat and <1 km/pixel nighttime lights imagery. Then, the final layers of the separate models were concatenated in a final fully connected layer and a ridge regression was used to fine-tune the model. To gather the dataset, the authors assembled data on asset wealth using the Demographic and Health Surveys (DHS) [9] conducted between the years 2009 and 2016, as well as an independent smaller set of household level panel data, the Living Standards Measurement Surveys (LSMS) [10]. With this, the authors were able to establish the villages of the African countries as ‘clusters’. Then, Landsat surface reflectance and nightlights images centered on each cluster location were obtained with 30m/pixel resolution, spanning 6.72km on each side. The model was able to explain on average 70% (R²= 0.7) of the variation in ground- 30 4 CASE STUDY: VALE DO RIBEIRA Vale do Ribeira (VR) was chosen as a case study due to its data availability and its ecological and economic importance. It is located in the southern region of Brazil, encompassing the states of Paran´a and S˜ao Paulo. With 28,306 km², it has 30 municipalities and it is the largest continuous area of preserved Atlantic Forest in Brazil, being declared as a Natural Heritage of Humanity by UNESCO (United Nations Educational, Scientific and Cultural Organization) in 1999. Finally, even though it has a high economic development, this area is considered economically poor and has the lowest human development index in the state of S˜ao Paulo [16], being very interesting to study. Figure 5: Vale do Ribeira and its municipalities. Source: Machicao et al. [17] However, due to the data image collection methodology proposed in this project, the 31 number of images obtained considering only this region was very small, as explained in section 6.1.2. Therefore, the other cities of S˜ao Paulo (SP) and Paran´a (PR) were also included in the training of the model, totalizing 1044 municipalities. PART V METHODOLOGY 33 5 OVERALL WORKFLOW The workflow used in this project to replicate the work of Yeh et al. [1] was divided in four steps: 1. Data Acquisition and Aggregation: overlay of a grid over the region of interest, with blocks of approximately 45km²in area. For each block, an average income indicator was calculated as the weighted average of the income of the municipalities, provided by the IBGE census of 2010, and their respective areas contained in the block. Multispectral daytime imagery (MS) and nighttime lights imagery (NL) were obtained through the Google Earth Engine API. 2. Feature extraction: pretrained ResNet-18 networks models were modified to adapt multi-band satellite imagery and used to extract the feature vectors of the images. The loss function used was mean squared error (MSE). The training followed a 5-fold cross-validation. 3. Feature Concatenation and Ridge Regression: concatenation of the feature vectors and refinement of the model using a ridge regression. The loss function used was MSE. The training followed a leave-one-group-out cross-validation. 4. Results analysis, in which performance metrics such as the coefficient of determination (R2) were calculated and graphs were drawn to better analyze the results. 34 Figure 6: Workflow diagram. Source: Author’s Compilation. 35 6 DATA ACQUISITION AND AGGREGATION 6.1 Data Acquisition As mentioned in section 2.5.4, in the preliminary study done in the paper SOUSA, I et al. (2022) [18], the data collection methodology consisted of obtaining the images centered on the census sectors of the municipalities of Vale do Ribeira. This approach resulted in a low performance of R2= 0.289, attributed to the considerably small size of the dataset and the possibility of overlap between the images. 6.1.1 Analysis of the Vale do Ribeira Region This section aims to describe the data analysis done in the Vale do Ribeira region in order to better analyze whether there is superposition between the images. The calculations were done using the geometric information of the territories and its rural/urban classification, provided by the IBGE in its official website. With that, the area of each census sector was calculated and the distribution of the census tract areas was obtained. Figure 7: Boxplot of census areas by urban or rural types. Source: Author’s Compilation. 36 It is possible to visualize that there is a major discrepancy between the size of the urban areas when compared to the size of the rural areas. Due to the proximity of the urban tracts and the area of 45 km²covered by one image in the model’s methodology (images resolution of 30m/pixel and CNN with input 224x224 pixels - as explained in section 6.1.4), the clusters selected in the model will probably be composed of multiple urban tracts. This implies that the same sector can be present in several images, harming the interpretability of the results. This overlap between images can be exemplified in Figure 8 and can also happen with less intensity for rural areas. Figure 8: Cluster boundaries of the municipality of Apiai, SP. Source: Author’s Compilation. 6.1.2 Grid with tiles of 45km² To overcome the problem of overlapping between images, a new methodology for the imagery collection is proposed. The idea is to overlay of a grid over the region of interest, with blocks of approximately 45km²in area. To adapt the socioeconomic indicator of income (see section 6.1.3) to the new methodology, for each block an average income indicator was calculated as the weighted average of the income of the municipalities and their respective areas contained in the block. With this, the total number of images for the Vale do Ribeira region was approximately of 500 images, which is very low. Therefore, to increase the number of the images, the other municipalities of SP and PR were also included in the training of the model, totalizing 9748 images. 37 Note that Albers Equal-Area Conic was chosen as the reference map projection for the grid. The Albers Equal-Area Conic projection parameters follow the cartographic standards defined by IBGE, with SIRGAS-2000, the official planimetric datum of Brazil, as the reference geodetic system. It is recommended for equal-area mapping, being best suited for land masses extending in an east-to-west orientation at mid-latitudes [19]. In this way, the blocks were automatically generated with the same corresponding sizes, always defined by Albers Equal-Area Conic plane coordinates. 6.1.3 Socioeconomic Indicators The socioeconomic proxies utilized in the model were obtained using data gathered from IBGE from the 2010 census, the latest census available. The IBGE and Atlas Brazil [20] publish a municipal Human Development Index (HDI) every decade based on the census data and other surveys. The HDI is a statistical index composed of income, longevity, and education indicators. It was developed and compiled by the United Nations to measure the various levels of social and economic development of countries. In this work, we focus only on the income indicator of this index at the municipal level, which can be found here. Note that since the number of images gathered only depends of the spacial extend of the area of interest, there is no longer the need of developing the model at the census sectors level. Instead, the municipality level is used, eliminating the necessity to adapt the indicators to a different granularity than the one provided by IBGE. 6.1.4 Satellite Imagery The satellite imagery chosen was multispectral daytime imagery from 30m/pixel Landsat and <1 km/pixel nighttime lights imagery. The motivation for choosing nighttime lights imagery was that earlier studies demonstrated that they can serve as an indicator of economic activities at night [21] and, when this indicator is compared across regions and over time, it can be used to measure economic performance. Regarding the choice of daytime imagery, earlier works demonstrated that high-resolution (<1 m/pixel) imagery from private-sector providers can be used to measure spatial variation in local economic outcomes in several developing and middle-income countries. Nevertheless, high-resolution imagery remains very expensive nowadays, so the idea is to focus on more coarser public imagery that can be found for free. 38 The methodology adopted by Yeh et al. (2020) [1] focuses on obtaining the Landsat surface reflectance and nightlights images centered on each cluster location. The same is done in this project, in which the clusters are the cells of the grid that was overlayed in the regions of SP and PR. Figure 9: Grid zoomed to the Vale do Ribeira region. Source: Author’s Compilation. Following the methodology of Yeh et al. (2020) [1], it was obtained a 3-year median composite for the Landsat surface reflectance. This composite was created by taking the median of each cloud free pixel available during the period of 3 years, which in this case is 2009 to 2011, since the Brazilian census survey was done in 2010. As the author’s describe in their original paper, the motivation for using three-year composites was twofold. First, multi-year median compositing has seen success in similar applications as a method to gather clear satellite imagery. Second, the outcome we are trying to predict tends to evolve slowly over time, and we similarly wanted the inputs to not be distorted by seasonal or short-run variation. The Landsat surface reflectance imagery was captured by the Landsat 5 and Landsat 7 satellites, with seven bands which we refer to as the multispectral (MS) bands: RED, GREEN, BLUE, NIR (Near Infrared), SWIR1 (Shortwave Infrared 1), SWIR2 (Shortwave Infrared 2), and TEMP1 (Thermal), all with a spatial resolution of 30 m/pixel. Those are all the surface reflectance bands available for both the satellites in question. See Table 1 39 for a description of the bands mentioned. Name Units Wavelength Description B1 0.45-0.52 µm Band 1 (blue) surface reflectance B2 0.52-0.60 µm Band 2 (green) surface reflectance B3 0.63-0.69 µm Band 3 (red) surface reflectance B4 0.77-0.90 µm Band 4 (near infrared) surface reflectance B5 1.55-1.75 µm Band 5 (shortwave infrared 1) surface reflectance B6 Kelvin 10.40-12.50 µm Band 6 surface temperature. B7 2.08-2.35 µm Band 7 (shortwave infrared 2) surface reflectance Table 1: Description of Landsat 5 and 7 surface reflectance bands. Note that in the paper of Yeh et al., Landsat 8 was also used since images taken from these satellites have been available since 2013 and in the original study the surveys were from several years (2009 to 2017). For comparability, it was also created 3-year median composites for the nightlights imagery (NL). Only DMSP imagery was used, since the first year of availability of the VIIRS dataset at GEE is 2014. The main difference here is that the authors used DMSP and VIIRS because no single satellite captured nightlights for all that period. Note that the nightlight imagery is much more coarse than the daytime imagery we acquired, with a resolution of 927.67 meters. The images are resized using nearest-neighbor resampling to cover the same spatial area as the Landsat images, this is done by default by GEE during reprojection. Finally, both MS and NL images were processed in and exported from GEE in 255 × 255 tiles, with a 30m/pixel resolution, then at the training stage center-cropped to 224 × 224 (the input size of the convolutional neural network architecture), spanning 6.72 km on each side (30 m Landsat pixel size ×224px). Note that the reason why 255 x 255 tiles are used is to have the flexibility of using ‘random crops’ as a form of data augmentation. 6.2 Data Aggregation The files were exported as TFRecords, a simple format for storing a sequence of binary records, including all the bands (MS + NL) but as well other relevant data such as the values of the socioeconomic data labels. Inspecting those records, I was able to obtain 46 error between fourteen values with order of magnitude ranging from 10−2to 109. Figure 14: Example of leave-one-group-out cross-validation for the concatenated model MS+NL. In this case, lowest MSE = 0.001 and best alpha = 16 Author’s Compilation. 47 8 RESULTS To analyze the results, the principal metric chosen was the coefficient of determination (R²). As mentioned in the section 2.4, this coefficient is used to identify the strength of the model, it captures how well the predictions match the observations, or how much of the variation in the observed data is explained by the predictions. Usually the R²varies between 0 and 1, being expressed as a percentage (the closer to 1, the better). This metric is very used to measure the performance of a regression method due to its facility of interpretation, being usually more informative than error metrics that have arbitrary ranges, such as MSE. To evaluate the performance in the same manner as Yeh et al. (2020) [1], the squared Pearson correlation coefficient (r²), the Spearman’s rank correlation coefficient (rank) and the MSE were also used. The r²is powerful to find patterns and relationships in data and the closer to 1, the higher the correlation between the variables. While the squared Pearson’s correlation (r²) assesses linear relationships, Spearman’s rank correlation assesses monotonic relationships (whether linear or not). Furthermore, heatmaps with the values of the real and estimated data were made to provide a visual aid of the results. 8.1 Model The coefficient of determination between the observed and predicted values was R²= 0.4016. Compared to the article of reference of Yeh et al. (2020) [1], where the authors obtained the best case as R2= 0.70, the result is not ideal. Nevertheless, from the plot produced it is possible to visualize that the real and predicted values are proportional. Furthermore, the value obtained is higher than the one found by Tri˜nanes et al. (2020) [14] (R²= 0.35) and than the values found by Tsuru et al. (2021) [23] for the GDP per capita of the states of Alagoas, Paraiba, Rio Grande do Norte and Sergipe. 48 In addition, compared to the preliminary results obtained in my previous study [18], the performance improved significantly and, with the new data collection methodology, the problem of overlapping between the images has been solved. Therefore, it can be concluded that the algorithm is very promising in the task of predicting socioeconomic indicators from satellite images for the municipalities of SP and PR. Figure 15: Regression plot for the combined MS+NL model. The dotted line corresponds to the line of best fit. Source: Author’s Compilation. For comparison, I also ran the ridge regression for the separated MS and NL models, and redid the whole training procedure considering just the RGB bands. See the metrics results in the Figure below. Figure 16: Final metrics results. Source: Author’s Compilation. 49 From this, it is possible to draw some conclusions: •The combined model MS + NL outperformed the other models in all metrics, proving itself as the best model for this task; •The CNN trained on the RGB bands performed poorly compared to the combined MS + NL model, which is expected and confirms the relevance of using multispectral daytime imagery and nightlight imagery when predicting socioeconomic data. Furthermore, for complementing the analysis, the correlation matrix for the results between the models were drawn. Figure 17: Correlation results for the income predictions of all models a) comparison of r²metric and b) comparison of rank metric Source: Author’s compilation. This allows us to make the following conclusions: •There is a high correlation between the MS and the RGB models, but the MS outperformed the RGB, which again is expected since the MS model includes more bands than the RGB one; •Even though the NL and RGB models have a similiar perfomance, the correlation between then is low (r2= 0.29), suggesting that they are good in predicting different features. 50 8.2 Results after calculating the income for each municipality In this section the regression results are analyzed considering the actual values as the income indexes provided by IBGE and the predicted ones as the indexes calculated using weighted average and the values obtained by the model. 8.2.1 S˜ao Paulo and Paran´a 8.2.1.1 Performance Figure 18: Regression plot after calculating the income for each municipality. The dotted line corresponds to the line of best fit. Source: Author’s Compilation. The coefficient of determination for the municipalities of SP and PR is R2= 0.394. This result is very similar to the one obtained for the model without any recalculations, which leads to the conclusion the methodology adopted for the calculation of the grid’s cells average income is well-suited for this project. 51 8.2.1.2 Visual Analysis Figure 19: Heatmaps for a) Real HDI Income indicator and b) Predicted HDI Income indicator. Source: Author’s compilation. To complement this visual analysis, the difference between the real income score y and the predicted ˆy value, defined as d=∥y−ˆy∥, was also calculated and plotted in the heatmap. From this figure it is possible to visualize that the majority of the municipalities yield a small difference between the real and estimated values, for instance 87,55% of them obtained d equal or lowest to 0.05. Figure 20: Differences between the real and predicted labels for the municipalities of SP and PR. Source: Author’s Compilation. 52 8.2.2 Vale do Ribeira 8.2.2.1 Performance Figure 21: Regression plot for the municipalities of Vale do Ribeira. The dotted line corresponds to the line of best fit. Source: Author’s Compilation. The coefficient of determination for the municipalities of Vale do Ribeira is R2= 0.267. This result is significantly smaller than the one obtained considering all the municipalities. This is expected and can be justified due to the small number of samples: only 30 municipalities are part of the Vale do Ribeira region. 53 8.2.2.2 Visual Analysis Figure 22: Heatmaps for a) Real HDI Income indicator for VR and b) Predicted HDI Income indicator for VR. Source: Author’s compilation. To complement this visual analysis, the difference between the real income score y and the predicted ˆy value, defined as d=∥y−ˆy∥, was also calculated and plotted in the heatmap. From this figure it is possible to visualize that the majority of the municipalities yield a small difference between the real and estimated values, for instance 73,33% of them obtained d equal or lowest to 0.05. Figure 23: Differences between the real and predicted labels for the municipalities from VR. Source: Author’s Compilation. PART VI CONCLUSION 55 9 CONCLUSION In this report, public satellite imagery was used to train a combined CNN model to estimate socioeconomic data over space in the states of S˜ao Paulo and Paran´a, with the goal of particularly analyzing the Vale do Ribeira region due to its economic and environmental importance. With an R2of 0.4016, the model yields a low performance when compared to the one of Yeh et al. (2020) [1], which attempts to replicate (R2= 0.70) . Nevertheless, it is still promising since the actual and predicted values are clearly proportional. Moreover, when compared to the best results of other studies that have analyzed socioeconomic data using machine learning and satellite imagery in Brazil, it presented a higher R2value by approximately 0.05. Finally, it successfully satisfies the goal of improving the performance of R2= 0.289 obtained in my preliminary study (SOUSA, I et al., 2022) [18]. In conclusion, although not accurate enough to replace field survey data, the model performed well and can be used to help analyze the socioeconomic situation of the Brazilian territory in years when official government survey data are not published, supporting the planning and evaluation of public policies. Furthermore, it is a step towards understanding how convolutional neural networks, daytime and nighttime multispectral imagery can be used in the task of forecasting socioeconomic data, and thus a stimulus for the development of future work in the field of predicting data through satellite imagery.