 Research
 Open Access
 Published:
Estimating the size of urban populations using Landsat images: a case study of Bo, Sierra Leone, West Africa
International Journal of Health Geographics volume 18, Article number: 16 (2019)
Abstract
Background
This is the third paper in a 3paper 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 multiband 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 crossvalidation analysis, suggest that masking nonurban areas in the Landsat section images prior to computing the candidate covariate regressors should further improve model generality.
Introduction
In resourcelimited 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 lowermiddleincome countries.
Such estimates are invaluable for health planning, refugee support [1], epidemiological modeling [2], and for state and municipalitysponsored 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 of interest. In our current study, we examine the possibility of using Landsat 5 thematic mapper (TM) data to estimate the population densities of sections in Bo, Sierra Leone, without the necessity of either explicitly estimating the number of individual residential structures present nor a requirement to extract and estimate rooftop areas.
Description of the study area
Bo is Sierra Leone’s second largest city, and its population and footprint has grown substantially over the past two decades. The city of Bo itself is approximately 30.10 km^{2} in area, and is divided into 68 mutuallyexclusive neighborhoods or sections [2]. These sections vary in size from 0.02 to 2.33 km^{2}. For 20 of the 68 sections, residential survey data collected in 2011 are available [3] as summarized in Table 1.
Our primary objective is to construct models for estimating the population densities \({\hat{d}}_{i = 1,\ldots ,20}\), from which estimates of the section populations \({\hat{p}}_{i = 1,\ldots ,20}\) will be derived. Fig 1 shows the 20 surveyed sections ordered by population density.
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), pixellevel 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 \(\langle \,d,Area\rangle \) of the estimated population densities \({\widehat{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.
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 bandspecific pixellevel 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 \(k1\) crossvalidation method (also referred to as “Leave one out crossvalidation,” or LOOCV), how effectively do these regression models generalize to estimating the population density of a section deliberately omitted from the LOOCV training set?
Methods and materials
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 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 openaccess literature.
The Landsat thematic mapper (TM)
Landsat 5 was an Earthobserving 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 opticalmechanical “whisk broom” (alongtrack) scanner [6, 7]. The scanner’s mirror system bidirectionally 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 northsouth by 183 km eastwest (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 LEDAPSprocessed dataset is available for the desired imagery [9]. The 3 major steps in LEDAPS processing are:

1
As a function of the bandspecific sensor gain and bias, convert the Landsat sensor outputs to sensor spectral radiances, the energy reaching the sensors.

2
As a function of the earthsun 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 wavelengthspecific 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.
Pixellevel 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, pixellevel 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 minmax 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 minmax pixel distributions defined the normalized variables nb [mean value of normalized LEDAPScorrected 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 nonzero 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 lowdimensional 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 is a color encoded map showing the population densities of the 20 surveyed sections.
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 loworder 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 LEDAPSprocessed imagery provided by USGS. We also used Bayesian MCMC (Markov chain Monte Carlo) sampling to find the variables for our regression models, rather than stepwise 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 LEDAPScorrected pixellevel intensity measurements. Table 5 summarizes the 3 datasets used in subsequent analysis: (1) nonspectral transforms, (2) spectral transforms, and (3) the total combined dataset. There are 379 total variables, with a subset of 304 spectral transforms and 75 nonspectral 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=\frac{Persons_i}{Area_{i = 1,\ldots 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.
Regression models
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. Second, if the Pearson correlation between a pair of covariates was .99 or greater, one of the covariates was dropped.

2.
Data transformation. Different candidate transforms for the dependent variable \(d_{i}=persons_i/km^2\) were evaluated to improve the linearity of the regressive estimator for \({\hat{d}}\). The square root transform \(\sqrt{d}\) was selected as the dependent variable to be estimated.

3.
Regression analysis. A Bayesian mixture analysis was run, using an MCMC (Markov chain Monte Carlo) MetropolisHastings sampler to evaluate the candidate regression equations [17,18,19]. A brief summary of the methods used is provided in Appendix 3. The best single equation found for estimating \(\widehat{\sqrt{d}}\) during the stochastic sampling was transformed to a conventional linear multiple regression equation.

4.
Backtransform \(\widehat{\sqrt{d}}\). The transformed estimated population density vector \(\widehat{\sqrt{d}}\) was backtransformed [20] into the original parameter space as \(\widehat{d_i}\). The goodnessoffit of the regression equation for estimating \({\widehat{d}}\) could then be evaluated. The population of each section was also estimated.

5.
Crossvalidation. “Leaveout one crossvalidation” (LOOCV) [21] was used to quantify how well the regression equation generalizes to estimating observations that were not included in the training set.
Results
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 trialanderror we found that all of the covariates selected were from the subset of Landsat variables defined for spectral (i.e. interpixel) 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.
Data transformation
Figure 3 shows the backtransformed estimated population density for (\(\hat{d_i}\) = persons\(_i\)/km\(^2\)), plotted as a function of the section population density for each transform of d. The regression model used was the top model in an ordered mixture of the 1000 bestfitting regressions found in the MCMC sample space. The green line is the true value of d. No transform was applied in plot (A), (B) is the backtransformed log transform (i.e. \(e^{\widehat{ln(d)}}\)), and (C) is the backtransformed square root transform (i.e. \([{\widehat{\sqrt{d}}}]^2\)). The square root transform \({\widehat{\sqrt{d}}}\) yielded the most linear estimation of the population density.
Regression analysis
Table 8 gives the parameters for the best regression model found for estimating \(\sqrt{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.
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.
Estimating section areas
Figure 4a shows the backtransformed estimates of the populations densities \({\widehat{d}}_i\), plotted as a function of the measured population densities. The regression equation in Table 8 was used to estimate \(\widehat{\sqrt{d}}\). The vector of estimates, and their .95 confidence intervals, were both backtransformed into the original parameter space: \({\hat{d}}=[\widehat{\sqrt{d}}]^2\) [20]. Panel (B) shows the estimate of the population obtained by multiplying the backtransformed estimate of \({\hat{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 \(\widehat{\sqrt{d}}\) and the estimated backtransformed population density \({\hat{d}}={[\widehat{\sqrt{(}d)}}]^2\). The RE can be positive or negative, and the Mean RE is the mean of the absolute values of RE.
Table 9 lists the REs by section, as well as \(\hat{d_i}\) and \(\hat{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 crossvalidation
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 crossvalidation. As an alternative, we used k1 crossvalidation, 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 (LeaveOut One CrossValidation) 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 \(n1\) 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 \(\sqrt{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 crossvalidation tests.
Although the median absolute value of RE for the backtransformed estimate is only 11.14%, the model failed to generalize (i.e. crossvalidate) 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 nonresidential 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 crossvalidation 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 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 finescale 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 bestfitting models that were retained by the sampler, as well as in the “top” model. (The number of top models tracked by MCMC sampler is userselectable.) 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 bestfitting 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 triband 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 manmade 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 nonhomogeneous 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 nonresidential buildings within an area may confound residential and nonresidential 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 nonresidential 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.
Landuse/landcover models
(LU/LC) modeling offers another approach to resolving the consequences of nonhomogeneous land use. A LU/LC model would differentiate between categories of land use in different areas of a section, distinguish between residential and nonresidential 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 groundtruth 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 highspeed 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 groundtruth 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 500\(m^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 crossvalidation decreases. It would also be possible to define training sets and test sets for conventional crossvalidation testing.
Even given highresolution 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 subsampled 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 nonhomogeneity of the population density within certain sections by masking (i.e. omitting) nonresidential subareas 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=(\text {Band 4}  \text {Band 3})/(\text {Band 4} + \text {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 crossvalidation error by improving compliance with the model’s assumptions.
Alternative approaches to crossvalidation
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.483475.00)/3475.00 \(\times \) 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 lowerthanobserved population and others had a higherthanobserved 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 (k1) crossvalidation. The results here were less satisfactory than for the population density \({\hat{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 crossvalidation 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 crossvalidation.
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 highlycorrelated variables executed first.
In this case, a 6 variate regressor equation was found, plus the nonzero 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 crossvalidation 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 crossvalidation. About half of the footprint for Reservation is bright green wetlands, and Moibawo Farm is heavily forested. The crossvalidation 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 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 nonhomogeneity of population density in several sections, including Reservation and Moibawo Farm, we have succeeded in this objective. An accurate 6covariate 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 highresolution 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) preprocessed 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 highresolution (<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 informationtheoretic 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.
Availability of data and materials
All data are fully available without restriction, with the relevant tabular data within the paper and its Appendices. GIS data are available on OpenStreetMap (http://osm.org/go/am_ZKeeU). Landsat imagery are available from the US Geological Survey (USGS).
Notes
 1.
A hyperparameter is a parameter on the \(\beta \) distribution, rather than a regression model parameter.
Abbreviations
 BIC:

Bayesian information criterion
 CV:

Coefficient of variation
 DOF:

Degrees of freedom
 EBL:

Empirical local Bayes estimator
 GIS:

Geographic information system
 LEDAPS:

Landsat Ecosystem Disturbance Adaptive Processing System
 LISA:

Local Indicators of Spatial Association
 LOOCV:

Leave one out crossvalidation
 LU/LC:

Landuse/landcover
 MCMC:

Markov chain Monte Carlo
 NDVI:

Normalized Difference Vegetation Index
 NIR:

Near infrared
 PCA:

Principal components analysis
 PIP:

Posterior inclusion probability
 R, G, B:

Red, green, blue
 SLC:

Scan line corrector
 TM:

Thematic mapper
 TOA:

Top of the atmosphere
 USGS:

U.S. Geological Survey
 VIF:

Variance inflation factor
 %RE:

Relative proportional error
References
 1.
Checchi F, Stewart BT, Palmer JJ, Grundy C. Validity and feasibility of a satellite imagerybased method for rapid estimation of displaced populations. Int J Health Geogr. 2013;12(4):12.
 2.
Ansumana R, Malanoski AP, Bockarie AS, Sundufu AJ, Jimmy DH, Bangura U, Jacobsen KH, Lin B, Stenger DA. Enabling methods for community health mapping in developing countries. Int J Health Geogr. 2010;9(1):56.
 3.
Hillson R, Alejandre JD, Jacobsen KH, Ansumana R, Bockarie AS, Bangura U, Lamin JM, Malanoski AP, Stenger DA. Methods for determining the uncertainty of population estimates derived from satellite imagery and limited survey data: a case study of Bo City, Sierra Leone. PLoS One. 2014;9(11):112241. https://doi.org/10.1371/journal.pone.0112241.
 4.
Hillson R, Alejandre JD, Jacobsen KH, Ansumana R, Bockarie AS, Bangura U, Lamin JM, Stenger DA. Stratified sampling of neighborhood sections for population estimation: a case study of bo city, sierra leone. PLoS One. 2015;10(7):0132850. https://doi.org/10.1371/journal.pone.0132850.
 5.
Harvey J. Population estimation models based on individual TM pixels. Photogramm Eng Remote Sens. 2002;68:1181–92.
 6.
Mika AM. Three decades of Landsat instruments. Photogramm Eng Remote Sens. 1997;63(7):839–52.
 7.
USGS: Landsat Thematic Mapper (TM). 2015. https://lta.cr.usgs.gov/TM.
 8.
NASA: Landsat 5. 2017. https://landsat.gsfc.nasa.gov/landsat52/. Accessed 13 Feb 2017.
 9.
USGS: LT52010542011001MPS00 Landsat TM Imagery. 2015. http://earthexplorer.usgs.gov.
 10.
USGS: Landsat Data Access. 2016. https://landsat.usgs.gov/landsatdataaccess.
 11.
USGS: Landsat Collections. 2016. https://landsat.usgs.gov/landsatcollections.
 12.
USGS: Landsat 7. 2016. https://landsat.usgs.gov/landsat7.
 13.
Masek JG, Vermote EF, Saleous N, Wolfe R, Hall FG, Huemmrich KF, Gao F, Kutler J, Lim TK. LEDAPS calibration, reflectance, atmospheric correction preprocessing code, version 2. ORNL Distributed Active Archive Center (2013). https://doi.org/10.3334/ORNLDAAC/1146.
 14.
USGS: What are the best spectral bands to use for my study?. 2016. https://landsat.usgs.gov/whatarebestspectralbandsusemystudy.
 15.
Wende, C.: An introductory Landsat tutorial. Technical report. 2004. http://www.ricercasit.it/public/documenti/clamSiTel/Materiali/Moduli%20Didattici%20I%20anno/Telerilevamento/esercitazione/Landsat_GeoCover_Tutorial.pdf.
 16.
SuarezAlvarez MM, Pham DT, Prostov MY, Prostov YI. Statistical approach to normalization of feature vectors and clustering of mixed datasets. Proc Math Phys Eng Sci. 2012;468(2145):2630–51. https://doi.org/10.1098/rspa.2011.0704.
 17.
Zeugner S, Feldkircher M. Bayesian model averaging employing fixed and flexible priors: the BMS package for R. J Stat Softw Articles. 2015;68(4):1–37. https://doi.org/10.18637/jss.v068.i04.
 18.
Hoeting JA, Madigan D, Raftery AE, Volinsky CT. Bayesian model averaging: a tutorial (with comments by M. Clyde, David Draper and E.I. George, and a rejoinder by the authors. Statist Sci. 1999;14(4):382–417. https://doi.org/10.1214/ss/1009212519.
 19.
Amini S, Parmeter CF. Bayesian model averaging in R. JESM. 2011;36(4):253–87. https://doi.org/10.3233/JEM20110350.
 20.
McDonald JH. Handbook of biological statistics. 3rd ed. Baltimore: Sparky House Publishing; 2014.
 21.
James G, Witten D, Hastie T, Tibshirani R. An introduction to statistical learning: with applications in R. New York: Springer; 2014.
 22.
Naimi N, Hamm NAS, Groen TA, Skidmore AK, Toxopeus AG. Where is positional uncertainty a problem for species distribution modelling? Ecography. 2014;37(2):191–203. https://doi.org/10.1111/j.16000587.2013.00205.x.
 23.
Collins R, Hanson A, Riseman E, Schultz H. Automatic extraction of buildings and terrain from aerial images. In: Gruen A, Kuebler O, Agouris P. (eds.) Automatic extraction of ManMade objects from aerial and space images. Monte Verità. Basel: Birkhäuser; 1995. p. 169–178. https://doi.org/10.1007/9783034892421_16
 24.
Blaschke T. Object based image analysis for remote sensing (review article). ISPRS J Photogramm Remote Sens. 2010;65:2–16.
 25.
Doppelhofer G, Miller RI, SalaiMartin X. Determinants of longterm growth: a Bayesian averaging of classical estimates (BACE) approach. Working Paper 7750, National Bureau of Economic Research (2000). https://doi.org/10.3386/w7750.
 26.
Schwarz G. Estimating the dimension of a model. Ann Statist. 1978;6:461–4. https://doi.org/10.1214/aos/1176344136.
 27.
Wilson C. Spectral analysis of civil conflictinduced forced migration on landuse/landcover change: the case of a primate and lowerranked cities in sierra leone. Int J Remote Sens. 2014;35(3):1094–125. https://doi.org/10.1080/01431161.2013.875633.
 28.
Getis A, Ord JK. The analysis of spatial association by use of distance statistics. In: Anselin L, Rey SJ, editors. Perspectives on spatial data analysis. New York: Springer; 2010. p. 127–45.
 29.
Tiecke TG, Liu X, Zhang A, Gros A, Li N, Yetman G, Kilic T, Murray S, Blankespoor B, Prydz EB, Dang H.AH. Mapping the world population one building at a time. 2017. arXiv:1712.05839v1.
 30.
Anselin L, Syabri I, Kho Y. Geoda: an introduction to spatial data analysis. Geogr Anal. 2006;38(1):5–22.
 31.
Klotz M, Kemper T, Geiß C, Esch T, Taubenböck H. How good is the map? a multiscale crosscomparison framework for global settlement layers: evidence from central europe. Remote Sens Environ. 2016;178:191–212.
 32.
Spiker JS, Warner TA. Scale and spatial autocorrelation from a remote sensing perspective. In: Jensen RR, Gatrell JD, McLean D, editors. Geospatial technologies in urban environments. New York: Springer; 2007. p. 197–213.
 33.
Weier J, Herring D. Measuring Vegetation (NDVI & EVI). 2000. http://earthobservatory.nasa.gov/Features/MeasuringVegetation/measuring_vegetation_1.php.
 34.
Liang F, Paulo R, Molina G, Clyde MA, Berger JO. Mixtures of g priors for Bayesian variable selection. J Am Stat Assoc. 2008;103(481):410–23. https://doi.org/10.1198/016214507000001337.
 35.
Forster BC. Urban residential ground cover using Landsat digital data. Photogramm Eng Remote Sens. 1980;46(4):547–58.
 36.
Forster B. Some urban measurements from Landsat data. Photogramm Eng Remote Sens. 1983;79(12):1693–707.
 37.
Forster BC. Coeffcient of variation as a measure of urban spatial attributes, using SPOT HRV and landsat TM data. Int J Remote Sens. 1993;14:2403–9.
 38.
Harvey J. Small area population estimation using satellite imagery. Stat Transit. 2000;4:611–33.
 39.
Wardrop NA, Jochem WC, Bird TJ, Chamberlain HR, Clarke D, Kerr D, Bengtsson L, Juran S, Seaman V, Tatem AJ. Spatially disaggregated population estimates in the absence of national population and housing census data. 2018;115:3529–37. https://doi.org/10.1073/pnas.1715305115.
 40.
Patel NN, Angiuli E, Gamba P, Gaughan A, Lisini G, Stevens FR, Tatem AJ, Trianni G. Multitemporal settlement and population mapping from Landsat using Google Earth Engine. Int J Appl Earth Obs Geoinform. 2015;35:199–208.
 41.
Trianni G, Lisini G, Angiuli E, Moreno EA, Dondi P, Gaggia A, Gamba P. Scaling up to national/regional urban extent mapping using Landsat data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 2015;8(7):3710–9. https://doi.org/10.1109/JSTARS.2015.2398032.
 42.
Son NT, Chen CF, Chen CR, Chiang SH. Mapping urban growth of the capital city of Honduras from Landsat data using the impervious surface fraction algorithm. Geocarto Int. 2016;31(3):328–41. https://doi.org/10.1080/10106049.2015.1047469.
 43.
Linard C, Kabaria CW, Gilbert M, Tatem AJ, Gaughan AE, Stevens FR, Sorichetta A, Noor AM, Snow RW. Modelling changing population distributions: an example of the Kenyan Coast ,1979–2009. Int J Digit Earth. 2017;10(10):1017–29. https://doi.org/10.1080/17538947.2016.1275829.
 44.
Esch T, Heldens W, Hirner A, Keil M, Marconcini M, Roth A, Zeidler J, Dech S, Strano E. Breaking new ground in mapping human settlements from spacethe global urban footprint. ISPRS J Photogramm Remote Sens. 2017;134:30–42.
 45.
Pesaresi M, Ehrlich D, Florczyk AJ, Freire S, Julea A, Kemper T, Syrris V. The global human settlement layer from Landsat imagery. In: Procedings IEEE Internatinal Geoscience and Remote Sensing Symposium (IGARSS); 2016 p. 7276–7279. https://doi.org/10.1109/IGARSS.2016.7730897.
 46.
Tiecke T. Open population datasets and open challenges (2016). https://research.fb.com/openpopulationdatasetsandopenchallenges/.
 47.
Deville P, Linard C, Martin S, Gilbert M, Stevens FR, Gaughan AE, Blondel VD, Tatem AJ. Dynamic population mapping using mobile phone data. Proc Natl Acad Sci. 2014;111(45):15888–93.
 48.
Feldkircher M, Zeugner S, and: Benchmark priors revisited:on adaptive shrinkage and the supermodel effect in Bayesian model averaging. IMF Working Paper 09, 1. 2009. https://doi.org/10.5089/9781451873498.001.
 49.
NASA: Make Your Own LandsatImage Tutorial. 2013. https://landsat.gsfc.nasa.gov/wpcontent/uploads/2013/05/MakeYourOwnLandsatImageTutorial.pdf.
Acknowledgements
Not applicable.
Funding
Funded by the Defense Threat Reduction Agency, Joint Science and Technology Office (http://www.dvidshub.net/unit/DTRACB#.UoUqZ9wo5zk ) via contract to the U.S. Naval Research Laboratory. There is no past, present or future Intellectual Property associated with the work described in the paper, and none of the authors have any financial interests or conflicts in the outcome of the study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Author information
Affiliations
Contributions
RH designed the study, wrote the manuscript, and developed the models. AC developed models and their accompanying portions of the manuscript. JDA, KHJ, RA, ASB, UB, JML, and DAS contributed to the acquisition and interpretation of data and to revising the manuscript. All authors read and approved the final manuscript.
Corresponding author
Correspondence to David A. Stenger.
Ethics declarations
Ethics approval and consent to participate
All data collection involving human subjects was approved, including for this study, by a total of three independent Human Subjects Research Institutional Review Boards: Njala University, George Mason University, and the U.S. Naval Research Laboratory. Written informed consent was obtained from each household representative who participated in the survey. Survey data were obtained as part of a broader study to determine not only population demographics but health metrics and health care utilization trends.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Roger Hillson: Independent consultant; retired from the Information Technology Division, Naval Research Laboratory, Washington, DC, USA.
Appendices
Appendices
Appendix 1: an overview of methods for satellitebased 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 nonresidential areas [35]. Applying Landsat data, Forster also developed a set of pixelbased 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 singleband 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 “bottomup” population estimation models using highresolution satellite imagery.
Satellite imagery has also been used to develop models to estimate the spatial footprint of cities in low and/or middleincome 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 multiyear 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 tenyear 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 highresolution satellite measurements using SAR (Synthetic Aperture Radar) imaging [44]. The highly automated GUF (Global Urban Footprint) utilizes SAR amplitude data from the TerraSARX and TanDEMX 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 flatroof houses may produce false negatives. An unsupervised classification system, based on Support Vector Data Description (SVDD), used these paired distributions to identify builtup 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 40year period.
HRSL
The HRSL (High Resolution Settlement Layer) [29, 46] effort is a fourth project for mapping the urban footprint. In this effort, manmade buildings were individually identified from highresolution 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 censusestimated 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 subSahara 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 \(\times \) 25 feet. Seven families, for a total of 62 inhabitants, used this residence as sleeping quarters. This serves as an example of extreme resourcepoverty, where a single family may be forced to live in one room, with multiple small rooms within a single residence [Rashid Ansumana, personal communication.]
Appendix 2: variable definitions
This Appendix contains the definitions for all of the candidate covariates used in this paper. See Table 5. There are a number of differences between the dataset in Table 6 and the list of candidate covariates summarized in [5]. For the normalized variable \(nb_i\), we used a minmax normalization, rather than a zscore transform. Based on preliminary simulations, the minmax transform is equally or more effective as a covariate than a zscore transforms. (The latter was calculated as the mean zscore of the pixels in a section, where the individual pixel zscores were calculated relative to the grand mean and \(\sigma \) for all sections combined.) A second difference was the introduction of a modified notation to explicitly differentiate between the nonspectral ratios of means \(r.sp_{ij}\) and the spectral ratio operator \(r_{ij}\).
Nonspectral variables
The nonspectral variables consist of the means, standard deviations, and variances of \(b_{i}\) within a single section, and the specified ratios of these values [Table 3].
The mean band values \(b_i\) and textures:
The square of \(b_{ij}\):
The nonspectral crossproduct of \(b_i\) and \(b_j\):
The nonspectral ratios of \(b_i\) and \(b_j\):
The difference/sum of \(b_i\) and \(b_j\):
Spectral transforms
The nonspectral transforms are the mean values \(b_i\) and associated measures of variation [standard deviation, variance, and coefficient of variation] of the elementary pixel measurements (\(b_i\)); as well as the squares \(s_i\), products \(p_{ij}\), ratios \(r.re_{ij}\) and sums/differences \(d_{ij}\) of the mean values \(b_i\).
The spectral transforms, which we will define next, are pixellevel functions. Unlike the nonspectral transforms, the spectral transforms are functions of the ratios, products, or differences of individual pixels, rather than functions of the mean values of band measurements.
Let p denote the number of pixels in section i, and let ij denote a specific pixel j in band i. Let \(Max_i\) be the maximum value of \(b_{ij}\) in section i, and \(Min_i\) be the minimum value of \(b_{ij}\). \(nb_{i}\) is defined as the mean value of the minmax normalized values of section i. Specifically:
From the preceding equation, each of the p normalized values in the summation term is a number between 0 and 1. The range of \(nb_{i}\), which is the mean of these values, will necessarily fall between 0 and 1. The standard deviation \(nbs_i\), variance \(nbv_i\), and CV \(nbc_i\) of \(nb_i\) are defined as:
For a given section, the spectral function \(r_{ij}\) is the mean ratio of the paired pixel values in bands i and j. Each is assumed to have p pixels (TM measurements). The variable \(bpx_{i_n}\) is the value of pixel n in band i, and \(bpx_{j_n}\) is the value of pixel n in band j.
The standard deviation of \(r_{ij}\) is denoted by \(rs_{ij}\), the variance by \(rv_{ij}\), and the coefficient of variation by \(rc_{ij}\). The variable \(ds_{ij}\), the mean of the paired ratios of the differences and sums of the pixel magnitudes in bands i and j, is computed similarly to \(rs_{ij}\).
Cylindrical coordinate transforms
The cylindrical coordinate transforms are denoted by \(ch_{ijk}\) and the rectangular coordinate transforms by \(rh_{ijk}\). In either case, 3 of the 7 available bands were assigned to the red (i), green (j), and blue (k) channels. The 3 normalized mean band values are then represented as the coordinates in a Cartesian cube. Via projection onto the chromatic plane, each color triplet is transformed into a single hexagonal number; this scalar value is the candidate covariate for the cylindrical coordinate transform.
The transform \(ch_{ijk}\) is calculated using Eqs. (12) and (13) in Hanbury [2008]. The rectangular coordinate transforms \(rh_{ijk}\) is calculated using the piecewise hue definition in Eq. (4), reference [Hanbury and Serra, 2003]. Prior to either transform, the data is normalized so that the maximum value of the data maps onto the red channel i was 255.
Appendix 3: Bayesian mixture analysis
The Bayesian Model Selection (BMS) R library (package) was used to develop and run the regression models [17,18,19]. BMS implements a Bayesian Model Averaging (BMA), based on the premise that a properly weighted mixture of Bayesian models can estimate a desired set of observations more accurately than a single regression model can. The number of best models for which the model parameters are saved is userdefined, and were conservatively set to 1000 in our simulations, rather than to the default value of 200. The sampler can be set to retain data for the best regression model only, but it would no longer be possible to graphically verify the convergence of the posterior model probabilities. or to obtain predictions based on an average mixture of models.
Define a regression model \(M_j\) for estimating the dependent variable y, \(\beta _{0}\) is the constant intercept term, and \(x_1\), \(x_2\), ...\(x_k\) are the k covariates. See [19], Eq (1).
Let \(\beta \) be the vector of parameters (regression coefficients) from model \(M_j\), and \(Pr(M_j)\) be the probability that \(M_j\) is the true model. The estimated posterior means of \({\widehat{\beta }}=(\widehat{\beta _0},\widehat{\beta _1},...,\widehat{\beta _k})\) are:
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 birthdeath 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 backtransforms; the backtransformed regression model can then be crossvalidated by using the k1 “Leave One Out” CrossValidation (LOOCV) [21] paradigm that will be discussed subsequently.
Selection of \(\beta \) priors
As will any Bayesian regression model, the choice of the priors for the \(\beta \) vector is critical, because the choice of priors determine model selection and shrinkage. In the BMS software, the Bayesian normalconjugate linear model defines the multiple regression model, and “Zellner’s gprior” defines the priors for the \(\beta \)s (i.e. the regression coefficients) [19, 48]. We set the hyperparameter g^{Footnote 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 gpriors, 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 gprior, 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 userspecified.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Hillson, R., Coates, A., Alejandre, J.D. et al. Estimating the size of urban populations using Landsat images: a case study of Bo, Sierra Leone, West Africa. Int J Health Geogr 18, 16 (2019) doi:10.1186/s1294201901801
Received
Accepted
Published
DOI
Keywords
 Remote Sensing
 Population Estimation
 Spatial Epidemiology
 Landsat
 Sierra Leone
 MCMC
 urban footprint