## 2 Materials and Methods
### Data
#### 2.2.2 Remotely Sensed Data  
This study used NTL data extracted from Visible Infrared Imaging Radiometer Suite Day-Night Band (VIIRS-DNB) products and Landsat 7 and 8 ETM images from 2012 to 2018. The time series data were processed on Google Earth Engine (GEE, [[https://earthengine.google.com](https://earthengine.google.com)]([https://earthengine.google.com](https://earthengine.google.com)) accessed on 24 June 2021), a cloud-based computing platform that allows users to conduct geospatial data analysis by providing various computation, analyzation, and visualization. GEE becomes an effective and constructive tool with various satellite images and geospatial datasets [26; 27]. GEE demonstrated high potential in the large-scale temporal urbanization analysis process. [PERSON] et al. [27] utilized Landsat time series data and computation algorithms in GEE to map urban areas and achieved an accuracy of 0.8 to 0.9. In this study, we adopted a method from previous research with processes based on the GEE platform. These include (1) extracting light intensity and area of lit from the annual mean composite of monthly VIIRS-DNB, (2) computing built-up indices from the annual median composition of cloud-masked Landsat 7 and 8 Surface Reflectance (SR) datasets, and (3) classifying the built-up areas from the derived indices using adaptive thresholding method ([PERSON]'s threshold) [28]. The built-up area extraction accuracy of 2012 and 2016 was assessed using 500 random sampling points from the national land-use dataset and virtual ground truth from Landsat 7 and 8 images each year. The overall accuracy of the build-up areas for 2012 and 2016 were 78% and 85% with kappa index values at 74% and 82%, respectively. The Landsat 7 (2012) and Landsat 8 (2013-2018) dataset was requested from the GEE platform, filtering by year, processing cloud mask, and selecting band 1 to 7. Each year of Landsat data has been reduced by median (derived from a histogram) for pre-compute the Landsat indices. We then adopted the built-up area extraction method from the previous paper, using the Modified Built-up Index (MBUI), which employs [PERSON]'s thresholding technique, and integrates with the nighttime light data [29]. The MBUI is calculated by Equation (1). MBUI = BUI - MNDWI,  
where:  
BUI = NDBI - NDVI,  
NDBI = \(\frac{(SWIR-NIR)}{(SWIR-NIR)}\),  
NDVI = \(\frac{(NIR-Red)}{(NIR+Red)}\), and  
MNDWI = \(\frac{(Green-SWIR)}{(Green+SWIR)}\). Considering the histogram of MBUI that represents some bimodal pattern, we choose [PERSON]'s method to classify the built-up area. [PERSON]'s thresholding is a fast and straightforward method to determine an optimal threshold value from an image histogram. The optimal threshold is given when the between-class variance is maximized [30]. [PERSON]'s thresholdinghas been used in several pieces of research [31; 32; 33]. The sample of the result is presented in Figure 3c. We employed GEE algorithms to request and export the monthly VIIRS-DNB dataset (2012-2018), which removes cloud cover using the VIIRS cloud mask product (VCM) and correct for stray light [34]. We calculated the annual average light intensity for each year, except the year 2012, which we calculated using the average light intensity from six-monthly datasets (June to December 2012) based on data availability. Then, we adopted the zonal statistic tool in ArcGIS Pro to identify the sum of light intensity in each district polygon from light intensity images. Zonal Statistics is a tool to calculate each zone defined by a zone dataset, based on values from another value raster dataset. A single output value (the sum of light intensity) was computed for every zone (district). Furthermore, the lit area was extracted using the thresholding method to differentiate each district's dark and bright areas. The threshold value was calculated for each year from 2012 to 2018 from the annual mean composite of VIIRS dataset by applying [PERSON]'s method. The threshold of the lit area was found to be between 19.25 to 20.93 nW cm\({}^{-1}\)sr\({}^{-1}\).