Estimating the size of urban populations using Landsat images: a case study of Bo, Sierra Leone, West Africa

Background This is the third paper in a 3-paper series evaluating alternative models for rapidly estimating neighborhood populations using limited survey data, augmented with aerial imagery. Methods Bayesian methods were used to sample the large solution space of candidate regression models for estimating population density. Results We accurately estimated the population densities and counts of 20 neighborhoods in the city of Bo, Sierra Leone, using statistical measures derived from Landsat multi-band satellite imagery. The best regression model proposed estimated the latter with an absolute median proportional error of 8.0%, while the total population of the 20 neighborhoods was estimated with an error of less than 1.0%. We also compare our results with those obtained using an empirical Bayes approach. Conclusions Our approach provides a rapid and effective method for constructing predictive models for population densities and counts utilizing remote sensing imagery. Our results, including cross-validation analysis, suggest that masking non-urban areas in the Landsat section images prior to computing the candidate covariate regressors should further improve model generality.


Introduction
In resource-limited environments, it is desirable to be able to rapidly estimate the density of local populations. The ability to estimate population sizes is important in places where population growth is relatively high and census data are relatively old. Many of these locations are in urbanizing areas of low-and lower-middle-income countries.
Such estimates are invaluable for health planning, refugee support [1], epidemiological modeling [2], and for state and municipality-sponsored allocation of public resources and services. Most commonly, such estimates are made using some combination of aerial imagery and local survey data. In two recent papers, we used groundtruth survey data from Bo, Sierra Leone, to model several different approaches for estimating section (neighborhood) population. As a function of sample size, comparisons were made between the uncertainty of the estimated population based on the mean occupancy of residential structures and the mean number of individuals per square meter of rooftop area [3,4].
Both studies required only a limited amount of survey data, in addition to estimates of the total number of residential structures in a region of interest. Methods that utilize rooftop area additionally require estimates of individual and total rooftop areas in regions

Analytical approach Estimating population densities
The use of satellite imagery for appraising land utilization, including population density estimation, is not novel. For a brief overview, see "Appendix 1". Our analysis uses selected TM measures of mean spectral reflectances (intensities), pixel-level spectral transforms, and diverse measures of spatial variability [that is, measures of texture] identified by Harvey [5] as candidate covariates. Because our population data are restricted to the measured populations of 20 sections in Bo, we test several different protocols for building and testing the regression models when sample sizes are small.

Estimating section populations
Given the estimated population densities for each section, the total population of the surveyed areas may be roughly estimated as the scalar product d, Area of the estimated population densities d and the measured section areas. The implicit assumption is that the population density is relatively homogeneous within each section. This assumption is not satisfied for some of the sections surveyed, although the regression models developed are still surprisingly accurate.

