Mapping Brick Kilns to Support Environmental Impact Studies around Delhi using Sentinel-2

: Cities lying in the Indo-Gangetic plains of South Asia have the world’s worst anthropogenic air pollution, which is often attributed to urban growth. Brick kilns, facilities for producing ﬁred clay-bricks for construction are often found at peri-urban region of South Asian cities. Although brick kilns are signiﬁcant air pollutant emitters, their contribution in under-represented in air pollution emission inventories due to unavailability of their distribution. This research overcomes this gap by proposing publicly available remote sensing dataset based approach for mapping brick-kiln locations using object detection and pixel classiﬁcation. As brick kiln locations are not permanent, an open-dataset based methodology is advantageous for periodically updating their locations. Brick kilns similar to Bull Trench Kilns were identiﬁed using the Sentinel-2 imagery around the state of Delhi in India. The unique geometric and spectral features of brick kilns distinguish them from other classes such as built-up, vegetation and fallow-land even in coarse resolution imagery. For object detection, transfer learning was used to overcome the requirement of huge training datasets, while for pixel-classiﬁcation random forest algorithm was used. The method achieved a recall of 0.72, precision of 0.99 and F1 score of 0.83. Overall 1564 kilns were detected, which are substantially higher than what was reported in an earlier study over the same region. We ﬁnd that brick kilns are located outside urban areas in proximity to outwardly expanding built-up areas and tall built structures. Duration of brick kiln operation was also estimated by analyzing the time-series of normalized difference vegetation index (NDVI) over the brick kiln locations. The brick kiln locations can be further used for updating land-use emission inventories to assess particulate matter and black carbon emissions.


Introduction
Background History of making fired clay-bricks in the Indian subcontinent is almost 4000 years old. Several brick structures have been excavated in the Indus Valley civilization archaeological sites where bricks were used for granaries, sewer system and other historic civic infrastructures. In most South Asia even today bricks are used as 'building blocks' in construction. In contrast to the reinforced inventory for Asia (REAS v2, [30]) due to lack of sufficient information. The latest version of REAS v3 [31] employs statistical estimates of brick production to infer their emission. Such assumptions in the emission inventory lead to uncertainties in simulating pollutant concentrations when using chemical transport models. For example, Zhong et al. [32] found simulated elemental carbon and SO 2 concentrations in the Kathmandu valley were underestimated compared to observations due to under-representation of brick kilns in the commonly used global emission inventories. A brick kiln location inventory is thus, beneficial for reducing emission uncertainties in South Asia and leading to accountability of their environmental impacts.
The research related to mapping brick locations itself is limited and recent [33]. In high spatial resolution (<3 m) satellite imagery brick kilns can be identified through three features: brick kiln chimney, bricks pile and slabs pile [20,33]. Furthermore, in nadir orthometric view these features are laid out in a ellipsoidal shape (or a circular shape depending on the type of brick kiln). Shadow of a brick kiln's chimney for imagery acquired at low sun elevation can be a useful visual indicator of a kiln in visible satellite bands. These features are shown in Figure 1d. Another reason that prevents identification of brick kiln through remote sensing is the prohibitive monetary cost associated with procuring high resolution imagery over a large area, specially where brick kilns are spread sparse. So far, Foody et al. [28] and Nazir et al. [34] have developed deep learning based methodology for brick kiln identification by using high resolution images available in Google Earth. However Google Earth's high resolution imagery suffers from the drawback that its satellite imagery is infrequently and non-uniformly updated, spanning intervals of several months over megacities to several years over rural areas. This is consequential for brick kiln detection due to the semi-permanent nature of brick kiln operations. New brick kilns are built to meet the brick demand or old ones are dismantled if fine soil becomes unavailable. Brick kiln operations are dynamic, responding to demand of bricks as well as availability of fine alluvial soil.
Thus a key challenge remains how can brick kiln be monitored regularly. Specifically, is it possible to use frequently updated publicly available dataset for detecting brick kilns. The objective of this study is investigate the use of Sentinel-2 imagery for detecting brick-kilns around Delhi. The originality of this paper is to demonstrate the use of open imagery for mapping brick kilns.

Methodology
The flowchart for mapping brick kiln locations followed in this research is shown in Figure 2.

