Population data sources

Our main sources of population density at the village level are the Census of Population and Housing, Sri Lanka, 2012, and the 2012-13 Sri Lankan Household Income and Expenditure Survey (HIES). The census universally covers all households in all of Sri Lanka’s 14, 000 villages. The HIES covers about 25,000 households in 2, 421 villages, spanning all the districts and sub-district divisions in the country. While the census is conducted decennially, the HIES is typically conducted once in three years.

The HIES follows a two-stage stratification process: the census blocks form the Primary Sampling Units (PSUs), and the households form the Secondary Sampling Units (SSUs) or Final Sampling Units (FSUs). The PSUs used in the HIES are census blocks, which are portions of villages. There were 62,571 PSUs versus 13,984 villages in the 2012 census. Therefore, the HIES can yield direct density estimates only at the PSU-level, and not at the village-level. Unfortunately, it is not possible to model population density at the PSU-level due to the lack of a PSU-level Geographical Information System (GIS) boundaries.

We therefore indirectly estimate the village population density using the HIES survey weights under the assumption that the survey weights reflect the inverse probability of the housing units being selected into the sample. This is approximately equivalent to scaling up estimated PSU population derived from the survey using that PSU’s share of the village population in the census. The steps used to estimate village level population density from the HIES are described in S1 Appendix. In order to apply this method, we limit the HIES sample to the 97 percent of villages that contain exactly one surveyed PSU in the village. The correlation in the HIES sample between the estimates of population density derived from the survey, and the actual village population density obtained from the census, is 0.91. This strong correlation indicates that this method, when applied in our context, generates credible village-level population density estimates.

Geo-spatial data

The satellite imagery based geo-spatial data we use are categorized into four types, based on their resolution, availability, and the nature of the extracted data. Most indicators are obtained around 2012 (with a gap of +/- 2 years), the year around which the models are estimated. First, we use low-resolution open-source indicators from three sources: (1) Night time lights from the Visible Infrared Imaging Radiometer Suite (VIIRS) at a resolution of 750 m per pixel. We use the maximum and mean intensity of two months, namely, March and September, 2014; (2) Global Forest Change data based on [16], from which we use mean tree cover in 2000, and gain and loss in forest area between 2000 and 2014, at a resolution of 30 meters per pixel; (3) Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)’s elevation and slope data at a resolution of 30 meters per pixel in 2009. Although elevation and slope data are slightly older, these measures do not typically change over time, and hence it might be still appropriate to use them.

Second, we use open-source built-up area measures based on higher resolution imagery. The Global Urban Footprint (GUF) (year 2012) and Global Urban Footprint plus (GUF+) (year 2015) provide built-up estimates at a resolution of 12 meters per pixel. While GUF was constructed using the TerraSAR-X/TanDEM-X satellites that are based on radar-detection, GUF+ additionally incorporates data from the Sentinel and Landsat satellites that are based on optical remotely sensing. We also use Global Human Settlement Layer’s (GHSL) built-up estimates that are at approximately 38 meters per pixel and Facebook’s High Resolution Settlement Layer (HRSL) built-up area measure, available at a resolution of 30 meters per pixel.

Additionally, we use geo-spatial indicators extracted by [13] from very-high resolution commercially available satellite imagery acquired from the company, Maxar (formally DigitalGlobe). These indicators cover about 1,178 to 1,275 villages within 55 sub-districts in the country. The list of these 55 sub-districts is available in S2 Appendix. These imagery were mostly acquired for the years 2011 and 2012, although some were also captured for the year 2010. Our third set of indicators include texture features extracted by [13], including Fourier transform, Gabor filter, Histogram of Oriented Gradients (HoG), Line Support Regions (LSR), Pantex, Normalized Difference Vegetation Index (NDVI), and Speed-Up Robustness Features (SURF). Fourth, we use object-based features extracted by [13], including the number of cars, building count and size, roof type, shadow pixels, road length and type, and type of farm land (paddy or plantations). The features are identified using deep learning-based convolutional neural networks (CNN) and object based image processing software. More details on these features can be found in [13]. The data purchased and indicators derived by [13] costed $125,000 in the year 2015, but these data are considerably cheaper now.