Table 1 Bo municipal survey data
Residential and household survey data for 20 municipal sections of Bo, ordered by population, showing the persons per municipal section, section area, and the total number of residential and non-residential structures [2,3] The min, mid and max values (1) Section (

Three questions to be addressed
After a preliminary discussion of materials and methods, we develop a regression model for estimating the population densities of the 20 sections. In particular, we address the following three questions: 1 Using the Landsat TM data to define a candidate set of independent variables, can we build one or more regression models for accurately estimating the measured population densities of the selected Bo City sections? The raw TM data consist of mean band-specific pixel-level intensity measurements for each section. 2 Can we then estimate the entire population of the sections in the dataset, given the estimated population densities of the individual sections, and the measured section areas?
3 Applying the k − 1 cross-validation method (also referred to as "Leave one out cross-validation, " or LOOCV), how effectively do these regression models generalize to estimating the population density of a section deliberately omitted from the LOOCV training set?

Survey methodology
The survey methodology is summarized in [3]. The data collection protocols for human subjects were approved by three independent Institutional Review Boards: Njala University, George Mason University, and the U.S. Naval Research Laboratory. Household data were collected from one adult representative of each participating household after obtaining written informed consent  18:16 from that individual. Most residential structures were home to multiple households. To be defined as a resident of a household, a child or adult had to use the structure as sleeping quarters most nights. Family members who usually worked in other locations or attended boarding schools were not considered to be residents. The total population of each section was calculated by adding up the total number of residents in each residential structure. The data for the 20 surveyed sections listed in Table 1 have already been published in open-access literature.

The Landsat thematic mapper (TM)
Landsat 5 was an Earth-observing satellite launched on March 1, 1984, into a near polar orbit at an altitude of 705 km, for collecting imagery of the Earth's surface. It was decommissioned in January 2013. Landsat 5 instrumentation included a Thematic Mapper (TM) with an optical-mechanical "whisk broom" (along-track) scanner [6,7]. The scanner's mirror system bi-directionally swept the TM's detectors along a line transverse to the northsouth path of flight. The archived Landsat 5 TM scenes have an area of 170 km north-south by 183 km east-west (i.e. 106 mi by 114 mi). [8].
All data used in this article were derived from the scene LT52010542011001MPS01 [9] with the indicated path (201), row (54), date and year (2011/1/1). Publication of this imagery is in full compliance with guidelines [10,11] authorizing the use and dissemination of USGS satellite imagery. The year 2011 was selected because the survey data for the population sections were collected in the same year [11]. Although Landsat 7 could have potentially provided more refined data, a failure of the TM scan line corrector (SLC) corrupted the scenes collected at the required dates (2011) and locations [12].

Correcting for atmospheric effects
The Landsat sensors capture reflected solar energy. The Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [13] is a software system for processing Landsat imagery to calculate the reflectance from the earth's surface. A LEDAPS-processed dataset is available for the desired imagery [9]. The 3 major steps in LEDAPS processing are: 1 As a function of the band-specific sensor gain and bias, convert the Landsat sensor outputs to sensor spectral radiances, the energy reaching the sensors. 2 As a function of the earth-sun geometry and the mean solar exoatmospheric spectral irradiances, convert the spectral radiances to the Top of the atmosphere (TOA) dimensionless reflectances. The latter is the dimensionless ratio of reflected energy to total energy. 3 Estimate the reflected energy measured at the earth's surface, rather than at the top of the atmosphere, by removing the interference imposed by the atmosphere itself on both the incoming and reflected solar radiation. This step requires correcting for wavelength-specific atmospheric scattering as well as masking and correcting for distortions imposed by cloud cover, shadows, and reflections from water.

TM data visualization
The TM data are multispectral, and each scene was captured in 7 different bands. Table 2 shows the bandwidth, resolution, and nominal utility for each of the 6 Landsat TM bands [14,15] used in this study. The data from the different bands are usually combined to create complex images that enhance specific features of the target region. By mapping each band onto the visible colors red (R), green (G), and blue (B), the individual Bo City band images can be combined into different composite images [15]. The mappings are specified by indicating the sequence of bands assigned to the visible composite colors R, G, and B. In the "NIR" (near infrared) (bands 4, 3, and 2) mapping shown in Fig. 2, Band 4 is assigned to composite color R. Because vegetation reflects brightly in the NIR band 4, the vegetation surrounding Bo City appears to be bright red.

Pixel-level section representations
Six of the seven Landsat 5 TM bands were utilized. Band 6 in the TM sensor is emittance (temperature), and not normally used in combination with reflectance data; omitting Band 6, pixel-level matrix representations of the surface reflectance from each section can be made for each band using the LEDAPS corrected data.
For example, the pixel magnitudes measured in Band 3 are shown in Tables 3 and 4 for Moibawo Farm and New London. The min-max normalization algorithm [16] was applied to rescale the sensor data between 0.0 and 1.0; in the two tables, these normalized values are multiplied by 100.0 to facilitate readability. Comparable visualizations could be made for each of the other bands. The area of New London is approximately 0.60 km 2 , and Moibawo Farm is 0.50 km 2 . The mean, standard deviation, and variance of the min-max pixel distributions defined the normalized variables nb [mean value of normalized LEDAPS-corrected pixel magnitudes], nbs [standard deviation], and nbv [variance] for these two sections in Band 3.
The resolution of the pixels for the 6 selected bands, including Band 3, is 30 m. There are 670 non-zero pixels in the New London section, and 559 pixels in the Moibawo Farm section. The areas estimated from these pixel distributions are consistent with the areas estimated from the shape files (i.e. map boundaries). Let NP i designate the number of pixels for each distribution, and PA the pixel area, which is always 900 m 2 . The ith section Area i is then: The mean value of nb i , the normalized pixels for Band i, is: The variances and standard deviations for these distributions can be similarly derived.

Estimation methods
The premise of this paper is that low-dimensional subsets of variables derived from Landsat data can be used to construct accurate regression equations for estimating the population densities of the 20 surveyed sections. In this section, we will describe the datasets, methods, and metrics that were used. Figure 1  (1)

The TM covariate dataset
In his study [5], Harvey proposed a large set of candidate Landsat TM covariates for estimating population densities in Australian census districts. He reduced this preliminary set of variables to a low-order set of covariates through a complex sequence of model testing. We used Harvey's full set of proposed candidate variables for our regression analysis. An obvious objection is that Harvey's models were tailored to estimate population densities in the urban northern latitudes of Australia, whereas we were fitting our models to estimate population densities in a region where much of the population resides in informal settlements. However, we utilized the full instantiated set of candidate variables, with some exceptions to be noted, as input to our model selection algorithms. At no point did we use the reduced sets of candidate variables or the specific regression models that were trained and tested against Australian census data. The regression covariates selected during model construction therefore reflected the unique attributes of our Bo City dataset.
Our methodology also used improved methods. Rather than implementing the TOA and atmospheric corrections manually, as Harvey was required to do given the technical restraints at the time, we were able to use the LEDAPS-processed imagery provided by USGS. We also used Bayesian MCMC (Markov chain Monte Carlo) sampling to find the variables for our regression models, rather than step-wise regression, although the latter remains a viable approach.

TM variable definitions
Multiple candidate variables were calculated for each of the 20 Bo City sections. To simplify the notation, the index term for the section [i.e. a number between 1 and 20] has been omitted, as there are no variables that are functions of more than one section. See Table 5.
Let p denote the number of pixels sampled in a given section and b i n denote the value of the Landsat thematic mapper (TM) sensor measurement of the nth pixel in band i. For each pixel, measurements were made in bands 1,2,3,4,5 and 7; i is restricted to these values. Additional candidate covariates were then derived from the LEDAPS-corrected pixel-level intensity measurements. Table 5 summarizes the 3 datasets used in subsequent analysis: (1) non-spectral transforms, (2) spectral transforms, and (3) the total combined dataset. There are 379 total variables, with a subset of 304 spectral transforms and 75 non-spectral transforms. The definitions and equations for all variables in Table 5 are given in "Appendix 2". The initial set of 379 candidate covariates was substantially reduced prior to initiating the regression analysis per se, using methods described below.

The TM data array
The 20 measured observations of persons per section, in combination with the measured section areas, yield the dependent variables d i = Persons i Area i=1,...20 . Our model estimates d i as a function of the Landsat TM measurements. The Landsat Thematic Mapper (TM) measurements and derived variables can be arranged in an array with 20 rows and 379 columns. Each row denotes a Bo City section, and each column corresponds to one of the 379 variables derived from the Landsat TM data. This array is shown schematically in Table 6. Two columns of demographic variables (section name and d = population density ) precede the 379 columns of TM data.

Software development
The regression simulations and auxiliary plotting functions were written in the programming language R by the first author. Support functions from multiple R libraries were used, particularly [17]. The second author developed additional R code for processing the Landsat imagery, and produced the 20 by 379 matrix of Landsat TM derived products.

Regression methods
We will now summarize the major steps: 1. Data reduction. We began with a data array containing 379 candidate regression covariates. This was reduced to an array of 159 covariates prior to conducting the regression analysis. First, the subset of 304 spectral transforms alone was found to yield a good solution.  [17][18][19]. A brief summary of the methods used is provided in Appendix 3. The best single equation found for estimating √ d during the stochastic sampling was transformed to a conventional linear multiple regression equation.  [20] into the original parameter space as d i . The goodness-of-fit of the regression equation for estimating d could then be evaluated. The population of each section was also estimated. 5. Cross-validation. "Leave-out one cross-validation" (LOOCV) [21] was used to quantify how well the regression equation generalizes to estimating observations that were not included in the training set.  Non-spectral ratio of the difference-to-the-sum of the

Data reduction
The original Landsat data array has 379 candidate regression covariates. Reducing the size of this dataset should increase the effectiveness of the MCMC sampling algorithm by reducing the size of the regression model search space. PCA (Principal Components Analysis) is often used to reduce a large dataset prior to subsequent analysis, but PCA transforms the original variable set by mapping combinations of variables onto a new coordinate system. We wanted to identify the individual Landsat variables which were most critical for estimating the population density, so PCA was not an appropriate method. Two preliminary steps were used to reduce the dataset prior to MCMC sampling. First, by trial-and-error we found that all of the covariates selected were from the subset of Landsat variables defined for spectral (i.e. inter-pixel) transforms (Table 5). Using only the spectral transform subset of variables reduced the size of the data array from 379 candidate covariates to 304 candidate covariates. Second, we removed a member of each pair of "identical" covariates whose Pearson correlation was 0.99 or greater [22]. The set of 304 covariates was reduced to a set of 159 covariates without any degradation on the quality of the regression models. See Table 7.     Table 8 gives the parameters for the best regression model found for estimating √ d using the sampling protocol summarized in Appendix 3. Given the low values of the VIF, there is no significant multicollinearity between the selected variables (col 7). The fit of the model is excellent: R 2 = .9951 and R 2 . adjusted = 0.9928 , on 6 and 13 degrees of freedom. See Table 9. The regression was run on the transformed population density variable (i.e. on the square root of the population density). The square root transform generated a more linear relationship between the Landsat sensor readings and the dependent variable of section population than a log transform or no transform, which contributes to the high values of R 2 for the transformed variable. See Fig. 3 for a comparison of regression plots made using log and square root transforms and no transform at all.

Regression analysis
One indication that a good solution has been found in the sample space is that the MCMC sampler frequencies and the analytical posterior marginal likelihoods both converged. For 10 7 iterations, the correlations were almost perfect (0.9657) between the empirical and analytical distributions. Figure 4a shows the back-transformed estimates of the populations densities d i , plotted as a function of the measured population densities. The regression equation in Table 8 [20]. Panel (B) shows the estimate of the population obtained by multiplying the back-transformed estimate of d by the measured section areas:

Relative proportional error RE
Harvey [5] recommends the Relative or Proportional Error as a measure of fit, rather than R 2 , and we will provide these values for the regression error.The Relative Error, which will be abbreviated here as the RE), is defined as: This measure is the same for both the population and the population density. It can be calculated for the estimated transform of the population density √ d and  Table 9 lists the REs by section, as well as d i and p i . The absolute value of the RE is shown in Fig 5. The fit is very good, and the median absolute RE is 8.0%. The REs for all sections is less than 20%, with the exception of Moibawo Farm, where the population density is underestimated by 24%.

