Evaluating nighttime lights and population distribution as proxies for mapping anthropogenic CO₂ emission in Vietnam, Cambodia and Laos

The ODIAC is a global high resolution (1 × 1 km) fossil fuel CO₂ emission data product (Oda and Maksyutov 2011, Oda et al 2018). The ODIAC is based on spatial disaggregation of CO₂ emission estimates made by the Carbon Dioxide Information Analysis Center (CDIAC) at the Oak Ridge National Lab (ORNL) (Boden et al 2016). CDIAC emissions are estimated by fuel type (solid, gas, and liquid fuels, bunker fuel, and gas flares) plus cement production, rather than the emission sector that is often used for the national inventory compilation (Marland and Rotty 1984). The ODIAC spatial disaggregation is done in two steps. First, emissions from point sources (mainly power plants) are estimated and mapped using the power plant emission estimates and geolocation taken from a global power plant database. The rest of the emissions (country total minus point source emissions), which we refer to a non-point source emissions, are distributed using the spatial distribution of satellite-observed nightlights (NTL) intensities (Oda and Maksyutov 2011, Oda et al 2018). Non-point source emissions are disaggregated to a 1 km × 1 km spatial resolution using Defense Meteorological Satellite Program (DMSP) calibrated radiance NTL data, with mitigated saturation effect, developed by NOAA's Earth Observation Group (Oda et al 2010). The calibrated radiance NTL data are a merged product of the regular DMSP NTL product and benefits from reduced gain observations (Ziskin et al 2012). Oda et al (2010) show an improved spatial emissions distribution from the original publication by Oda and Maksyutov (2011) due to the use of the calibrated radiance data.

Globally, this emission disaggregation is done for 65 individual countries and 5 aggregated geographical region groups (Oda and Maksyutov 2011). Vietnam, Cambodia and Laos are a part of the Asia Pacific geographical region group in the ODIAC country and region categorization. Since the country emissions for those three countries are once aggregated to the regional total before the emissions disaggregation, the country total emissions do not exactly match with the original CDIAC estimates. However, the geographic aggregation for Vietnam, Cambodia and Laos is necessary for an emission modeling framework that aims at quick, timely updates using the latest fuel statistical data from companies such as BP. Further details of the ODIAC approach and methodology are described elsewhere (Oda and Maksyutov 2011, Oda et al 2018). In this study, we excluded point source emissions and only use non-point (diffuse, area) source emissions (total minus - point source emissions) and specific to the study region of interest, focusing on the remotely-sensed NTL, nonpoint source ODIAC emissions. The version of ODIAC used in this study (ODIAC 2017) covers 2000 to 2016 but we use data for the years 2000, 2005 and 2010. The data product is available from http://db.cger.nies.go.jp/dataset/ODIAC/ (Oda and Maksyutov n.d.).

The gridded population data for this study is from the WorldPop Project (www.worldpop.org 2018). The underlying method is one of the more advanced techniques for gridding population distributions which relies on a hybrid approach of using a statistical weighting layer for estimating population density constrained by a dasymetric redistribution of the administrative unit population counts (Stevens et al 2015). Such an approach compares favorably or outperforms other techniques for producing gridded population maps (Sorichetta et al 2015, Stevens et al 2015, Reed et al 2018) while recognizing results are dependent on the spatial fidelity and accuracy of the ancillary data and census data used in the model.

The statistical component of the WorldPop model involves training a random forest (RF) model (Breiman 2001) to create a population density layer that will be used to disaggregate the administrative unit-based population counts to a regular grid of fixed spatial resolution. The RF-model is a non-parametric, non-linear statistical machine-learning approach that combines a set of decision trees into an 'ensemble' learner of multiple trees for a stronger output prediction. The ensemble of models is combined with random sampling of both training observations (bagging) for individual trees and covariate selection during tree growth. The individual trees and their prediction estimates can then be validated against withheld observations (out-of-bag errors, OOB), which can then be used across all trees in the RF to assess both variable importance and prediction performance (Breiman 2001, Liaw and Wiener 2002, Strobl et al 2009). The OOB error represents a robust and unbiased measurement of the prediction accuracy of the RF model (Breiman 2001) and a reliable proxy of the accuracy of the final gridded population datasets produced using the RF-based approach (Sorichetta et al 2015, Stevens et al 2015).