Data
Sentinel-2: Surface reflectance of European Space Agency's Sentinel-2 Level 2 bands (B4: red, B3: green, B2: blue, B8: near-infrared, B11: shortwave infrared-1 and B12: shortwave infrared-2) was used. Sentinel-2 images acquired between 1 January 2018 and 1 June 2018 with cloud pixel percentage lower than 20% were cloud-masked and median-composited. The period of January to June was chosen for compositing images as brick kilns operate during these dry months. Due to the operation and preparation of fired red bricks, brick kilns visibly stand apart from the surrounding land-cover. In other months, due to seasonal suspension of operation, kilns are sometimes covered under small grass. The visible bands (red, green and blue) each have 10 m resolution making them suitable for delineating the shape of brick kilns from space and were used for object detection. The other bands (shortwave infrared-1 and shortwave infrared-2) have coarser resolution of 20 m and were used (along with visible bands and near-infrared band) for pixel classification. The images were composited over Google Earth Engine (GEE) [35] platform, where atmospheric and geometrically corrected Sentinel-2 dataset is available, and downloaded for further processing.
Landsat 7, Landsat 8: Since most of the brick kilns locations were originally agricultural fields prior to their conversion to brick kiln, a time-series of normalized difference vegetation index (NDVI = (nearinfrared −red/nearinfrared + red)) over brick kiln location can indicate when the agricultural field was transformed to brick kiln. Due to availability of a long time-series of USGS's Landsat 7 imagery since 1999, they were used for estimating NDVI at brick kiln locations (further discussed in Section 3.3) after atmospheric correction and cloud-masking. To discuss the relationship of brick kilns with urban growth of built-up areas in and around Delhi, imagery from 2001 and 2017 was classified using surface reflectance from Landsat 7 and Landsat 8 respectively. These datasets were prepared over Google Earth Engine (GEE) [35] platform.
Vertical construction: Digital building height (DBH) data was prepared by extracting the normalized digital surface model from publicly available Japan Aerospace Exploration Agency's 30 m resolution ALOS World 3D digital surface model (AW3D30) [36]. Misra et al. [37] have already demonstrated the possibility of extracting digital building heights with a root mean square error of 6.92 m. The DBH represents the building heights for structures built upto 2010, since the AW3D30 digital surface model is a composite of the images acquired between 2005 to 2010.
GIS thematic layers: Road and river network vector shapefiles were obtained from Hijmans [38] to explore drivers for the location of brick kilns.

Location
The location for this research was the region peripheral to the Delhi state in Indo-Gangetic plains of India. This region lies in the floodplains of the river Yamuna and its tributary river Hindon. The floodplains are covered with rich alluvial top soil. The land-use of the peri-urban area varies from sub-urban to rural settlement and agricultural lands. The last two decades have seen rapid construction activity in the region surrounding Delhi and its satellite cities. This makes it a suitable region for setting up of brick kilns due to abundance of raw material (fine top soil) and close proximity to the market (building construction sites). Additionally this region experiences severe air pollution episodes during the dry winter season coinciding with the period of brick kiln operations. Fine particulate matter (PM 2.5 ) concentrations rise to over 250 µg/m 3 , more than 10 times the World Health Organization specified daily PM 2.5 concentration guideline [39]. Since an earlier research has also explored manual count of brick kilns for preparing emission inventory around Delhi [13] as part of its emission inventory preparation effort, selection of the same region allowed for validation of our brick kiln detections.