LOOCV cross-validation
In our current study, the number of aggregated population observations is 20. This is insufficient to divide the observations into training groups (sets) and test groups (sets), as is normally done for cross-validation. As an alternative, we used k-1 cross-validation, where k = 20 . Assume that a regression equation with x independent variables has been found for estimating d i for all k observations, where each observation is the measured population density d i . Referring to Table 8, x = 6 and n = 20 . There are k LOOCV (Leave-Out One Cross-Validation) tests that can be constructed and executed. An obvious drawback is that there is only a single observation available for estimation on each trial. In each of the n LOOCV trials, a single observation d j was omitted from the dataset. Using the same x independent variables, a new regression model was fitted to the remaining n − 1 population density observations d i . The reduced model was then used to estimate the single omitted population density d j . This process was repeated for all n trials. A different regression equation was parameterized for each of the n trials, but the same set of x independent variables was always used. Figure 5 shows the absolute value of the percentage relative error RE for each section. Table 10 shows the details of the calculation, as well as the RE for the transformed population density √ d . In Fig. 5, the bar charts show both the relative error (RE) for the estimation of the population density by section and the absolute values of the RE for the cross-validation tests.
Although the median absolute value of RE for the backtransformed estimate is only 11.14%, the model failed to generalize (i.e. cross-validate) well in at least 3 cases. d for Moibawo Farm was underestimated by almost 50%, New York was underestimated by over 26%, and Roma was overestimated by about 41%. It is difficult to discern a simple pattern in the outliers. Moibawo Farm, like Reservation, has large open non-residential areas. But if this caused the underestimation in population density, the estimate for Reservation should have been similarly affected.