Predictions of population density from the RF model at the pixel level are inherently biased because the model is parameterized at the administrative unit level. Therefore, the pixel level estimates of population density are used as weights to anchor a dasymetric redistribution approach (Mennis 2003) of administrative census counts redistributed to a regular grid of fixed spatial resolution (100 m). Full details on the combined machine-learning and dasymetric modeling technique is found in (Gaughan et al 2015, Sorichetta et al 2015, Stevens et al 2015).

The WorldPop model is parameterized for each year of interest using an aggregated set of administrative units (Laos n = 10035, Cambodia n = 1621, Vietnam n = 688) and by maintaining the same consistent boundary definitions over time (Gaughan et al 2016). By relying on census data to inform the modeling process, the gridded population outputs effectively represent residential population counts which is in contrast to data such as Landscan which models ambient population counts (Dobson et al 2000). In addition, the 100 m gridded population datasets are eventually spatially aggregated (by summing population per grid cell) and coregistered to match the 1 km spatial resolution outputs from the ODIAC model. Population grids were done for 2000, 2005 and 2010. In terms of the RF model parameterization, subnational population counts for 2000, 2005 and 2010 were either extracted, interpolated, and extrapolated using two census dates (i.e., ${t}_{0}$ and ${t}_{1}$) for all three countries (i.e., 2005 and 2010 for Laos, 1998 and 2008 for Cambodia, 1999 and 2009 for Vietnam). This was done by calculating the growth rate of each administrative unit within each country and applying it to the corresponding population counts (Doxsey-Whitfield et al 2015). For each administrative unit, the corresponding exponential growth rate (r) was calculated using the following formula:

$$\begin{eqnarray}&&r=\displaystyle \frac{1}{t}\,\mathrm{ln}\left(\displaystyle \frac{{P}_{1}}{{P}_{0}}\right),\end{eqnarray} \tag{ 1 }$$

where r is the growth rate of a given administrative unit between ${t}_{0}$ and ${t}_{1},$ ${P}_{0}$ is its total population at time ${t}_{0},$ ${P}_{1}$ is its total population at time ${t}_{1},$ and t represents the number of years between ${t}_{0}$ and ${t}_{1}.$