The other geo-spatial variables we employ include the logarithm of village-area, and locational indicators such as binary indicator variables for the districts and the sector. There are 25 districts in Sri Lanka. The official sector classification includes the urban, rural, and estate sectors. Although 91% of the villages have a unique sectoral assignment, there are multiple sectors associated with each village in the remaining 9%. To simplify, we first create the following modified sectoral definition: (1) rural or estate; (2) urban, and then assign villages with multiple assignments to the urban sector.

Table 1 concisely presents the key variables used in the study, along with their source, time frame, resolution type, and geographical coverage. Note that all data sources are Tables 2 and 3 present the summary statistics for the population variables and geo-spatial indicators for the national sample and the 55 sub-districts, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Satellite indicators by resolution and availability.

https://doi.org/10.1371/journal.pone.0237063.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary statistics for villages in the national sample.

https://doi.org/10.1371/journal.pone.0237063.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Summary statistics for villages in the 55 sub-districts sample.

https://doi.org/10.1371/journal.pone.0237063.t003

Modeling

We use Poisson regressions to model population density at the village level. This takes the following form: (1) where P_v is the population density (persons/km²) of village v. X_v is a vector of satellite imagery-based indicators defined for village v. Z_v is a vector that includes the natural logarithm of village area (Ln area_v), an indicator variable for urban areas, and indicator variables for Sri Lanka’s districts.

To prevent model over-fitting, we employ the least absolute shrinkage and selection operator (LASSO) regularization which estimates a regression model with an added constraint that enforces parsimony. We follow a two-step procedure. First, a full set of variables is included, on which LASSO regularization is conducted to choose variables. Second, the final Poisson model is estimated using the chosen variables. Eq 2 presents the objective function for LASSO regularization. Let τ represent the vector of all coefficients, β and γ, put together, of size p. (2) where l is the value of the log likelihood function of the Poisson model using parameters α and τ, and λ is a non-negative regularization parameter. While setting λ = 0 yields unconstrained Poisson regression estimates, a large λ penalizes the absolute values of the elements of τ. Fivefold cross validation is applied to choose the value of λ that minimizes the root-mean squared error (RMSE) across the folds. The in-sample R² from this step indicates the goodness of the model fit. Since LASSO also generally shrinks the magnitude of all coefficients towards zero [17], we avoid biased predictors by running a simple Poisson model based on the Lasso selected variables. From this, we obtain the out-of-sample R², mean absolute error (MAE), and RMSE, all of which indicate out-of-sample predictive accuracy. For estimations, we use the glmnet package in R which we invoke from Stata, using the wrapper, rlasso [18].

First, we implement the above model using the village level census population density as the dependent variable (in the census year, 2012) to probe how well geo-spatial data can predict population density. We conduct this exercise for a model in the entire country separately using different indicators: (i) district and urban-area binary indicators; (ii) including log village-area to i; (iii) including low-resolution open-source indicators to ii; (iv) including higher resolution open-source indicators to iii. Doing this enables us to compare prediction accuracy between higher versus lower resolution data in the national sample. In the 55 sub-districts where the features from [13] are available, we implement models using: (v) texture-features plus indicators in iv; and (vi) object-features plus the indicators used in v. For the 55 sub-districts, we compare prediction accuracy successively from models in (i) through (vi). We also repeat the analysis separately for urban areas, and for rural/estate areas.

We conduct two robustness checks. Focusing on the 55 sub-districts, we reduce the training sample size to one-half and to one-quarter of the original sample, to examine if the prediction accuracy is sensitive to the size of the training sample. We also implement a flexible random forest algorithm (using the randomForest package in R) using population density as the response variable, and employing the different sets of geo-spatial data described above, and report out-of-sample accuracy measures.