Discussion
The model used in our research was specific to the 20 sections that we studied. The cross-validation study demonstrates that the six covariates in the regression model could be used to construct 19 separate regression equations for estimating the population density d of an omitted section, although there were several outliers noted. The model has not yet been tested in other urban areas with different patterns of residential structures, building materials, roads, or other characteristics, and it is likely that adaptation to the model and variables would be required.
Because the MCMC sampling of the solution space is stochastic and incomplete, the regression model Table 10 This table summarizes the

results from the k − 1 'LOOCV' cross-validation analysis
For each section, col (1) is the section population, Col (2) is the area in km 2 , and column (3) is the population density d = persons/area . Col (4) is the transform √ d , col (5) is the estimate √ d , and (6) is the back-transformed estimate of d . Cols (7) and (8) are the error (Er) and RE) of the estimated variable √ d ; and cols (9) and (10)  summarized in Table 8 is not unique, although the "top model" solution was very effective for predicting d. A fixed random number seed was used in the simulations to enable the replication of results between simulations. Given different initial random number seeds, or alternative numbers of sampler iterations, alternate solutions could have been found. All six of the selected regression variables are measures of covariate spatial variation (variance, coefficient of variation, and standard deviation), as can be seen in Table 8. These measures denote spatial variations in brightness between relatively large 30 m pixels. A typical Bo residential structure is smaller than a single 30 m pixel, and these measures of spatial variation cannot capture fine-scale modulations in reflectance within individual rooftops. The TM resolution is also insufficient for the application of feature extraction algorithms for explicit capture of rooftops or other structural boundaries [23,24].

Statistical significance of individual regressors
The stochastic nature of the simulation does not, however, diminish the significance of the variables selected with respect to their relative importance in the sample space as good candidate predictive variables (i.e. regression covariates) for estimating d. Four of the PIP (posterior inclusion probability)values are were close to 1.0. It is highly likely they would be included in any of the 1000 best-fitting models that were retained by the sampler, as well as in the "top" model. (The number of top models tracked by MCMC sampler is user-selectable.) Four covariates (nb7v, r_sp37 , nb1v, and ch245c) out of the six in the regression equation have PIPs close to 1.0. The high PIP values indicated that all four variables were included in almost every one of the 1000 best-fitting models tracked by the Bayesian MCMC sampler, which implies that the selection of these four variables was robust. The PIPs of the remaining two covariates were 0.47 and 0.54; each was retained in about half of the 1000 best regression models. The PIP is also proportional to Schwarz's Bayesian information criterion (BIC) [25,26].
Another advantage of our approach is that each of the six regression covariates was calculated directly from Landsat imagery, rather than as a transform of multiple Landsat variables. In data reduction methods such as PCA (Principle Components Analysis), the significance of the individual Landsat variables may be obscured by the complex mapping of the individual variables into the transform space.