Transfer Learning Based Object Identification
Deep learning based object detection approaches have risen in popularity due to significantly better performance than conventional approaches such as hierarchical segmentation or multi-resolution segmentation. However deep learning methods require huge training dataset to learn weights and biases for the artificial neurons. Concurrently, lack of a sufficient annotated training data is a common concern for application of deep learning based object detection to remote sensing. Transfer learning is a promising approach for those applications which are constrained by small training dataset. Transfer learning facilitates application of weights learned on a general domain dataset to be fine-tuned for a specific domain dataset. Domain adaptation involves two steps, (1) freezing the layers of a pre-trained model (trained on a general domain) and adding new trainable layers on top of the frozen layers and training them over the target domain dataset; (2) partially unfreezing the higher layers and fine-tuning the unfrozen layers by training them at a very low learning rate. Here freezing a layer refers to not updating the layer parameters while training. Imagenet [40] is a repository of over 14 million natural color images with more than 20,000 human labeled categories. Several convolutional neural network architectures pre-trained on the Imagenet are available such as VGG-16 [41], ResNet-50 [42] and others. Due to the ability of transferring weights learned from a general repository like Imagenet, transfer learning is becoming popular within remote sensing space for object detection tasks such as ship detection in synthetic aperture radar imagery [43,44], airplane detection in very high resolution imagery [45] and slum detection in medium resolution imagery [46,47].
Training data preparation: For the present study, first the training dataset was prepared over which transfer learning was performed. Brick kilns locations were identified visually with Google Earth Pro desktop software, to prepare an inventory of the locations. Unaided identification of brick kilns in the Sentinel-2 imagery is difficult due to the relative coarse pixel size with respect to the size of the kilns. Care was taken to ensure that a co-located brick kiln is humanly identifiable in both the Google Earth as well as the Sentinel-2 image. For example, as can be seen in Figure 3b, a single brick kiln can have 4 to 10 Sentinel-2 pixels along its one dimension. Brick kiln dimensions also vary at different locations. After determining the locations of the brick kilns, 200 64 × 64 pixel-sized RGB (natural color composite) training patches from Sentinel-2 were extracted for each of the 4 classes-brick kilns, fallow land, built-up and vegetation. 80 additional patches for each class were prepared for validation and testing. The patches were selected in the RGB bands for transfer learning as these bands have the highest resolution among other Sentinel-2 bands and because the existing models are pre-trained on the Imagenet which only has RGB images. Band reflectance of the patches was scaled from 0 to 1 using a 'min-max' scaler. Mean of 1 and 99 percentile reflectance of the brick kiln patches was used for scaling the respective bands (red: 0.227 to 0.074 sr −1 , green: 0.187 to 0.083 sr −1 and blue: 0.148 to 0.059 sr −1 ) of the training dataset. Examples of scaled patches for the 4 classes are shown in Figure 4. To enhance training sample size data and prevent overfitting, an augmentation procedure was performed randomly on 40% training patches. Data augmentation included following transformations: translation, rotation by 50 o and vertical and horizontal flipping.  (a) brick kiln Transfer learning and prediction: Keras Python library [48] with Tensorflow [49] backend was used for transfer learning. VGG-16 network was used due to the ease of training the model (due to less number of hyperparameters) and early convergence of training. VGG-16 has 16 layers in total, 12 convolutional layers with 3 × 3 kernel-sized filters and 4 four fully-connected layers. VGG-16 weights pre-trained on the Imagenet dataset is also available as part of the Keras package. The network takes in a 224 × 224 RGB patch as input and outputs the probability of that patch belonging to each of the 1000 classes on which the VGG-16 network was trained. Transfer learning of brick kilns patches using the pre-trained VGG-16 network was performed in the two steps mentioned earlier: adding layers to VGG-16 and training added layers, and fine-tuning the network. Firstly, the output layer of the VGG-16 network was replaced by two successive fully connected layers with the penultimate layer having 512 nodes with a ReLU (rectified linear) activation and the final layer having 4 nodes with a 'softmax' activation. The 4 nodes of the final layer output the classification probability of the target 4 classes (that is, brick kiln, fallow-land, built-up, and vegetation). The modified architecture is shown in the Table S1. The weights and bias of the fully connected layers were trained for 300 epochs with a batch size of 64 with the 'Adam' (Adaptive Moment Estimation) optimizer using a moderate learning rate (0.001) and 'categorical cross entropy' loss function. The goal of the optimizer is to minimize the loss function. Categorical cross entropy loss is used for multi-class classification tasks. Secondly, fine tuning was performed wherein the initial 12 convolutional layers were frozen. Unfreezing more initial layers always runs the risk of overfitting. This is because the initial layers learn general features (such as edges or lines) while the layers near the output learn features specialized for the target domain. The weights of the unfrozen layers were fine-tuned for 1000 epochs with stochastic gradient descent (SGD) optimizer using a smaller learning rate (0.0001) with momentum (0.9) and a 'categorical cross entropy' loss function. Providing a momentum prevents the SGD optimizer from getting stuck in local minima despite a low learning rate. Finally, by giving a RGB patch as an input, this network outputs the probabilities of that patch belonging to each of the classes, brick kiln, fallow-land, built-up, and vegetation.
As VGG-16 itself is not a localization algorithm, brick kiln location within the patch was further identified using a sliding-windows approach. The 64 × 64 pixel sized patch window was slid with a 75% overlap with its neighbors in each directions. Thus any given patch partially overlapped with 15 other patches (excluding itself). In presence of a ground truth brick kiln, all the partially overlapping patches containing that brick kiln must identify the brick kiln with the highest probability. Consequently when at least 5 partially overlapping patches detected brick kiln with the highest probability, the overlapping region was assigned to contain a brick kiln.