Next, we implement the bottom-up method focusing on the 1,178 villages in the 55 sub-districts sub-sample for which commercially-available texture and object features were obtained. Out of these villages, the HIES was conducted in 414 villages. In this sub-sample, we estimate Eq 1 using the HIES density as a function of all the geo-spatial indicators, including the commercially-procured satellite indicators. Note that the dependent variable, population density, is now measured for the year 2012-2013, the year the HIES was conducted. Since this is a survey-based model, we use the inverse of village population for village v as weights (). In addition, there could be factors correlated with geo-spatial indicators that affect the probability of selection of the village into the HIES. We also, therefore, adjust the weights based on the predicted probability of the village being sampled into the HIES based on geo-spatial indicators available for all villages [19, 20]. This probability is obtained from the following probit model: (3) where IN_v is a binary indicator for whether village v is sampled into HIES, and lasso_v represents all LASSO-selected variables from the Poisson model. We estimate Eq 3 using data for all the 1,178 villages in the 55 sub-districts. The corrected weights based on the probability of the village selection into the HIES, , is given by .

Using the estimated model, we predict out-of-sample densities for the remaining 946 non-HIES villages and report out-of-sample accuracy measures by comparing them with the actual census density. To provide additional context, we compare the accuracy of this model with the accuracy of estimates from a similar model that uses census-density as the dependent variable, and with the accuracy of the estimates of other bottom-up estimates in the literature, and with the accuracy of open-source top-down population products. Further, we calculate population count estimates by multiplying these density estimates with village-level area, and report their accuracy with respect to the census population.

Predictive power of Geo-spatial data

Poisson regression results using the census data based on the variables selected from LASSO regularization indicate that geo-spatial indicators have strong predictive power in predicting village level population density (Table 4). The results are presented in panel A. At both the national level and in the 55-sub-districts sample, simple locational indicators do not explain much of the variation in population density. The out-of-sample estimation R² is only about 0.36. However, adding log village area as an explanatory variable tremendously improves the out-of-sample R² to 0.65 in the national sample and to 0.589 in the 55-sub-districts. At the national level, conditional on the locational indicators and village size, the value added in using low-resolution open-source indicators (out-of-sample R² is 0.702), and additionally using high resolution open-source indicators (out-of-sample R² is 0.75) is moderate but non-trivial. In the 55 sub-districts, adding commercially-procured texture features to open-source imagery-based models does not improve the prediction accuracy (out-of-sample R² remains at around 0.7). However, adding object features to the model improves the out-of-sample R² by 10 points (0.83). For the interested reader, the marginal effects of the full set of LASSO-selected variables from the model using all geo-spatial indicators are provided in S1 Table.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Accuracy of population density estimates.

https://doi.org/10.1371/journal.pone.0237063.t004

The results separately for the rural/estate, and urban sectors are presented in panels B and C respectively in Table 4. Log village area remains a strong predictor of population density separately in both sectors. Using the full national sample, the open-source satellite indicators explain more variation in population density in the rural/estate sectors (0.808) than the urban sector (0.569). In the sample with 55 sub-districts, the predictive power of urban population density increases tremendously when adding the commercially-procured texture and object features, relative to only including open-source indicators (0.804 compared to 0.525). Adding commercially-procured indicators leads only to a limited improvement in prediction accuracy in rural areas (0.869 compared to 0.838). This result is consistent with greater uniformity in the relationship between built-up area and population density in rural areas, where simple measures of built-up area are sufficient to predict population density. In urban areas, due to the complex nature of the relationship between population density and buildings, where the latter can be both commercial and residential, advanced measures such as object classifiers and texture measures are required for better prediction.