Interpreting spectral signatures
The variables and combinations of variables that were selected for the regression model are consistent with our understanding of the natural world. Within this scene, one can see that the unpopulated areas are heavily vegetated whereas the populated areas surveyed are a combination of tarpaulin and zinc/aluminum roofs, paved and unpaved driving/walking surfaces, as well as bare earth and vegetation between structures. The interpretation of why specific combinations of variables were selected is somewhat conjectural.
With the exception of r_sp37 , all of the covariates are measures of spatial variation ("texture"), rather than measures of brightness. For the band 7 covariate nb7v, a high variance is negatively associated with d; this band can aid in the differentiation between soil types and minerals, and is also sensitive to water content. ch245c is the coefficient of variation (CV) for a cylindrical transform of bands 2, 4, and 5; this tri-band mapping onto a single value constitutes a form of data compression. All 3 bands reflect vegetation brightly, but it is the CV that appears to be positively associated with the population density.
A characteristic of regional statistics, like the ones we used, is that each region has different fractional amounts of the previously stated ground cover materials. Manmade materials often reflect more in the infrared portion of the spectra (e.g. NIR, SWIR1, and SWIR2) as compared to vegetation, and vegetation absorbs more light in the visible portion of the spectra (e.g. blue, green, red) as compared to soil and man-made materials. Armed with this knowledge, we can infer that the multiple variables used in the regression analysis are differentiating the natural, vegetated areas from the built up regions to deduce population density in the region.
The inclusion of the blue band is present in three of the variables: nb1v, r_sp15s , and r_sp14c in Table 8. This seems noteworthy, given the interaction between blue light and Rayleigh scattering as well as Mie scattering. Particulates of various sizes in the atmosphere can either selectively scatter shorter wavelengths (e.g. blue and violet via Rayleigh scattering) or scatter light over a broader wavelength range (e.g. Mie scattering). As part of our future research, we would like to examine how blue light is scattered as a result of particulates in the atmosphere over urban areas as compared to that of densely forested areas, and to see if this is a critical factor for interpreting spectral signatures.

Correcting for non-homogeneous population density
An implicit assumption of this approach is that the population density is relatively homogeneous within a section. This assumption can be problematic in at least 3 ways: 1 If an area (section) is primarily wild vegetation or barren soil, it violates the assumption that the population density is relatively uniform within an area. If so, the spectral statistics for a section may primarily be a function of an "empty" region on the ground, rather than being representative of an area populated (although perhaps sparsely) with built structures and associated property. The Bo City section Reservation provides an extreme example of both issues. This section is essentially a large swamp, with a small number of buildings at the perimeter [3] that were originally constructed for government use. 2 A predominance of non-residential buildings within an area may confound residential and non-residential regions. Other than collecting additional survey data, or utilizing local knowledge to annotate the section imagery, there is no obvious way to differentiate between residential and non-residential structures in the Landsat imagery. 3 Variation in rooftop materials can confound the sensor interpretation within a given area. However, in the 20 sections surveyed, we would not expect great variation in sensor readings attributable to differences in roofing materials. Of the 1165 residential structures surveyed in the 20 sections, 1156 had zinc roofs, 8 had tarpaulins, and one was "other. " For this reason, it is unlikely that we confounded residential rooftops with bare earth or cultivated land.

Land-use/land-cover models
(LU/LC) modeling offers another approach to resolving the consequences of non-homogeneous land use. A LU/ LC model would differentiate between categories of land use in different areas of a section, distinguish between residential and non-residential structures, and allow for differences in rooftop construction. Wilson and his collaborators have developed LU/LC models for Bo City that dramatically illustrate the changes in LU/LC as consequence of forced migration during civil conflict between 1998 and 2002. For example, in [27] see Figures 7 and 8, and the accompanying tables. Although elegant, this approach would require a level of ground-truth data, data fusion, and model development that cannot be achieved within the scope of our approach.

Spatial autocorrelation and image resolution
Spatial autocorrelation methods [28] are not new, but the power of these statistical methods has been enhanced by the advent of high-speed computers, the availability of large GIS datasets [29], and the development of custom software packages that facilitate the work of the analyst [30]. The interactions between spatial entities are usually modeled as a function of adjacency (i.e. contiguities between polygonal representations) and/or distance. The links denoting distance can also be weighted. Both global (e.g. Moran's I) and local (e.g. LISA, Local Indicators of Spatial Association; Geary's c) measures of spatial autocorrelation have been developed [28,30]. There is a significant interaction between spatial autocorrelation patterns and map resolution [31]. As a concise example, Spiker and Warner [32] derived autocorrelation measures for a satellite image of Morgantown, WV, at three different pixel sizes: 0.7m, 15m, and 60m. The local value of Moran's I is sensitive to buildings and other features of the urban infrastructure at high resolution, while at 60m resolution, geographical features (the river primarily, and secondarily land usage with respect to urban vs rural) dominate. The local values of Geary's c support a similar trend.
Since the resolution of the Landsat sensor data is 30 m, we cannot readily analyze the the accuracy of our population estimation methods as a function of image resolution. We also cannot construct and evaluate complete contingency or distance maps for spatial autocorrelation analysis, because our survey data is limited to 20 of 68 sections of Bo City. Given the findings discussed above, it would useful to repeat our analysis using sensor data at different levels of resolution, using more complete survey data. For example, the interaction between spatial autocorrelation patterns for housing (i.e. structure) density, the ground-truth population density, and the estimated population density could all be examined.