Random Forest Based Pixel Classification
Although deep learning methods are quite adept at object identification, the relative coarse resolution of Sentinel-2 imagery with respect to the brick kiln features (such as chimney or presence of bricks) gives rise to several false positive identification. False positives occur due to those land-use objects which share geometric features similar to brick kilns in coarse resolution imagery such as red rooftops of round-shaped buildings or playgrounds. To remove such false positive detections, pixel-wise spectral classification was performed for the same target area using random forest classification. Random forest algorithm was chosen as it is non-parametric, not prone to over-fitting and good at dealing with any outliers in training data [50]. Sentinel-2 bands (red, green, blue, near-infrared, shortwave infrared-1, and shortwave infrared-2) and derived NDVI was used for random forest based classification (with 100 decision trees) to the class brick kiln soil, fallow-land, built-up, vegetation and water.
If a region that was proposed by the transfer learning algorithm as containing a brick kiln also contained at least 5 pixels classified as brick kilns by random forest classifier, then that region was accepted as the brick kiln location. Regions of brick kiln locations overlapping more than 70% were removed during post-processing through 'faster non-max suppression' algorithm [51].

Brick Kiln Identification
Initial transfer learning was able to achieve a training and validation accuracy as 65.42% and 66.33%, and the categorical cross entropy loss as 0.84 and 0.86. After fine tuning, there was further improvement to training and validation accuracy as 85.12% and 86.67%, and the categorical cross entropy loss as 0.36 and 0.52. The fine tuning learning is shown in Figure 5. This was further enhanced by taking intersections of the identified regions with the spectrally classified pixels. Figure 6 shows the brick kilns identified by the current methodology in the north-eastern portion of Delhi. It can also be seen that patches of bare-land are located near brick kilns. These patches are used for drying molded clay bricks before firing them in the kilns. The final confusion matrix is shown in Table 1. Precision (true positive/(true positive + false positive)) and Recall (true positive/(true positive + false negative)) of brick kiln detections in the north-eastern portion of the study domain was 0.99 and 0.72 respectively. This signifies that underestimation of brick kiln counts by about 28%. In all 1564 brick kilns were detected around Delhi. If we account for 28% underestimation, the true count of brick kilns would be about 2172 brick kilns.
The overall F1 score (2 × (Recall × Precision)/(Recall + Precision)) of our hybrid object detection and pixel classification approach is 0.83. This score is considerably lower than that obtained by using very high resolution imagery. For example using Google Earth imagery for brick kiln detection Nazir et al. [34] and Foody et al. [28] obtained F1 score of 1 and 0.94 respectively. Possible reasons for a lower F1 score for detections using a coarser resolution Sentinel-2 imagery were explored. Brick kilns operate in the dry season from November to July as mentioned earlier.
Brick kilns are detected with high accuracy if the kilns are surrounded by vegetation or soil such that the geometric shape of the kiln is discernible. However accuracy of the detection is low when the soil is spectrally similar to the brick kilns. This is true for kilns operating in arid soils lacking any vegetation such that the kilns are not discernable to the object detection algorithm. An example of one such false negative is shown in Figure 7 where the brick kilns could not be identified in the Sentinel-2 imagery. On the other hand, some temporarily inactive brick kilns had grass growing on them, leading to difference in spectral properties compared to other active brick kilns and causing their pixels not being classified as brick kilns. Brick kilns that were located completely adjacent to each other in the ground truth, were not detected as it was impossible to segregate or identify them as two distinct brick kilns in the Sentinel-2 imagery. Our count of 1564 kilns (as shown in Figure 8a) is almost 1.5 times the number of brick kilns discussed earlier by Guttikunda and Calori [13] over the same area in 2013. Although their data collection was performed 5 years prior to the Sentinel-2 data used in the present study, several brick kilns locations in the south-western portion of the study domain (around the city of Gurgaon as shown in Figure 8c) were missed by their manual tagging approach. This suggests that the calculated contribution of brick kiln emission to Delhi's urban air pollution (PM 2.5 ) was underestimated by the authors of [13].

Drivers of Brick Kiln Locations
Brick kilns around Delhi are located in clusters situated between urban agglomeration of satellite towns (Figure 8a). The brick kilns are located along the flood-plains (river Yamuna) and along secondary roads (Figure 8b). Flood-plains provide access to the raw materials such as fine alluvial clayey soil and water, for producing bricks. Network of secondary roads is used for transporting bricks, possibly due to the current guideline of the Central Pollution Control Board (CPCB) that bars brick kilns from being located next to highways. Construction of buildings requires bricks, and brick kilns are located close to wherever dense horizontal expansion has taken place (Figure 8c). Brick kilns exist at the built-up edge of city or in the rural-urban zone where urban growth is taking place. Similarly, brick kilns are also located less than 10 km from multi-storey residential buildings (Figure 8d).