As a robustness check, we show that Random forest (RF) models provide similar results in terms of out-of-bag R², RMSE, and MAE, overall, and across resolution-types and sector-types (S2 Table). These results collectively indicate that LASSO selection performs as well as random forest models in this setting. While random forest models are known to outperform LASSO models in some contexts [21], these conclusions are based on linear models. Our usage of Poisson regressions may negate the predictive advantage of random forest models, which are traditionally better able to account for non-linear relationships. Further, studies predicting population counts and density have specifically shown mixed results on comparative performances between random forest and LASSO-based estimates [6].

We also show that the accuracy of our estimates is robust to changing the size of the training data (S3 Table). Using a sample of 1,178 villages in the 55 sub-districts, the out-of-sample R² increases marginally to 0.843 if we reduce the training sample by half, with a mild reduction to 0.832 if we further reduce the sample size by another one-half. There is only a small increase in the MAE and RMSE with the reduced sample sizes.

The bottom-up method using survey and Geo-spatial data

While it is encouraging to see that geo-spatial data accurately predict census population density at the village level, these results do not show that surveys combined with geo-spatial data can be used to approximate population density. We now turn to examining this by testing the accuracy of a survey-based model in the 55 sub-districts sample. We estimate a Poisson model of HIES population density in the HIES-villages, and then obtain out-of-sample density estimates in the non-HIES villages.

The results from a model using all geo-spatial indicators, including open-source and commercially-procured imagery-based indicators are provided in panel A in Table 5. For the model using the inverse of village population as weights, the R², Spearman Rank Correlation (SRC), and the MAE between the density estimates and the census density are 0.79, 0.91, and 664, respectively (columns 1 to 3). If the weights are corrected based on the probability of the village being selected into the HIES, we obtain similar results: 0.78, 0.92, and 665 for the same measures. In panel B, we present the same results based on a model using only open-source satellite imagery. The out-of-sample prediction accuracies are very similar in this case (0.75, 0.92, and 663, and 0.77, 0.92, and 657 for the three measures using the two types of weights). From both panels A and B, we can see that the HIES-based model performs equally well compared to a similar model that uses census-density as the dependent variable.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Out-of-sample accuracy of HIES-based village population density estimates using open-source and commercially-procured satellite data, 55 sub-districts.

https://doi.org/10.1371/journal.pone.0237063.t005

We multiplied these density estimates with village-level area to obtain population count estimates. The accuracy measures based on population counts are reported in columns (4)-(6). The R² and SRC are generally lower for population counts compared to density estimates, but the MAEs are also lower. Further, we estimated a similar model using population count as the dependent variable, from which we directly obtained population estimates (S4 Table). Population count estimates derived indirectly from a density-based model are more accurate than those obtained directly from a population-count-based model.

The mean and the median relative errors are reported in columns 7 and 8, respectively. The relative errors are the same for population density and count estimates because the only difference between the two measures is the multiplicative factor, village-area. Using all geo-spatial indicators, the mean and median errors are 37% and 28%, and 35% and 29%, without and with the additional weight correction. Using only open-source imagery-based indicators, the mean and median REs are 34% and 28% for models without and with weights correction.

Commercially-available versus open-source imagery

Table 5 predictionaccuracy provides the model accuracy using both open-source indicators and commercially-available indicators for the 55-sub district sub-sample in panel-A, and using the open-source imagery based indicators for the national sample (all sub-districts) in panel-B. It is striking that the accuracy of the model involving open-source imagery is marginally better than the accuracy of the model including the commercial-indicators. Therefore, while commercial indicators may be expensive, unavailable, or unobtainable in some instances in specific countries, our model without necessarily relying on them, works almost equally well with open-source imagery based indicators.

We estimated the model by rural and urban areas separately. These results are not reported in the paper, but are available upon request. Mirroring Table 5, in rural areas, the predictive power of the models do not change with and without high-resolution commercial imagery. However, in urban areas, where there is a complex nature of the relationship between population density and buildings as elaborated earlier, commercial imagery does increase the predictive power much more compared to a model that only has low and medium resolution open-source imagery.