Future research Simulated subsampling
One approach to studying the relationship between resolution, spatial autocorrelation, and model accuracy would be to simulate resampling of the surveyed population using a fixed grid size, perhaps with grid squares as small as 500m 2 . The grid size must still sufficiently large to ensure that the population within each grid square is too large to be mapped onto specific dwellings that are within the square. Population maps at diverse resolutions could then be constructed by combining the populations of 2, 3, or 4 adjacent grid squares into single cells. The smaller the cell, the finer the sample granularity would be.
The Landsat measurements, which are currently averaged over the area of each section, would also have to be recalculated for each of the grid squares for each of the grid resolutions. In the bands used, the Landsat sensor (i.e. pixel) resolution of 30 m would still be significantly smaller than the sizes of the reduced sample grid squares. (A pixel resolution of 30 m is still larger than a typical residential dwelling.) The independent variable would be the number of persons per grid cell, and both global and local measures of spatial autocorrelation could be computed. This approach should disclose regions that are locally clustered and spatially correlated, as a function of grid resolution. The Landsat sensor values would also have to be recomputed, roughly matching the resolution of the resampled grid squares. Given a finer grid resolution, we could determine if the relative error RE for the LOOCV cross-validation decreases. It would also be possible to define training sets and test sets for conventional cross-validation testing.
Even given high-resolution subsampling, it would still not be possible to construct a complete adjacency or distance matrix for the current dataset, because only 20 of 68 sections were surveyed. But within contiguous subregions of Bo City, the following two questions could also be clarified: (1) Do patterns of autocorrelation in the sub-sampled ground truth population data present and/ or vary as a function of resolution? (2) If so, do these patterns modify the estimated population density distributions using the Landsat data?

Masking section imagery
A strategy for improving model generalization would be to partially mask the imagery for each section prior to calculating the values of the covariates. The objective is to correct for the non-homogeneity of the population density within certain sections by masking (i.e. omitting) non-residential sub-areas of a section. This requires omitting pixels corresponding to areas of vegetation. This could done manually as proof of concept. Alternatively, the NDVI (normalized difference vegetation index) could be calculated for each section, and pixels that have relatively high positive values [33] could be omitted from further consideration. (Given rasters for Band 3 and Band 4, the NDVI = (Band 4 − Band 3)/(Band 4 + Band 3) ). A limitation of this approach is that it may not mask nonresidential areas that are either barren or dominated by unhealthy vegetation, but the distribution of included and excluded pixels will also be a function of the exclusion threshold selected. The index values range between − 1.0 and + 1.0. An NDVI value of zero or less means that no vegetation is present, and a maximum value of +1.0 is the strongest possible indicator of healthy vegetation at the pixel location. Here again, the objective is to demonstrate a decrease in the cross-validation error by improving compliance with the model's assumptions.

Alternative approaches to cross-validation
The median absolute value of the relative proportional error RE, as defined in 4 and enumerated in Table 9, is about 8.0%. For example, referring to Table 9 for section Roma, RE = (3818.48-3475.00)/3475.00 × 100% = 9.88%. The median absolute value of the 20 values of RE is 8.85%. Conversely, the sum of the estimates of the section populations in Column 6 is very close to the measured value of the total population. While some sections had a lower-than-observed population and others had a higher-than-observed population, the estimated total sum across all sections (25,856) was very close to measured population size (25,954), an error of less than 1.0%.
The generality of the model was tested using LOOCV (k-1) cross-validation. The results here were less satisfactory than for the population density d estimations.
Although the median absolute relative error was only 11.14%, the RE errors of over 40% for 2 of the 20 sections and over 20% for two additional sections. A limitation of the LOOCV cross-validation paradigm was that only a single observation was available for each trial. Extending the training set would reduce the limitations imposed by the small number of 20 observations available. A larger dataset could be partitioned into multiple training sets and test sets; this would provide a far more robust approach to cross-validation.

Alternative estimators
Finally, there is an additional consideration for which we have conducted a preliminary test. The empirical local Bayes estimator (EBL) can provide a useful and effective benchmark, but it is a controversial one [34]. As Zeugner [17] succinctly states, "It does not constitute a real prior since it involves 'peeking' at the data in order to formulate a prior. " Allowing for these limitations, we developed an EBL model using the data set already described. This was done using the BMS package for R [17], as was the preceding work; the spectral data subset was used, with a reduction in highly-correlated variables executed first.
In this case, a 6 variate regressor equation was found, plus the non-zero intercept. See Tables 11 and 12 for details. In Fig. 6, the EBL bar charts show both the relative errors (RE) for the estimations of the population density and the absolute values of the REs for the crossvalidation tests. A comparison of Figs. 5 and 6 show that the EBL is far more effective than the conventional Bayesian model developed within. Specifically, the RE for the estimated population density is much lower (compare Figs. 5a, 6a). The cross-validation RE (Fig. 6b) is greatest for Moibawo Farm (270 persons/km 2 ) and Reservation (273 persons/km 2 ), the two sections with the lowest population densities (Table 1) and the greatest RE underestimations for cross-validation. About half of the footprint for Reservation is bright green wetlands, and Moibawo Farm is heavily forested. The cross-validation RE for Salina, which has a large industrial area surrounding the main road (the "old railroad line"), is overestimated by almost 25% . An interesting research question is which  Fig. 5a, b, except a Local Empirical Bayesian (EBL) estimator was used model will be better generalize to data sets that were not used to condition either model.