Age of Brick Kiln Operations
Based on field surveys, it was noticed that active brick kiln sites were covered with soil and red dust whereas abandoned brick kilns were covered under vegetation of grass and shrubs. This can provide additional information about a brick-kiln's operational status from satellite. Figure 9 shows an abandoned brick kiln and newly constructed brick kilns between 2006 and 2018 in northern Delhi. Time series of NDVI, an index of vegetation greenness, was calculated using cloud-free Landsat 7 imagery over brick kilns to determine their period of activity. NDVI varies between −1 to 1, with agricultural vegetation showing NDVI higher than 0.4 while sparse vegetation has NDVI between 0.2 to 0.4. NDVI of bare soil is lower than 0.2. NDVI over sample brick kiln locations are shown in Figure 10. An agriculture land is identifiable by its crop calendar of seasonal planting and harvesting of the crops. This is reflected as periodic variation of NDVI values. On the other hand if a low NDVI persists for a duration longer than the usual crop calendar, then that suggests the onset of brick kiln construction at that location. Figure 10a shows NDVI over a brick kiln whose operations have continued since prior to 1999 until present. Figure 10b shows NDVI over a crop field that was converted to brick kiln in 2008. Some brick kilns were also abandoned as shown by Figure 10c, where the location was undergoing agricultural activity until 2008, followed by brick kiln construction and then possible abandonment in 2017. The sharp increase of NDVI post 2017 suggests that brick kiln operation was suspended at that location.

Future Work
Future directions involve preparing a comprehensive regional dataset of brick kiln locations as well as their operational characteristics. Our methodology demonstrates the utility of publicly available Sentinel-2 imagery for brick kiln detection. This shall be further useful for mapping brick kilns over the whole South Asian region. Such a dataset would be useful for assessing drivers and impacts of brick kiln in different cities. As demonstrated in this paper, sudden and persisting shift in NDVI could indicate the duration of brick kiln operation. This could be further employed for assessing the net historical coal consumption and emission of trace gases from each brick kiln. However brick kilns do not operate all year round and not all brick kiln start or stop their operations at the same time. To accurately assess the contribution of brick kiln emission to air quality, exact operation months are needed. Such brick kiln dynamics cannot be inferred from NDVI trends. Availability of high spatio-temporal resolution imagery, for example Planet's daily 3 m spatial resolution imagery, can help in identifying the operation duration each year by tracking presence or absence of molded clay bricks outside the brick kiln (shown in Figure S1). Such an information will be useful for inferring activity status of individual brick kilns. Recent guidelines from India's CPCB have instructed kiln operators to upgrade to the relatively less polluting 'zig-zag' type brick kilns, having a rectangular shape and a taller chimney stack height for better dispersal of pollutant plume. The methodology proposed in the present work cannot distinguish between kiln types due to the limited resolution of the Sentinel-2 imagery. However changes in time series of surface backscatter from synthetic aperture radar (SAR) imagery, such as Sentinel-1, over the brick kiln locations may be used for inferring whether any upgradation of the brick kilns took place. Thus, information of brick kiln location identified in Sentinel-2 optical imagery can be combined with their SAR backscatter to assess which operators are complying with emission guidelines.

Conclusions
Brick kilns are known to be major emitters of black carbon, sulfur dioxide and other fine particulate matter, yet their spatial distribution remains largely unknown. We successfully demonstrated the first use of publicly available dataset, Sentinel 2, with deep learning methodology to detect brick kiln locations on a large spatial scale with a recall of 0.72, precision of 0.99 and F1 score of 0.83. 1564 brick kilns were detected around Delhi, higher than previously reported. Brick-kilns lie at urban-rural interface and are usually less than 10 km away from urban boundaries. The location of brick kiln is dependent on resource availability, such as fine clay from river-beds and access to road network. Brick kiln clusters were located near urban areas with tall buildings or areas where outward urban expansion has take place in last 16 years. Brick kilns are non-permanent structures whose duration of operation were estimated by a time-series of NDVI values over their location. As several current emission inventories under-represent the contribution of brick kilns to air pollution in South Asian cities, this dataset and methodology will be useful to adequately reveal their impact on air quality.