Cost effectiveness

The data purchased and indicators derived by [13] cost $125,000 in the year 2015, although these data are considerably cheaper now. Procuring satellite imagery to update population density in the between-census years would be considerably cheaper than repeating a whole census at frequent intervals. To provide a sense of how expensive a whole census is, we draw upon a study that examined the cost of conducting Ghana’s housing census [22]. Comparing with Ghana, if anything, underestimates the per person cost of census data collection in Sri Lanka because Ghana is comparatively less developed. The study estimates the cost of conducting the census in Ghana at $2.5/person. Given Sri Lanka’s total population of 20 million, the total cost of conducting the census in Sri Lanka works out at least $50 million dollars.

We make this cost comparison with the census not to indicate that our method can replace the decennial censuses, but to indicate that it cheaper to update population density estimates in the between-census years using our method, rather than to repeat full censuses more frequently. Furthermore, if we used only open-source imagery, for which the accuracy of the model estimates are similar, the modelling costs are even lower.

Regular updating of satellite imagery

As much as the survey, the model relies crucially on the availability of updated geo-spatial data. It might be worthwhile to inspect if these data and indicators are available in regular intervals. Remarkably, many of the data sets in Table 1 are being regularly updated. For instance, the GHSL and the tree cover data are being updated, and GUF data has been updated with the release of the World Settlement Footprint which is in the process of being updated again now [23–25]. Further, all the commercial data used in this study can all be updated with new imagery as it is collected. While one of the open-source datasets, the Facebook data, may not be scheduled to be updated regularly, remarkably, the prediction accuracy does not fall if we drop Facebook built-up area from the model (these results are not presented in the paper, but available upon request).

Moreover, the costs for acquiring these data are dropping as more companies are able to produce the data and more satellite providers collect the data. Possibly better data sets will be created in the future that could be used in similar approaches described in this manuscript. The Gates Foundation is funding building-footprint extraction for most of Africa from very high spatial resolution imagery and plans to provide the data to research such as this [26]. We have seen some output from this work, and used them in other studies and found them to be very accurate; these data would be useful for modeling such as the present study. Additionally, similar work is being done by Microsoft in the US, that has also released these data to the public. Some of us have been in discussions with them in expanding this to other countries. Other sources of data, such as the Open Street Map are getting better and more complete every day [27]. Some of the inputs used such as slope and elevation, while important for modeling, do not change over time. Taken together, while most of the data sets used in this study are being regularly updated, there will be plenty of data (and possibly better) to do this type of work as time progresses.

Comparison with other bottom-up estimates

The out-of-sample R²s for our bottom-up estimates range from 0.75 to 0.79 (Table 5). This is higher than the out-of-sample R² of 0.739 (the square of the reported correlation coefficient of 0.86) between original and estimated population density obtained for the case of Australia at the level of the Australian census collection districts [9]. Our median RE at 28% is lower than those reported for the Netherlands using geo-spatial data of similar dimension and resolution as ours, which ranges from 42% to 85.4% [8]. These higher errors could be attributed to the study’s usage of fewer satellite indicators, linear rather than non-linear regressions, and simple rather than stratified random sampling in the survey. Although [8] report lower errors (18.3% or 18.5%) by utilizing floor-space or volumetric (3-Dimensional) satellite data and functional information about buildings, these data are costly to obtain periodically in low-income countries.

Comparison with open-source top-down estimates