Conclusions
The objective of our study was to demonstrate that it is possible to rapidly develop a predictive model for estimating the population density, and the contingent population count, for local neighborhoods in an urban environment using Landsat data. Although some limitations are imposed by the non-homogeneity of population density in several sections, including Reservation and Moibawo Farm, we have succeeded in this objective. An accurate 6-covariate linear multiple regression model was developed for estimating the population density d. Methodological improvements are also suggested, including NDVI masking of section imagery prior to variable calculation, and higher resolution subsampling of the original survey data. Although our approach will probably not be as accurate as methods using high-resolution satellite imagery, if offers a number of advantages with respect to speed and simplicity for the estimation of local populations: 1 It uses LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) pre-processed Landsat sensor data for deriving variable values. 2 It is not necessary to manually (or automatically) extract residential structure outlines or to define GIS layers or geographical features that correlate with residential areas. 3 Only 30 m LandSat data resolution is required, not high-resolution (<10m) imagery. 4 Each of the six regression covariates selected was derived directly from Landsat sensor imagery, rather than being a composite variable, as in principal components analysis. 5 The posterior inclusion probability (PIP), calculated for each covariate, provides a measure of the variable's information-theoretic significance within the top 1000 candidate regression models. 6 The calculations are also relatively speedy, requiring only a few minutes to run 10 6 Markov chain Monte Carlo (MCMC) iterations and less than 30 min to execute 10 7 iterations. All results discussed in this article are from simulations run with 10 7 iterations, following exploratory simulations with 10 6 iterations.
Potential strategies were discussed that will maintain the above advantages while potentially improving the accuracy and generality of the models.