Geospatial covariates used as input to the RF have been standardized, harmonized, and co-registered, to the CIESIN Gridded Population of the World v4 archive of administrative-boundaries (http://sedac.ciesin.columbia.edu/downloads/docs/gpw-v4/gpw-v4-country-level-summary-rev10.xlsx) across all three countries of interest. We followed the approach outlined in Gaughan et al (2016) and used both temporally-invariant and temporally-explicit covariates (table 1) while excluding NTL data. Temporally-invariant data include topography and slope, average annual precipitation and temperature (representative of the current conditions), presence of main roads and their intersections, waterways, water bodies, and coastlines. Temporally-explicit data, for 2000, 2005 and 2010, include the European Space Agency (ESA) Climate Change Initiative (CCI) land cover layers and presence of protected areas (WDPA 2017). For all temporally-invariant categorical data we produced the corresponding 'distance-to' covariate dataset, while for the temporally-explicit categorical data we calculated the distance-to-edge covariates, where distances inside the edge are negative and distances outside the edge are positive. See (Lloyd et al n.d.) for the production methodology for how the various covariate datasets were assembled and harmonized. Figure 1 shows the entire study region with land cover, roads, and rivers displayed along with panels of main city area gridded population patterns for 2010.

Zoom In Zoom Out Reset image size

To determine the agreement between remotely-sensed NTL and gridded population, we examine the relationship between non-point source CO₂ emissions and population at multiple administrative levels over the three time points. We then compare the NTL-disaggregated CO₂ emission estimates (i.e. the ODIAC model) to a disaggregation of the CO₂ emission estimates by the gridded population (i.e. WorldPop model) for each year. By comparing disaggregations based on NTL versus a simple, per capita disaggregation we can directly interrogate how disassociation between nightlights and population estimates might affect the utility of either approach. Informing the underlying gridded population model is a set of geospatial covariates (table 1) that provide some insight as to the disaggregation of population counts, over space and time, making the population-driven disaggregation approach useful to assess similarities and differences in how NTL, as a singular covariate, redistributes the CO₂ emission estimates.

To disaggregate by population, we summed the total population count to the national level and created a population proportion layer for each year (e.g. 2010 people per pixel/total population by country). Next, we multiplied the proportion of population through by the total CO₂ by country for a given year. We did this for the three years of interest (2000, 2005, and 2010) resulting CO₂ emissions disaggregated by the population model described in section 2.2. We present different metrics to highlight differences in NTL-driven versus population-driven model outputs and identify spatial distributions related to the underlying data informing respective CO₂ models for each year. In addition, as done by previous studies (Hogue et al 2016), we examine the agreement of the two CO₂ emission maps at different spatial resolutions.

The NTL model indicates expanding source areas for all three countries in a spatially-explicit manner from the 2000–2005–2010 period (figure 2). Vietnam, which is a more urbanized country, shows greater spatial distribution in the estimated CO₂ emissions as the number of lit grid cells in Vietnam are proportionally higher than in Laos and Cambodia.

Zoom In Zoom Out Reset image size

Similar to the NTL maps, the annual gridded population data maps show a general increase in population counts over time (figure 3). The percentage of variance explained by the RF model and the prediction error associated to it are equal to 79% and 0.72, respectively, for all three years. The relative importance of a covariate in the RF model is captured by percent increase in mean squared error (%MSE). This means, in the ensemble model prediction phase, that the lower the %MSE, the more important a given covariate will be relative to other covariates when randomly selected and included in a node prediction. Results from the model in this study show that distance to roads and rivers are important to the model in terms of random effects for all three years (figure 1). Agriculture is also important in 2000 and 2005 but becomes slightly less so in 2010 (Supplemental (S1)). Full plot overviews for the model are provided in S1.

Zoom In Zoom Out Reset image size

Figure 4 highlights some important differences in per pixel estimates of CO₂ emissions from the NTL-driven (i.e. ODIAC) model versus the POP-driven model. The most noticeable difference between figures 4(a) and (b), is that the population-driven model allocates estimates of CO₂ emission to every grid cell while NTL does not. The gridded population modeling approach may estimate fractional population in a given grid cell (Stevens et al 2015) and will not predict zero in grid cells unless the administrative unit contains zero people. The type of approach performs favorably, and in many cases outperforms, other types of gridded population modeling approaches (Reed et al 2018), which may translate into better gridded residential CO₂ products. The NTL-driven model includes zero as a possible grid cell value (i.e. grey-shaded pixels in figure 4) and thus CO₂ emissions are more concentrated into certain areas for the NTL-driven model rather than spread across the entire study area.

Zoom In Zoom Out Reset image size

When we difference the NTL-driven model estimates from the POP-driven model for disaggregating CO₂ emission estimates, the regional biases in the CO₂ emission from the NTL-driven model are more apparent. In figure 5, positive differences indicate grid cells where the NTL-driven model estimates higher emissions and negative differences indicate areas where population-based CO₂ emissions estimates are higher than NTL-based estimates. Thus, figure 5 illustrates that through time, nighttime lights-based emissions disaggregation tends to concentrate those estimates in different ways than population-based estimates. It is clear that population-driven model of emissions estimates, which excludes nighttime lights from ancillary data sources, tends to estimate higher emissions outside of the highly urbanized and developed areas of Vietnam, Cambodia, and Laos (shown by the orange-red tones in figure 5). These tend to be areas along major roads, intersections, and where population has been counted in census data but where nighttime lights may not be present. However, across 2000–2010 we see an expansion in the concentration of CO₂ source area estimates due to an increase in nighttime lighted areas most likely associated with development, as noted by the increase in 'purple' grid cells from 2000–2005–2010.

Zoom In Zoom Out Reset image size

Figure 6 shows the level of disagreement between the two gridded CO₂ datasets by spatially aggregating at multiple spatial resolutions (i.e. by coarsening the two 1 km gridded datasets to 2, 3, 5, 10, 20, 50 and 100 km). As done in Oda et al (2019), we calculated the sum of the absolute differences at grid levels and defined the initial difference at 1 km as 100%. The results provide a means to assess how quickly the various datasets for each country and year converge on similar patterns regardless of what data informs the CO₂ disaggregation. At the country level, while decreasing the spatial resolution of the two datasets especially at the first 25 km or so, there is a sharper gradient in the decrease of disagreement for Vietnam than for Laos and Cambodia. This suggests that, in general, there is a stronger correlation between the NTL and population spatial pattern in Vietnam than in the other two countries (with a better correlation in Laos than Cambodia), and that this correlation is stronger at lower spatial resolutions. In other words, the coarser the spatial resolution of the gridded CO₂ datasets, the more similar the outputs that are obtained by aggregating the NTL-based and the POP-based CO₂ grid.

Zoom In Zoom Out Reset image size

Temporally, by decreasing the spatial resolution of the two gridded CO₂ datasets for Vietnam, the disagreement is lowest for 2000 than for 2005 and 2010. Quite the opposite is observed for Laos and Cambodia, with the disagreement between the two corresponding gridded CO₂ datasets generally lower for 2010 than for 2000 and 2005. In other words, while for Vietnam the outputs obtained by aggregating the NTL-based and POP-based CO₂ grids decrease 'faster' overall compared to Cambodia and Laos, the pattern across the different years suggests Vietnam has a different pattern of development and land use, which is reflected by the associations of NTL to population across time. We also note that the differences did not converge to zero. At a coarser resolution, the sum of differences is not largely impacted by the differences in small spatial patterns. We interpret the remaining differences at coarse spatial scale as a proxy bias due to the use of the NTL.

To examine different spatial aggregations of the CO₂ emission estimates for 2000, 2005, and 2010, we sampled the NTL-based ODIAC and WorldPop population distribution outputs using different administrative unit level boundaries (i.e., n = 11,163 (level 3), n = 678 (level 2), and n = 63 (level 1) administrative units) and show the correlation coefficients for each level, by year, in table 2. The results highlight expected patterns for Vietnam, with improved correlation over time and improved correlation for coarser administrative unit levels (1st level is better than 2nd level, and better than 3rd level). This is consistent with results shown in figure 6 noting less disagreement when data sets are spatially coarsened. Laos and Cambodia show some deviation from this rule which we discuss further below. In addition, we plot the percent shares in the county total population to the NTL-based emissions from the ODIAC model which is located in the S1.

We also examined the numbers of the administrative units that NTL failed to allocate CO₂ emissions for (administrative units with zero emissions) (S1). At the finest administrative level for each country, the numbers of the administrative units consistently decreased especially from 2005 to 2010. For example, 2,847 administrative units in Vietnam in 2000, which accounted for 54% of the total source area (approximated by the sum of admin unit areas) and housed 15% of the total population, were not identified by NTL (zero emissions allocated). In other words, using NTL as the only proxy data for disaggregation failed to identify those 15% of the population and the their CO₂ emissions were misallocated somewhere else. However, we believe these emissions might not be significant as this 15% people's activity do not seem to be CO₂ intensive (no lights, or not developed). In 2010, 31% of the source region remains zero emissions, but that only accounts for 6% of the population. The numbers gradually decreased to 42% in 2005 and then 31% in 2010. Cambodia showed a drastic change from 92% in 2000, 90% in 2005 and then 46% in 2010. Laos also showed the change in the similar way (62% in 2000, 63% in 2005 and 30% in 2010). This analysis is meaningful as this is not impacted by the differences in spatial patterns and population is constrained at these administrative levels. We present this as a loose estimate of error associated with NTL-based emission downscaling at the administrative unit level.