We compare our estimates with the following five “top-down” population products: (1) Landscan for the year 2010; (2) WorldPop for 2010 and 2015; (3) Facebook HRSL for 2015; (4) Gridded Population of the World v3 (year 2010) and v4 (year 2015), and (5) GHSL for 2015. All sources employ a top-down approach using a combination of areal interpolation techniques, including basic dasymetric approaches in conjunction with ancillary data, and statistical modeling methods to distribute census population to smaller grids. While GPW redistributes data based on an equal weighting technique, other sources use measures such as built-up area, road count and density, elevation, slope, and light intensity (“covariates”) to proportionally redistribute population. Facebook, WorldPop 2015, GPW 2015, and GHSL estimates are based on the 2012 Sri Lankan census; WorldPop 2010 and GPW 2010 are based on the 2001 census; and Landscan uses mid-year population of the country in the past year calculated by the Geographic Studies Branch, US Bureau of Census [28]. Detailed information on the input population, ancillary data, and the redistribution methodology for each source is presented in S5 Table.

At first glance, it may seem odd that we use built-up area measures from Facebook and GHSL, and compare our population density estimates with the same sources. It is important here to distinguish between estimates of built-up area, top-down estimates of population density, and bottom-up estimates of population density. Currently, built-up area estimates are publicly provided by many companies including the Global Human Settlement Layer and Facebook. These estimates of built-up area can be produced regularly solely from satellite imagery, and do not require a ground-truth measure of population density. Our model uses these built-up area as explanatory variables to predict population density. These companies also independently combine their built-up area measures with census data to estimate local population density. They utilize top down methods that redistribute census population counts according to built-up area, and therefore become stale over time. In contrast, our method is a “bottom-up” method that calibrates the model to survey data and therefore can be implemented in between census years. The idea behind Table 6 is to compare our estimates with other top-down estimates in a census year.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Accuracy of existing top-down population products.

https://doi.org/10.1371/journal.pone.0237063.t006

Using R², SRC, and mean REs as metrics, our model estimates are better than all top-down estimates except Facebook and Worldpop 2015 (Table 6). Our median REs are higher than only Facebook and Worldpop 2015, and comparable to GHSL and GPW 2015, and much lower compared to GPW, Worldpop, and Landscan (all in 2010). The relative performance of our population count estimates in Table 5, in comparison to the top-down population estimates in Table 6 is similar to that of the population density estimates.

Table 6 accuracy shows that the top-down products using the most recent census for redistribution, and using higher resolution imagery-based covariates are performing better than the others. For example, GPW and WorldPop 2010 estimates use the 2001 census for redistribution, and hence it is not surprising that they perform poorly in comparison to the 2012 census. Similarly, an earlier unpublished version of Facebook’s estimates we worked with only yielded a correlation coefficient of 0.5 because Facebook initially used the 2001 Sri Lankan census for calibration. The lower accuracy of products that use the older census might seem intuitive. By the same argument, it is important to foresee that even the best performing products seen in Table 6 may become inaccurate with the passage of time because their source data, the census, eventually becomes outdated. In other words, while the product that uses the 2012 census for redistribution could estimate population density accurately in 2012, it might not predict density as accurately in the later years, until the next census is available (likely in 2022) and utilized. This reinforces our motivation for using bottom-up methods that combine up-to-date and periodic surveys with geo-spatial data to estimate population density, without having to rely on the census that is only available at most decennially in most countries, including Sri Lanka.

Coefficient of variation of the estimates

We calculated the coefficient of variation (CV) of predicted density estimates for each sub-district based on the bottom-up method, and compared them with the full-survey based density estimates for the entire country. The CV from the predicted density values is based on the average of a Poisson model, with errors clustered at the sub-district level. It was calculated using the Poisson command in Stata followed by the margins, pred(n) noesample vce(unconditional) command. The average CV of population density calculated directly from the full survey, across 329 sub-districts, is 50.0. This is calculated using the Horvitz-Thompson approximation, which allows for positive correlation in population density across GN Divisions within sub-districts. The average CV of population density calculated from the model predictions, for the same 329 sub-districts, is 6.7. Therefore, in the full sample, moving from the survey to the model-based predictions generates an approximate 85 percent reduction in the size of the CV of sub-district population density.