Appendix 1: an overview of methods for satellite-based estimations of population and urban footprint
The use of satellite imagery for appraising land utilization, including population density estimation, is not novel. A study using urban Landsat data for Sydney, Australia, demonstrated that three distinct classes of surface reflectance data corresponded to vegetation, urban residential areas, and urban non-residential areas [35]. Applying Landsat data, Forster also developed a set of pixel-based regression models for predicting the type of land cover [buildings, concrete, roads, trees, grass, water and soil] in Sydney. Correlations of up to 70% were achieved, and up to 85% for tree and road cover over extended areas [36]. More recently, Forster demonstrated that the coefficient of variation of single-band SPOT (Satellite Pour l'Observation de la Terre) HRV (High Resolution Visual) imagery could be used to estimate housing density and average house size in the area of Sydney, Australia [37]. Harvey developed regression models for estimating the population densities and total population of Australian Census Districts (CDs) from Landsat thematic mapper data [5,38]. The first study area included 225 CDs, and the second 138 CDs, which enabled the partitioning of the data into "training" and "test" datasets. See [39] for a review of more recent "bottom-up" population estimation models using high-resolution satellite imagery.
Satellite imagery has also been used to develop models to estimate the spatial footprint of cities in low-and/or middle-income countries, and their evolution over time.
Wilson and his collaborators have developed a series of temporal models for Sierra Leone demonstrating how Land Use and Land Cover (LU/LC) evolve as a function of civil conflict. The spatial LU/LC models developed for the capital city of Freetown, and the metropolitan cities of Bo and Kenema, encompass three multi-year intervals preceding, during, and after a period of civil strife [27]. In a study of the Indonesian island of Java, Patel et al [40] established that Landsat imagery is valuable for tracking urbanization-that is, the growing footprint of cities, rather than the rise in the number of inhabitants-in growing cities. In this study, the Landsat database was accessed via the Google Earth Engine (GEE). Trianni [41] also used GEE to generate estimates of human population size using Landsat data in Brazil, China, and Indonesia. Son [42] and his collaborators used Landsat data to track the expanding footprint of the capital city Tegucigalpa; while Linard et al. [43] tracked changes in the human population over a ten-year period on the Kenyan coast, integrating data from multiple sources, including Landsat imagery and historical census data.

GUF
Aerial imagery, including hyperspectral datasets, can now be supplemented by high-resolution satellite measurements using SAR (Synthetic Aperture Radar) imaging [44]. The highly automated GUF (Global Urban Footprint) utilizes SAR amplitude data from the TerraSAR-X and TanDEM-X satellites, in combination with texture maps derived from the SAR amplitude data [31]. Strong primary and secondary backscattering from vertical structures result in characteristic bright spots in the texture maps. Conversely, natural phenomena such as tall trees and rock formations can result in false positives, while dense groups of low flat-roof houses may produce false negatives. An unsupervised classification system, based on Support Vector Data Description (SVDD), used these paired distributions to identify built-up regions characterized by vertical buildings, with a spatial resolution of 12 m.

GHSL
The GHSL (Global Human Settlement Layer) [45] was derived by processing archived Landsat data. Data collected by the MS (Multispectral Scanner), TM (thematic mapper), and ETM (Enhanced thematic mapper). (In our current study, only TM data was used.) The resulting sequence of maps shows dramatic changes in patterns of urban growth over a 40-year period.

HRSL
The HRSL (High Resolution Settlement Layer) [29,46] effort is a fourth project for mapping the urban footprint. In this effort, man-made buildings were individually identified from high-resolution satellite imagery. Census data was then proportionally mapped onto the urban areas of interest. This proportional allocation was validated by comparison with census data for Ghana, Malawi, and Vietnam. In addition, the predicted spatial incidence of Malawi household residences was verified using household survey data, although no attempt was made to estimate the number of occupants in specific residences.

A novel approach to estimating the urban footprint
Deville [47] and his collaborators analyzed two large datasets of mobile phone (MP) calls made in France and Portugal during 2007 and 2008. Using the MP data, MP usage was estimated as a function of time of day, season, and location. Usage was scaled by the census-estimated national population to estimate the local population densities. The estimated population densities compared favorably to estimates derived independently from multiple sources.

Accuracy and image resolution
In a comprehensive review of low resolution (LR) and high resolution (HR) mapping techniques, [31] found broad agreement between the corresponding maps. However, HR methods including GUF and GHSL were superior with respect to "completeness and precision, " including more accurate representation of the boundaries of large urban areas where the density of residential structures are decreasing. In rural sub-Sahara Africa, these HR methods were the only systems that successfully classified at least some of the small, disconnected settlements in rural areas.

The challenge of improvised settlements
Unplanned improvised settlements [3], temporary refugee camps [1] and other informal communities pose special challenges to almost any population mapping method. First, the clusters of improvised structures may be difficult to delineate and extract from aerial imagery, although our own methodology does not utilize feature extraction. Second, the presumed relationship between rooftop area and occupancy levels may be violated; this can make proportional allocation of the population more difficult, a method used in several of the preceding studies [29,47]. As an example of the latter, a specific residence in Bo was included that had a measured rooftop area of 8 meters on a side, or about 25 feet × 25 feet. Seven families, for a total of 62 inhabitants, used this residence as sleeping quarters. This serves as an example of extreme resource-poverty, where a single family may be forced to live in one room, with multiple small rooms within a single residence [Rashid Ansumana, personal communication.] Hillson et al. Int J Health Geogr (2019) 18:16 Given k covariates in each candidate regression model M j , there will be 2 k models in the model space. For values of k exceeding 14, an exhaustive enumeration of candidate regression models is computationally challenging, and a heuristic algorithm is usually required. Our simulations were run using the birth-death MCMC sampler, and the random number generator seed was explicitly specified in order to permit replicability. The sampler continually draws random models, and moves to the latest model if its marginal likelihood exceeds that of the current model. Let X be a specific subset of k covariates, y the fixed set of N observations. The number of times each candidate regression model is retained converges to the distribution of the Posterior Model Probability p(M i |y, X) [19].
At the conclusion of the stochastic simulation, the single best model can be identified and fitted using the point estimates of the specified coefficients as a linear multivariate regression model. This simplifies the use of scale transforms for the dependent variable (log, inverse, square root) and their associated back-transforms; the back-transformed regression model can then be crossvalidated by using the k-1 "Leave One Out" Cross-Validation (LOOCV) [21] paradigm that will be discussed subsequently.

Selection of β priors
As will any Bayesian regression model, the choice of the priors for the β vector is critical, because the choice of priors determine model selection and shrinkage. In the BMS software, the Bayesian normal-conjugate linear model defines the multiple regression model, and "Zellner's g-prior" defines the priors for the β s (i.e. the regression coefficients) [19,48]. We set the hyperparameter g 1 equal to the Unit Information Prior (UIP), the default value, which proved to be effective. Because the number of candidate regressors greatly exceeded the maximum sample size of 20, the prior for the model size was set equal to k = 20/2 = 10 . See [34] for a discussion of the rationale for the UIP prior and others g-priors, and [17] pp. 16-18 for additional discussion on the choice of priors.

Determination of best single model
At the conclusion of the stochastic simulation, the single best model (the "top model" [17]) can also be exported as a linear model estimated under Zelner's g-prior, and used for subsequent analysis. This model will yield a limited number of significant regression coefficients that can be analyzed further using conventional OLS analysis. The significance of the individual covariates is quantified by the PIP, which is equivalent to the Bayesian information criterion (BIC) developed by Schwarz [26]. A dynamic list of the best 1000 regression models found by the MCMC sampler was maintained, and this permitted the analysis of convergence during the iterative process. The number of models tracked is user-specified.